Antisolar differential rotation of slowly rotating cool stars

TL;DR

This paper investigates how slow stellar rotation can lead to antisolar differential rotation, with angular momentum transported towards the poles, through new simulations and analysis of angular momentum transport mechanisms.

Contribution

It introduces a new explanation for antisolar rotation in slow rotators based on off-diagonal eddy viscosity and simulation results, advancing understanding of stellar rotation profiles.

Findings

Simulations show angular momentum transport towards poles in slow rotation.

Negative radial gradient of angular velocity causes antisolar rotation.

Correlation between radial and latitudinal velocities affects temperature profiles.

Abstract

Rotating stellar convection transports angular momentum towards the equator, generating the characteristic equatorial acceleration of the solar rotation while the radial flux of angular momentum is always inwards. New numerical box simulations for the meridional cross-correlation $\langle u_\theta u_\phi\rangle $, however, reveal the angular momentum transport towards the poles for slow rotation and towards the equator for fast rotation. The explanation is that for slow rotation a negative radial gradient of the angular velocity always appears, which in combination with a so-far neglected rotation-induced off-diagonal eddy viscosity term $\nu_\bot$ provides "antisolar rotation" laws with a decelerated equator. Similarly, the simulations provided positive values for the rotation-induced correlation $\langle u_r u_\theta\rangle $, which is relevant for the resulting latitudinal temperature profiles (cool or warm poles) for slow rotation and negative values for fast rotation. Observations of the differential rotation of slowly rotating stars will therefore lead to a better understanding of the actual stress-strain relation, the heat transport, and the underlying rotating convection.

Equations

$$\frac{\partial}{\partial t}(\rho r^{2}\sin^{2}\theta{\it\Omega})+\nabla\cdot\left\{\rho r^{2}\sin^{2}\theta{\it\Omega}{\mbox{\boldmath$U$}}+\rho r\sin\theta\langle u_{\phi}\mbox{\boldmath$u$}\rangle\right\}=0,$$

$$\begin{array}{rl}{Q}_{r\varphi }=-{\nu }_{\Vert }& \sin \theta \frac{r\partial \Omega }{\partial r},Q_{\theta \varphi }=-{\nu }_{\Vert }\sin \theta \frac{\partial \Omega }{\partial \theta },\end{array}$$

$$\displaystyle Q_{ij}=\langle u_{i}(\mbox{\boldmath$x$},t)u_{j}(\mbox{\boldmath$x$},t)\rangle,$$

$$\begin{array}{rl}{Q}_{r\varphi }=-{\nu }_{\Vert }& \sin \theta \frac{r\partial \Omega }{\partial r}+\\ & +{\nu }_{\perp }{\Omega }^2{\sin }^2\theta \cos \theta \frac{\partial \Omega }{\partial \theta }+{\nu }_{\Vert }V\sin \theta \Omega \end{array}$$

$$\begin{array}{rl}{Q}_{\theta \varphi }=-{\nu }_{\Vert }& \sin \theta \frac{\partial \Omega }{\partial \theta }+\\ & +{\nu }_{\perp }{\Omega }^2{\sin }^2\theta \cos \theta \frac{r\partial \Omega }{\partial r}+{\nu }_{\Vert }H\cos \theta \Omega \end{array}$$

$$\displaystyle{\cal N}=r\left(\begin{array}[]{lll}\sin\theta\ \nu_{\|}&-\cos\theta\sin^{2}\theta{\it\Omega}^{2}{{\nu_{\bot}}}&0\\ -\cos\theta\sin^{2}\theta{\it\Omega}^{2}\ {\nu_{\bot}}&\sin\theta\ \nu_{\|}&0\\ 0&0&0\end{array}\right)$$

$$V=\sum _{l=0}{V}^{(l)}{\sin }^{2l}\theta \sin \theta ,H=\sum _{l=1}{H}^{(l)}{\sin }^{2l}\theta \cos \theta$$

$$\frac{\mathrm{d}\log \Omega }{\mathrm{d}\log r}=V^{(0)}$$

$$\frac{\delta \Omega }{\Omega }=-\frac{1}{2}\sum _{l=1}\left(d{V}^{(l)}+\frac{H^{(l)}}{l}\right)$$

$$\displaystyle\rho\left(\frac{\partial\mbox{\boldmath$u$}}{\partial t}+(\mbox{\boldmath$u$}\cdot\nabla)\mbox{\boldmath$u$}\right)=-\nabla p+\nabla\cdot\tau+\rho\mbox{\boldmath$g$}-2\rho\mbox{\boldmath${\it\Omega}$}\times\mbox{\boldmath$u$},$$

$$\displaystyle{\partial\rho}/{\partial t}+\nabla\cdot(\rho\mbox{\boldmath$u$})=0,$$

$$\displaystyle\frac{\partial e}{\partial t}+\nabla\cdot\left(\left(e+p\right)\mbox{\boldmath$u$}\right)=\nabla\cdot\left(\mbox{\boldmath$u$}\cdot\tau-\mbox{\boldmath$F$}_{\rm cond}\right),$$

$$\displaystyle\tau_{ij}=\rho\nu\big{(}u_{i,j}+u_{j,i}-\frac{2}{3}(\nabla\cdot\mbox{\boldmath$u$})\delta_{ij}\big{)},$$

$$U=\frac{\mathcal{R}}{\gamma -1}\rho T$$

$$Q_{ij}=Q_{ij}^{(0)}-\epsilon (2{\Omega }^2{\delta }_{ij}-{\Omega }_i{\Omega }_j)$$

$$\begin{array}{rl}\epsilon & =\frac{2}{15}{\int }_0^{\infty }{\int }_{-\infty }^{\infty }\frac{{\nu }^2{k}^4-3{\omega }^2}{({\omega }^2+{\nu }^2{k}^4{)}^2}E\mathrm{d}k\mathrm{d}\omega \\ & =-\frac{2}{15}{\int }_0^{\infty }{\int }_{-\infty }^{\infty }\frac{\omega }{({\omega }^2+{\nu }^2{k}^4{)}^2}\frac{\partial }{\partial \omega }\left(({\omega }^2+{\nu }^2{k}^4)E\right)\mathrm{d}k\mathrm{d}\omega \end{array}$$

$$⟨u_r^2⟩=⟨u_r^{(0)2}⟩-\epsilon {\Omega }^2(2-{\cos }^2\theta ).$$

$$Q_{r\theta }=-\epsilon {\Omega }^2\sin \theta \cos \theta ,$$

$$Q_{xz}=-{\nu }_{\Vert }\frac{\text{d}{U}_z}{\text{d}x}+{\nu }_{\Vert }V\sin \theta \Omega$$

$$Q_{yz}={\nu }_{\perp }{\Omega }^2\cos \theta \sin \theta \frac{\text{d}{U}_z}{\text{d}x}+{\nu }_{\Vert }H\cos \theta \Omega .$$

$$Q_{xy}=-\epsilon {\Omega }^2\sin \theta \cos \theta -{\nu }_{\Vert }\frac{\text{d}{U}_y}{\text{d}x}.$$

$${\nu }_{\Vert }={\alpha }_{\mathrm{vis}}{\tau }_{\mathrm{corr}}{u}_{\mathrm{rms}}^2,$$

$${\nu }_{\Vert }=\frac{2}{15}{\int }_0^{\infty }{\int }_{-\infty }^{\infty }\frac{{\nu }^3{k}^6}{({\omega }^2+{\nu }^2{k}^4{)}^2}E\mathrm{d}k\mathrm{d}\omega$$

$${\nu }_{\perp }=\frac{48}{105}{\int }_0^{\infty }{\int }_{-\infty }^{\infty }\frac{\nu k^2{\omega }^2(5{\nu }^2{k}^4-3{\omega }^2)}{({\omega }^2+{\nu }^2{k}^4{)}^4}E\mathrm{d}k\mathrm{d}\omega$$

$$\begin{array}{rl}{\nu }_{\perp }=\frac{48}{105}{\int }_0^{\infty }{\int }_{-\infty }^{\infty }& \frac{\nu k^2{\omega }^2}{({\omega }^2+{\nu }^2{k}^4{)}^2}\\ & \left(\frac{E}{{\omega }^2+{\nu }^2{k}^4}-2\frac{\partial }{\partial \omega }\frac{\omega E}{{\omega }^2+{\nu }^2{k}^4}\right)\mathrm{d}k\mathrm{d}\omega .\end{array}$$

$${\int }_{-\infty }^{\infty }\frac{{\nu }^3{k}^6{\omega }^2(5{\nu }^2{k}^4-3{\omega }^2)}{({\omega }^2+{\nu }^2{k}^4{)}^4}\mathrm{d}\omega =\frac{\pi }{8},$$

Full text

1 Introduction

Fast stellar rotation together with turbulent convection leave an imprint on the stellar surface in the form of starspots (Solanki 2003; Strassmeier 2009). Nevertheless, numerous puzzles remain for the quantitative description of the concerted action of stellar rotation and magnetic-field amplification in cool late-type stars. Surface differential rotation in its solar and antisolar form is one of them. Direct numerical simulations and mean-field models deal with the impact of Reynolds stresses and thermal energy flows on angular momentum transport in rotating convection, which are thought to be responsible for the observed meridional flows and differential rotation (e.g., Kitchatinov & Rüdiger 1999; Käpylä et al. 2011; Warnecke et al. 2013; Gastine et al. 2014; Featherstone & Miesch 2015). Kitchatinov & Olemskoy (2011) showed that the meridional flow is distributed over the entire convection zone in slow rotators but retreats to the convection zone boundaries in rapid rotators. Mean-field models had already predicted antisolar differential rotation for stars with fast meridional circulation (Kitchatinov & Rüdiger 2004). Tracking sunspot and starspot migration from spatially resolved solar and stellar disk measurements (e.g., Künstler et al. 2015) provided us with direct observations of differential rotation also in its antisolar form (polar regions rotate faster than the equatorial regions) meaning that models can now be tested against observations.

Differential rotation from spot tracking is no longer confined to solar observations (for a summary of solar differential rotation tracing measurements see, e.g., Wöhl et al. 2010). Tracking starspot migration from spatially resolved stellar disk measurements (Doppler imaging) is a method to directly study stellar surface differential rotation. After pioneering work on image cross-correlation and smeared Doppler imaging by Donati & Collier Cameron (1997), Barnes et al. (2005) concluded that differential rotation decreases with effective temperature and rotation. We have now a number of (active) stars where differential rotation has been detected directly by means of Doppler imaging. A differential rotation versus rotation-period relationship from Doppler-imaging results was suggested by Kővári et al. (2017) in the form $\delta{\it\Omega}\simeq\pi/100$ rad/day, where $\delta{\it\Omega}$ is the equator-pole difference of the angular velocity. Such a very weak dependence of $\delta{\it\Omega}$ on the stellar rotation rate of dwarf stars and giants was first predicted by Kitchatinov & Rüdiger (1999).

Space-based ultra-high-precision time-series photometry allowed confirmation of the temperature dependency of surface differential rotation for stars over a wide range in the Hertzsprung-Russel diagram (Reinhold et al. 2013). The observed $\delta{\it\Omega}$ only varied from 0.079 rad/day for cool stars ($T_{\rm eff}=3500$ K) to 0.096 rad/day for $T_{\rm eff}=6000$ K, which is a rather weak variation of the observed differential rotation on the rotation periods. Further, the hotter stars show stronger differential rotation, peaking at the F stars where there is still a significant convective envelope but only comparably weak magnetic activity. On the contrary, despite their small differential rotation M stars appear to be highly dynamo-efficient (Gastine et al. 2013). The overall rate of rotation plays an observationally biasing role here because smaller stars rotate much faster than bigger ones. Comparison of the large set of differential rotation measurements from the Kepler mission with the theoretical predictions by the $\Lambda$ effect theory showed a fair agreement and gives us confidence in applying the method to a particular star. However, the photometric data do not contain information on the sign of differential rotation (but see Reinhold & Arlt (2015) for possible exceptions) and further spectroscopic time-series data are needed.

Benomar et al. (2018) report the asteroseismic detection of surface rotation laws of solar-type stars with rather large equator-pole differences of the angular velocity. Among the sample of 40 stars there are up to 10 candidates for antisolar differential rotation with a weak anticorrelation to rapid rotation. Asteroseismology for the two solar analogs 16 Cyg A and B (which rotate slightly faster than the Sun) provided positive equator-pole ${\it\Omega}$ differences only slightly larger than the solar value (Bazot et al. 2019). We are therefore tempted to study the possibility of antisolar rotation mainly for slowly rotating main sequence stars. Also, the results of numerical simulations by Gilman (1977) and later by Gastine et al. (2014); Brun & Palacios (2009); Guerrero et al. (2013); Gastine et al. (2014); Käpylä et al. (2014) suggest that stars with slow rotation possess antisolar rotation laws. Most recently, Viviani et al. (2018) reported 3D simulations of turbulent convection for rotation rates of the solar value and faster. A transition of antisolar to solar-like differential rotation happened for about 1.8 solar rotation rates (their Fig. 5). Simultaneously, the geometry of the dynamo-excited large-scale magnetic field became nonaxisymmetric. Even more important for the understanding of the differential rotation problem is that the radial rotation shear also simultaneously changed from negative (“subrotation”) to positive (“superrotation”). We demonstrate in the present paper that indeed the subrotation ${\it\Omega}$-profile generically belongs to the antisolar-rotation phenomenon of decelerated equators.

We present numerical 3D box-simulations of outer stellar convection zones subject to slow rotation with a fixed Prandtl number. We then check if the resulting azimuthal cross-correlations generate solar-type or antisolar-type rotation laws. Our analytical differential rotation model is described in Sect. 2 while the simulations of the rotating convection boxes are presented in Sect. 3. The results in terms of the eddy viscosity tensor for the angular momentum flow are given in Sect. 4 and they are finally discussed in Sect. 5.

2 Differential rotation

In the following our theory of differential rotation in shellular convection zones in the mean-field hydrodynamic approach is briefly reviewed. This is mainly the theory of angular momentum conservation including meridional flow and Reynolds stress, that is,

[TABLE]

where $\rho$ is the mass density, ${\it\Omega}$ the angular velocity, $U$ is the ensemble average of the fluid velocity, and $u$ are the fluctuating parts of the flow. Equation (1) describes the contributions of the two main transporters of angular momentum in (unmagnetized) rotating convection zones. The cross-correlations Q_{r\phi}=\langle u_{r}(\mbox{\boldmathx},t)u_{\phi}\rangle(\mbox{\boldmathx},t) and Q_{\theta\phi}=\langle u_{\theta}(\mbox{\boldmathx},t)u_{\phi}\rangle(\mbox{\boldmathx},t) describe the radial and latitudinal turbulent transport of angular momentum (Rüdiger 1989). In the simplest case they can be parametrized via the diffusion approximation,

[TABLE]

with $\nu_{\|}$ being the positive eddy viscosity111$\nu_{\|}=\nu_{1}\equiv\nu_{\rm T}$ in the notation of Kitchatinov et al. (1994). In this approximation the two cross-correlations would vanish for uniform rotation. If, therefore, under certain circumstances the cross-correlations for uniform rotation do not vanish, the Boussinesq formulation (2) can no longer be true and uniform rotation cannot form a solution of Eq. (1) in rotating turbulence fields. Indeed theory, simulation, and observation suggest that large-scale stellar convection produces finite values for the two mentioned cross-correlations, this phenomenon being referred to as the “$\Lambda$ effect”. For the Sun as a rapid rotator (compared with the typical correlation times) Hathaway et al. (2013) indeed reported positive latitudinal cross-correlations for the northern hemisphere and negative latitudinal cross-correlations for the southern hemisphere in contradiction to the simple diffusion approximation (2) which would provide opposite signs.

The symmetry properties of the cross-correlations $Q_{r\phi}$ and $Q_{\theta\phi}$ differ from those of all other components of the one-point correlation tensor,

[TABLE]

of a rotating turbulence field. While $Q_{r\phi}$ and $Q_{\theta\phi}$ are antisymmetric with respect to the transformation ${\it\Omega}\to-{\it\Omega}$, all other correlations are not. The turbulent angular momentum transport is thus odd in ${\it\Omega}$ while the other two tensor components – the cross-correlation $Q_{r\theta}$ included – are even in ${\it\Omega}$. It is easy to show that $Q_{r\phi}$ is symmetric with respect to the equator if the averaged flow is also symmetric. In this case, the component $Q_{\theta\phi}$ is antisymmetric with respect to the equator. These rules can be violated if, for example, a magnetic field exists whose amplitudes are different in the two hemispheres.

One can also show that isotropic turbulence even under the influence of rotation does not lead to finite values of $Q_{r\phi}$ and $Q_{\theta\phi}$. With a preferred (radial) direction $g$, a tensor $(\epsilon_{ikl}g_{j}+\epsilon_{jkl}g_{i})g_{k}{\it\Omega}_{l}$ linear in ${\it\Omega}$ can be formed that has nonvanishing $r\phi$ and $\theta\phi$ components. Rotating anisotropic turbulence is therefore able to transport angular momentum. The spherical coordinates $(r,\theta,\phi)$ are used in this paper if the global system is concerned while $(x,y,z)$ represent these coordinates in a Cartesian box geometry.

For the zonal fluxes of angular momentum we write

[TABLE]

for the radial flux and

[TABLE]

for the meridional flux. Here the first terms come from the Boussinesq diffusion approximation with $\nu_{\|}$ as the eddy viscosity while $V$ and $H$ form the components of the $\Lambda$ tensor describing the angular momentum transport of rigidly rotating anisotropic turbulence. The terms with ${{\nu_{\bot}}}$ follow from the fact that a viscosity tensor connects the Reynolds stress with the deformation tensor222${\nu_{\bot}}=\nu_{2}$ in the notation of Kitchatinov et al. (1994). The nondiffusive terms in the zonal fluxes (4) and (5) can be written by means of the stress-strain tensor relation $Q_{i\phi}=-{\cal N}_{ij}\nabla_{j}{\it\Omega}$ with

[TABLE]

(Rüdiger 1989). The positive standard eddy viscosity $\nu_{\|}$ is quenched by fast rotation. On the other hand, Kitchatinov et al. (1994) showed that for isotropic and homogeneous turbulence ${\nu_{\bot}}$ is positive. We use ${{\nu_{\bot}}}$ to refer to the rotation-induced off-diagonal viscosity term. The off-diagonal viscosity does not contribute at the poles or at the equator. We note that this term in Eq. (5) transforms a positive (negative) radial ${\it\Omega}$ gradient into positive (negative) cross-correlations which – as the solution of the equation for the angular momentum conservation – finally leads to accelerated (decelerated) equators. Hence, the rotation law of a convection zone can never be only radius-dependent. After the Taylor-Proudman theorem the isolines of the angular velocity ${\it\Omega}$ tend to become cylindrical so that (say) slower rotation in the depth of the convection zone is transformed to polar deceleration (solar-type rotation). If, on the other hand, the inner parts rotate faster than the outer parts then automatically the polar regions rotate faster than the more equatorial regions (“antisolar rotation”)333In the linear-in-${\it\Omega}$ approximation by Kippenhahn (1963) a very similar mechanism is realized via meridional flow..

The expansion

[TABLE]

is used for the normalized $\Lambda$ effect as in earlier papers. The coefficients $V^{(l)}$ and $H^{(l)}$ describe the latitudinal profile of the $\Lambda$ effect. Quasilinear theory of rapidly rotating anisotropic turbulent convection in the high-viscosity limit leads to $-V^{(0)}=V^{(1)}=H^{(1)}>0$ (with $V^{(l)}=H^{(l)}=0$ for $l>1$) which implies that $Q_{r\phi}$ vanishes at the equator. For rigid rotation and in cylindric coordinates $-H\sin\theta\cos\theta$ is the angular momentum flow in axial direction while the radial flux of angular momentum vanishes. The function $H=H({\it\Omega})$ is positive definite, meaning that the angular momentum is exclusively transported from the poles to equator parallel to the rotation axis. $V^{(0)}$ is always negative (Rüdiger et al. 2005b).

Either of the above-mentioned theories and simulations of the radial $\Lambda$ effect lead to results of the form $V\propto-\cos^{2}\theta$; that is, $V^{(0)}<0$ and $V=0$ at the equator. For slow rotation, $V^{(l)}$ and $H^{(l)}$ with $l>0$ become so small that a radial rotation law with

[TABLE]

results. Negative $V^{(0)}$ values generally lead to radial ${\it\Omega}$ profiles with negative shear. In this case a meridional circulation is driven by the centrifugal force which at the surface transports angular momentum towards the poles (“counterclockwise flow”). Hence, the equator rotates slower than the mid-latitudes which automatically leads to an antisolar rotation law with $\cos\theta\partial{\it\Omega}/\partial\theta<0$ at the surface.

If neglecting meridional flow and ${\nu_{\bot}}$, the Reynolds stress (10) maintains a latitude-dependent surface rotation law ${\it\Omega}={\it\Omega}(\theta)$ with

[TABLE]

for the pole-equator difference of ${\it\Omega}$ and with the normalized thickness $d$ of the convectively unstable layer with stress-free boundary conditions. Under the assumption that $V^{(l)}=H^{(l)}=0$ for $l>1,$ antisolar rotation would only be possible for (formally) negative $H^{(1)}$.

Figure 1 illustrates the consequences of (12). The equation of angular momentum is solved for a negative $V^{(0)}$. For this demonstration, meridional circulation and the off-diagonal viscosity ${\nu_{\bot}}$ have been neglected. The latitudinal $\Lambda$ effect represented by $H^{(1)}$ is varied from 1 to $-1$. Not surprisingly, a negative pole-equator difference of the surface rotation law (solar-type differential rotation) originates from $H^{(1)}=1$ (left panel). For $H^{(1)}=0,$ a shellular rotation profile results. Moreover, for $H^{(1)}=-1$ the solar-type surface rotation law changes to an antisolar-type surface rotation law with positive pole-equator difference (right panel). At the same time the isolines of the angular velocity of rotation change from disk-like (left panel) to cylinder-like (right panel). This phenomenon is due to Reynolds stress rather than to Taylor-Proudman theorem. In the left plot the angular momentum is transported by the $\Lambda$ effect along cylindrical planes to the equator causing the ${\it\Omega}$-isolines to become disk-like. In the middle plot the transport is radial, meaning that the ${\it\Omega}$-isolines become shellular and in the right plot the angular momentum is transported toward the rotation axis generating cylindrical ${\it\Omega}$-isolines.

The question remains as to how the meridional circulation would modifiy these results if it were included. The above finding that formally negative $H$ easily produces antisolar rotation profiles remains true if the meridional circulation due to radial shear is also taken into account. The vorticity of the circulation depends on the sign of the radial shear $\partial{\it\Omega}/\partial r$. At the surface it flows towards the equator for superrotation and towards the poles for subrotation (Kippenhahn 1963). On the other hand, a circulation which flows towards the equator at the surface of the convection zone (“clockwise flow”) produces differential rotation with an accelerated equator while it produces a polar vortex if it flows towards the poles (“counterclockwise flow”). One takes from Fig. 2 (top) that for all choices of $H^{(1)}$ the meridional circulation flows counterclockwise in the northern hemisphere, reducing the equatorial acceleration (left panel) or amplifying equatorial deceleration (right panel). Indeed, with circulation included, the bottom left plot of Fig. 2 with $H^{(1)}>0$ represents a model for convection zones with solar-type rotation laws, while the right panel of Fig. 2 with $H^{(1)}<0$ represents an antisolar rotation law. As the middle plots of Fig. 2 show, a meridional circulation towards the poles even produces a weak antisolar rotation without any $\Lambda$ effect. Typically, as a result of the Taylor-Proudman theorem the isolines of the angular velocity ${\it\Omega}$ (with meridional flow included) become cylindrical. This effect appears in all plots of the bottom row of Fig. 2 but it is most prominent for $H^{(1)}<0$ which even without circulation generates cylinder-like ${\it\Omega}$-isolines. Simultaneously, for rotation laws with small $\partial{\it\Omega}/\partial z$ the amplitude of the circulation sinks.

Several analytical studies of the $\Lambda$ effect led to positive $H^{(1)}$, that is, $\cos\theta Q_{\theta\phi}>0$, for rigid rotation. Also, numerical simulations of rotating convection (Hupfer et al. 2006) or driven anisotropic turbulence under the influence of solid-body rotation (Käpylä 2019a) provide transport of angular momentum towards the equator. Earlier, Chan (2001) found transport towards the equator only for fast rotation while for slow rotation occasionally the opposite result appeared. Simulations by Rüdiger et al. (2005b) of rotating turbulent convection with much higher resolution provided very small $Q_{\theta\phi}$ for slow rotation and large positive $Q_{\theta\phi}$ for fast rotation. There seemed to be no hope, therefore, of explaining antisolar rotation laws for rotating stars with a hydrodynamical theory of turbulent rotating flows. In this paper numerical simulations of rotating convection in boxes are presented providing transport of angular momentum towards the poles for slow rotation. This transport towards the poles, however, does not result from the $\Lambda$ effect but is due to the rotation-induced off-diagonal viscosity term ${\nu_{\bot}}$ in (5) in connection with a subrotation law $\partial{\it\Omega}/\partial r<0$ which appears for slow rotation (Viviani et al. 2018). For solar-like convection zones with rotation profiles quasi-uniform in radius and rotating with the present-day solar rotation rate the ${\nu_{\bot}}$ term does not play any role.

3 Rotating convection

We perform simulations for convection with a fixed ordinary Prandtl number $\mbox{Pr}=\nu/\chi$ with $\chi$ being the thermal diffusion coefficient. For stellar material the heat conductivity $\chi$ strongly exceeds the other diffusivities. Also the Roberts number ${\rm q}={\chi}/{\eta}$ with $\eta$ as the microscopic magnetic resistivity is therefore much larger than unity. For numerical reasons we must work with the approximate surface value $\mbox{Pr}=0.1$.

The simulations are done with the Nirvana code by Ziegler (2002), which uses a conservative finite volume scheme in Cartesian coordinates. The length scale is defined by the depth of the convectively unstable layer. Periodic boundary conditions are formulated in the horizontal plane. The upper and lower boundaries are impenetrable and stress-free. The initial state is convectively unstable in a layer that occupies half of the box. Convection sets in if the Rayleigh number exceeds its critical value. In the dimensionless units the size of the simulation box is $2\times 6\times 6$ in the $x$, $y$, and $z$ directions, respectively. The lower and upper boundaries of the unstably stratified layer are located at $x=0.8$ and $x=1.8$, respectively. The numerical resolution is $128\times 384\times 384$ grid points. The stratification of density, pressure, and temperature is piecewise polytropic, similar to that used in Rüdiger et al. (2012). The density varies by a factor five over the depth of the box, hence the density scale height is 1.2.

The code solves the momentum equation,

[TABLE]

where $\rho$ is the mass density, $u$ the gas velocity, $p$ the gas pressure, $g$ gravity, and ${\it\Omega}$ the rotation vector, in a corotating Cartesian box under mass conservation,

[TABLE]

together with the energy equation,

[TABLE]

where \mbox{\boldmathF}_{\rm cond}=-\kappa\nabla T being the conductive heat flux with the heat conduction coefficient $\kappa$. The viscosity tensor is

[TABLE]

and the total energy e=U+\rho\mbox{\boldmathu}^{2}/2 is the sum of the thermal and kinetic energy densities. An ideal gas with a constant mean molecular weight $\mu=1$ is considered, hence

[TABLE]

for the thermal energy density with $\cal R$ the gas constant and $\gamma=c_{\rm p}/c_{\rm v}=5/3$. The gas is kept at a fixed temperature at the bottom and a fixed heat flux at the top of the simulation box. More technical details including the boundary conditions have been described in Rüdiger et al. (2012).

Figure 3 gives the auto-correlations of the one-point correlation tensor (3) for convection that is subject to slow and fast rotation. The colatitude is $\theta=45^{\circ}$. As expected, the horizontal turbulent intensities are identical for slow rotation while the vertical intensity $\langle u_{r}^{2}\rangle$ has the dominating value. The latter is strongly suppressed by faster rotation (the radial turbulence intensity is reduced by more than a factor of two for ${\it\Omega}=30$) while there is almost no visible suppression of the horizontal components.

For driven turbulence in a quasilinear approximation we have

[TABLE]

with

[TABLE]

(Rüdiger 1989). The sign of the expressions follows from the form of the positive-definite spectrum $E(k,\omega)$ and is obviously negative for white-noise spectra but is positive for spectral functions which are sufficiently steep in $\omega$. Also for the maximally steep spectrum such as $\delta(\omega)$ the $\varepsilon$ is positive, describing a rotational quenching of the turbulence intensities.

The vector of rotation at the colatitude $\theta$ is \mbox{\boldmath{\it\Omega}}={\it\Omega}(\cos\theta,-\sin\theta,0), so that Eq. (18) gives

[TABLE]

The rotational quenching of the radial turbulence intensity shown in Fig. 3 can be described by Eq. (20) with $\varepsilon>0$.

Another direct consequence of (18) is the existence of the cross-correlation of radial and latitudinal fluctuations, that is,

[TABLE]

which vanishes at the poles and the equator by definition. In a sense, the correlation $Q_{r\theta}$ mimics the turbulent thermal conductivity tensor. If a radial temperature gradient exists, a negative cross-correlation $Q_{r\theta}$ organizes heat transport to the poles resulting in a meridional circulation towards the equator at the surface (Rüdiger et al. 2005a).

For positive $\varepsilon,$ the rotation-induced cross-correlation $Q_{r\theta}$ after (21) becomes negative. This theoretical result complies with the results of numerical simulations (Käpylä 2019a). The expression (18) only describes the rotational influence on isotropic turbulence. Rüdiger et al. (2005b) also considered the rotational influence on turbulence fields which are anisotropic in the radial direction with the general result that $Q_{r\theta}$ is always less than zero (northern hemisphere) for steep spectra $E$ and for all rotation rates.

Hereafter we switch to the Cartesian box coordinates $x$ (representing the radial coordinate $r$), $y$ (representing the colatitude $\theta$), and $z$ (representing the azimuth $\phi$) hence $Q_{r\theta}\to Q_{xy}$, $Q_{r\phi}\to Q_{xz}$ and $Q_{\theta\phi}\to Q_{yz}$. The shears $r\partial{\it\Omega}/\partial r$ and $\partial{\it\Omega}/\partial\theta$ translate into ${\textrm{d}}U_{z}/{\textrm{d}}x$ and ${\textrm{d}}U_{z}/{\textrm{d}}y$, respectively. After averaging over the horizontal ($yz$) plane the relations (4) and (5) turn into

[TABLE]

and

[TABLE]

The cross-correlation (21) completed by the viscosity term becomes

[TABLE]

The attention is focused here on the influence of the diffusion terms in Eqs. (23) and (24) in order to probe the existence of the viscosities $\nu_{\|}$ and ${\nu_{\bot}}$ by simulations. To this end the rotation rates are assumed to be so small that the nondiffusive terms in the relations (23) and (24) for the cross-correlations are negligible. The second term in Eq. (24) is positive for outwards decreasing meridional flow $U_{y}$, for example. Consequently, the cross-correlation $Q_{xy}$ should be positive for slow rotation and negative for rapid rotation, changing the sign at a certain value of the parameter ${\it\Omega}$ (which denotes the angular velocity ${\it\Omega}$ of the rotation in code units). One finds such a transition from positive to negative values for $5<{\it\Omega}<10$ in the simulations given in Fig. 4. The coincidence suggests that indeed the influence of viscosity terms in the expressions of cross-correlations may lead to direction reversals of transport terms as a function of rotation.

Without rotation, all cross-correlations vanish. For the slow-rotation models with ${\it\Omega}=1,$ already finite values appear (left plots in Figs. 4 and 5). At the radial boundaries the correlations $Q_{xy}$ and $Q_{xz}$ vanish by the boundary condition ($u_{x}=0$) but the horizontal cross-correlation $Q_{yz}$ remains finite; it is always positive at the top and bottom of the unstable box which indicates $H>0$ if a possible mean circulation $U_{z}$ were maximal or minimal at the top or bottom of the convection box (as it is, see Fig. 6).

The correlation $Q_{xy}$ for slow rotation is positive so that heat is transported towards the equator. At the same time the horizontal correlation $Q_{yz}$ assumes negative values. Independent of the rotation rates and for both hemispheres we find the general result that always $Q_{xy}Q_{yz}<0$. The simulations therefore show that angular momentum flux to the equator (poles) is always accompanied by heat transport to the poles (equator). Somewhere between ${\it\Omega}=5$ and ${\it\Omega}=10$ the cross-correlations $Q_{xy}$ and $Q_{yz}$ change their signs becoming negative ($Q_{xy}$) and positive ($Q_{yz}$) for fast rotation. These signs are well-known from the analytical expressions derived for driven turbulence for fast rotation. In this case, the angular momentum is transported inward as well as toward the equator by the convection; in other words, it flows along cylindric surfaces.

For increasingly fast rotation the amplitudes of the negative $Q_{xy}$ are increasing, contrary to $Q_{yz}$ which decrease. This is a basic difference for the two cross-correlations. We note that for the transformation ${\it\Omega}\to-{\it\Omega}$ the correlations $Q_{yz}$ change their sign which is not the case for $Q_{xy}$. The reason is that $Q_{xy}$ is even in the rotational rate ${\it\Omega}$ while the horizontal cross-correlation $Q_{yz}$ is odd.

4 The eddy viscosities

Equation (21) neglects the influence of a possible radial shear ${\textrm{d}}U_{y}/{\textrm{d}}x$ of a meridional flow. The question is whether large-scale mean flow characterizes the simulation box as in the simulations of Chan (2001) and Käpylä et al. (2004), which could be used to calculate the eddy viscosities after relations (23) and (24). The mean flows in the box have been calculated for slow rotation. The top panel of Fig. 6 gives flows in the meridional direction and the bottom panel gives zonal flows in the azimuthal direction. The basic rotation there has been varied from ${\it\Omega}=1$ to ${\it\Omega}=10$ and the Prandtl number is fixed at $\mbox{Pr}=0.1$. The following estimate concerns the first example with ${\it\Omega}=1$ with the cross-correlation $Q_{xy}\simeq 0.01u^{2}_{\rm rms}$ (from Fig. 5) and the shear $\delta U_{y}/\delta x\simeq-0.2$ (from Fig. 6). The standard eddy viscosity $\nu_{\|}$ in code units is 1.5 for slow rotation. The microscopic viscosity of the model in the same units is $\nu=6\cdot 10^{-3}$, meaning that for the Reynolds number $\nu_{\|}/\nu\simeq 240$. One finds very similar values for ${\it\Omega}=3$.

The dimensionless eddy viscosity $\alpha_{\rm vis}$, following

[TABLE]

may also be introduced which is often assumed in turbulence research to be $\alpha_{\rm vis}\simeq 0.3$. To find the correlation time $\tau_{\rm corr}$ an auto-correlation analysis as done by Küker & Rüdiger (2018) is necessary. The result is $\tau_{\rm corr}\simeq 0.1$ in code units, hence \alpha_{\rm vis}\lower 1.72218pt\hbox{;\buildrel<\over{\scriptstyle\sim};}0.5 which indeed is of the expected order of magnitude.

Figure 6 also demonstrates that the negative shear ${\textrm{d}}U_{y}/{\textrm{d}}x$ is always accompanied by a negative shear ${\textrm{d}}U_{z}/{\textrm{d}}x$ of the zonal flow, hence ${\textrm{d}}U_{y}/{\textrm{d}}x\cdot{\textrm{d}}U_{z}/{\textrm{d}}x>0$. Transformed to the global system, this means that a subrotating shell generates a counterclockwise circulation where the fluid drifts towards the poles at the surface. This type of flow pattern has already been described in the text below Eq. (11). On the other hand, the existence of $U_{z}$ allows us to estimate the off-diagonal viscosity after Eq. (23) with $H\approx 0$ for ${\it\Omega}=1$ or ${\it\Omega}=3$ to ${\nu_{\bot}}\simeq 0.05$, hence ${\nu_{\bot}}/\nu_{\|}\simeq 0.03$ in code units or ${\nu_{\bot}}/\nu_{\|}\simeq 3\tau^{2}_{\rm corr}$ in physical units.

Following Eqs. (22) and (23), the shear ${{\textrm{d}}U_{z}}/{{\textrm{d}}x}$ contributes to the cross-correlations $Q_{xz}$ and $Q_{yz}$. Because of its negativity, the viscosity term in (22) is positive, meaning that the actual value of $|V|$ is even larger than indicated by $Q_{xz}$ from the simulations.

4.1 Rotation-induced off-diagonal eddy viscosity

The same radial velocity gradient appears in the expression for the horizontal cross-correlation in a higher order of the rotation rate. The two viscosities in (23) can be expressed by the integrals

[TABLE]

and

[TABLE]

(Rüdiger 1989). Not surprisingly, the first integral is positive definite. On the other hand, the second integral is positive for all other monotonously decreasing spectra in line with the following expression:

[TABLE]

The integral in (27) is even positive for spectra of the white-noise-type. To demonstrate this point, we evaluate the frequency integral with uniform $E$. As

[TABLE]

here the ${\nu_{\bot}}$ is also positive. We note however that very steep spectra such as $\delta(\omega)$ lead to vanishing ${\nu_{\bot}}$ which explains the absence of antisolar rotation laws in the calculations based on that turbulence model (Kitchatinov et al. 1994). It therefore seems likely that observations of slowly rotating stars with a decelerated equator question the application of turbulences with $\delta$-like frequency spectra in stellar convection models. The steepest frequency spectra describe fluids in the high-viscosity limit while spectra with a quasi-white-noise behavior belong to the inviscid approximation.

As the corresponding integral for the standard eddy viscosity (26) is $\pi/2$ one obtains, for rather flat spectra, ${\nu_{\bot}}/\nu_{\|}\ \simeq(l^{2}_{\rm corr}/2\nu)^{2}\simeq\tau_{\rm corr}^{2}$. The latter relation requires that for the background viscosity $\nu$ the standard expression $\nu\simeq 0.5l^{2}_{\rm corr}/\tau_{\rm corr}$ be used. Because of ${\nu_{\bot}}>0$ and ${{\textrm{d}}U_{z}}/{{\textrm{d}}x}<0$ (see Fig. 6) negative contributions to the horizontal cross-correlation are produced. The existence of the off-diagonal viscosity ${\nu_{\bot}}$ therefore explains the resulting negativity of the cross-correlation $Q_{yz}$ for slow rotation (see Fig. 5). The result confirms the above finding that all values of $Q_{yz}$ are positive at the top and bottom boundaries where the mean shear vanishes. For faster rotation, the increasing positive values of $H$ overcompensate the negative contribution from the radial shear of $U_{z}$ which becomes increasingly unimportant (see Fig. 6). The transition from negative to positive $Q_{yz}$ happens in Fig. 5 for ${\it\Omega}\simeq 5$. Figure 6 also demonstrates that the shear ${\textrm{d}}U_{z}/{\textrm{d}}x$ grows for ${\it\Omega}<5$ and overcompensates the positive $H$ term producing the obtained negative cross-correlations.

As the mechanisms of the null crossings of $Q_{xy}$ and $Q_{yz}$ are different, the values for the critical ${\it\Omega}$ should not be identical for both cases.

4.2 Antisolar rotation

Models are considered of such slow rotation that $H\simeq 0$ and a (negative) $V^{(0)}$ exists in addition to the rotation-induced off-diagonal viscosity ${{\nu_{\bot}}}$. The existence of the latter has been indicated by the simulations for the horizontal cross-correlation $Q_{yz}$ shown in Fig. 5. This leads to functions of $H$ that are formally negative (for slow rotation) for which Fig. 1 demonstrates the appearance of antisolar rotation laws. The off-diagonal element ${\nu_{\bot}}$ is varied in Fig. 7 from ${{\nu_{\bot}}}=-0.5$ (left panel) through zero (middle panel) to ${{\nu_{\bot}}}=0.5$ (right panel). For the models shown in the top row of the plots the resulting meridional circulation is artificially suppressed. One finds that negative $V^{(0)}$ always produces rotation laws with negative radial shear, which in combination with positive ${{\nu_{\bot}}}$ leads to antisolar rotation and vice versa. Hence, for $V^{(0)}{{\nu_{\bot}}}<0$ the equator rotates slower than the polar regions and just this condition is the result of the given numerical simulations. One could also demonstrate that the equator is accelerated for $V^{(0)}{{\nu_{\bot}}}>0$ (not shown).

If the associated meridional flow is also allowed to transport angular momentum, as done in the models in the second row of Fig. 7, the ${\it\Omega}$ isolines become cylindrical in accordance with the Taylor-Proudman theorem and the antisolar rotation law is only modified but not destroyed. We note how in the middle panels (only) the meridional circulation changes the type of the rotation law from uniform on spherical shells to cylindrical with respect to the rotation axis. The circulation cells always flow counterclockwise, that is, towards the poles at the surface.

From the comparison of Figs. 1 and 7 one finds that antisolar rotation profiles result both for positive ${\nu_{\bot}}$ in common with subrotation and/or for negative $H^{(1)}$. The latter can be excluded with numerical experiments where $U_{y}$ and $U_{z}$ are artificially suppressed mimicking the existence of strict solid-body rotation. This has been confirmed by Pencil Code444http://github.com/pencil-code simulations which use the same setup as in Käpylä (2019b); see Fig. 8.

Here we present results from a slowly rotating run with and without mean flows. The used Coriolis number ${\rm Co}={\it\Omega}d/\pi u_{\rm rms}\approx 0.25$ corresponds to ${\it\Omega}\simeq 4$ in Nirvana code units. We find that the signs of $Q_{xy}$ and $Q_{yz}$ change when the mean flows are suppressed. This is consistent with a dominating contribution from turbulent viscosity in $Q_{xy}$ and the rotation-induced off-diagonal viscosity term in $Q_{yz}$ in accordance with theory. Finally, we note that the magnitude of $Q_{xz}$ also increases clearly in the case where the shear is suppressed, indicating a strong contribution from the term involving $\nu_{\|}$.

5 Discussion

The cross-correlations $Q_{r\theta}$, $Q_{r\phi}$, and $Q_{\theta\phi}$ of the fluctuating velocities in a rotating turbulence are playing basic roles in understanding the rotation laws of stars with outer convection zones. Briefly, the tensor component $Q_{r\theta}$ transports thermal energy in the latitudinal direction (”warm poles”) while $Q_{r\phi}$ and $Q_{\theta\phi}$ transport angular momentum in the radial and meridional directions. For driven turbulence of a uniformly rotating density-stratified medium (stratified in the radial direction) the correlations fulfill simple sign rules independent of the rotation rate: it is $Q_{r\theta}<0$, $Q_{r\phi}<0$ and $Q_{\theta\phi}>0$, taken always in the northern hemisphere. In physical quantities this means that thermal energy is transported to the poles while the angular momentum is transported inward and towards the equator. The resulting warm poles drive a clockwise meridional circulation (northern hemisphere) which together with the positive $Q_{\theta\phi}$ transports angular momentum towards the equator. Hence, if the convection zone can be modelled by driven turbulence under fast global rotation then the resulting surface rotation law will always be of the solar-type. The simultaneous solution of the Reynolds equation and the corresponding energy equation provides rotation profiles, meridional circulation patterns, and pole–equator temperature differences that are relatively close to observations (Kitchatinov & Rüdiger 1999; Küker et al. 2011). The results are similar for all rotation rates and the number of tuning parameters is relatively low.

There are, however, increasing observational indications for the existence of antisolar rotation laws where the equator rotates slower than the polar regions. The results of observations (Metcalfe et al. 2016) as well as of numerical simulations (Gastine et al. 2014; Viviani et al. 2018) suggest that stars with slow rotation possess antisolar rotation laws. Based on photometric Kepler/K2 data from the open cluster M67, Brandenburg & Giampapa (2018) also argue in favor of the appearance of antisolar rotation laws for slow rotators with large Rossby numbers.

We assume that because of the often-stated positivity of the function $H$ the horizontal $\Lambda$ effect always transports angular momentum towards the equator in favor of an accelerated equator. If, however, a rotation law with a negative radial gradient exists then the rotation-induced off-diagonal components of the eddy viscosity tensor such as ${\nu_{\bot}}$ transport angular momentum towards the poles in favor of a polar vortex (“antisolar”). This also implies that a strictly radius-dependent rotation law ${\it\Omega}={\it\Omega}(r)$ can never exist in (slowly) rotating convection zones. After the Taylor-Proudman theorem, the rotation will tend to produce $z$-independent rotation laws, and therefore a negative radial shear of the angular velocity is always accompanied by slightly accelerated polar regions.

We compute the cross-correlations $Q_{r\theta}$, $Q_{r\phi}$, and $Q_{\theta\phi}$ from numerical simulations of convection in a rotating box. The averaging process concerns the horizontal planes, hence only radial shear can influence the cross-correlations (see Eqs. (4), (5), and (24)). For fast rotation, the well-known findings of positive $Q_{\theta\phi}$ and negative $Q_{r\theta}$, $Q_{r\phi}$ are reproduced (northern hemisphere). For slow rotation however, the signs of both $Q_{r\theta}$ and $Q_{\theta\phi}$ change almost simultaneously, meaning that the angular momentum is now transported to the poles and the heat is transported to the equator. The resulting warm equator leads to a counterclockwise meridional circulation which transports the angular momentum to the poles. The new signs of the quantities $Q_{r\theta}$ and $Q_{\theta\phi}$ may thus lead to antisolar differential rotation.

This behavior, however, is due to the appearance of large-scale flows in zonal and in meridional directions. The zonal flow $U_{z}$ mimics differential rotation with negative radial gradient. Via the off-diagonal viscosity in Eq. (5), for positive ${\nu_{\bot}}$, a negative contribution to $Q_{\theta\phi}$ leads to overcompensation of the positive but small values of $H$. A similar effect happens for $Q_{r\theta}$ where the negative radial gradient of $U_{y}$ combined with the positive eddy viscosity $\nu_{\|}$ provides positive contributions to the negative cross-correlation (21). From the associated numerical values, with Eq. (25) $\nu_{\|}$ can also be calculated leading to $\nu_{\|}/\nu\simeq 240$, or to the (reasonable) dimensionless quantity \alpha_{\rm vis}\lower 1.72218pt\hbox{;\buildrel<\over{\scriptstyle\sim};}0.5. The simulations lead to the relation ${\nu_{\bot}}/\nu_{\|}\simeq 3\tau^{2}_{\rm corr}$ for the ratio of the rotation-induced off-diagonal viscosity term and the standard (diagonal) eddy viscosity.

We also studied the influence of the radial shears of the large-scale flows $U_{y}$ and $U_{z}$ on the cross-correlations $Q_{xy}$ and $Q_{yz}$ with numerical experiments where the $U_{y}$ and $U_{z}$ can artificially be suppressed. In these cases, for all rotation rates, our calculations led to negative $Q_{xy}$ and positive $Q_{yz}$. Without large-scale flows the analytical results are confirmed, namely that $Q_{r\theta}<0$ and $Q_{\theta\phi}>0$ for the northern hemisphere formulated in spherical coordinates. For slow rotation, both signs are changed if the large-scale shear flows are allowed to back-react.

In summary, we show that a so-far neglected rotation-induced off-diagonal eddy viscosity term combined with rotation laws with a negative radial gradient (subrotation, as existing in slowly rotating stars) is able to produce differential rotation of the antisolar type. This result complies with the (numerical) findings of Viviani et al. (2018) that negative radial ${\it\Omega}$-gradients and antisolar differential rotation are closely related. Therefore, if new observations confirm the existence of decelerated equators at the surface of slowly rotating stars then we shall better understand the nature of the active eddy viscosity tensor and the underlying turbulence model.

Acknowledgements.

PJK acknowledges the computing resources provided by CSC – IT Center for Science, who are administered by the Finnish Ministry of Education; of Espoo, Finland, and the Gauss Center for Supercomputing for the Large-Scale computing project “Cracking the Convective Conundrum” in the Leibniz Supercomputing Centre’s SuperMUC supercomputer in Garching, Germany. This work was supported in part by the Deutsche Forschungsgemeinschaft Heisenberg programme (grant No. KA 4825/1-1; PJK) and the Academy of Finland ReSoLVE Centre of Excellence (grant No. 307411; PJK).