Astron. Nachr. / AN 999, No. 88, 789 – 795 (2011) / DOI please set DOI!

Transport of angular momentum and chemical species by anisotropic mixing in stellar radiative interiors

L.L. Kitchatinov^{1,2,3}^{,}^{⋆} and A. Brandenburg^{3,4}

1 Institute for Solar-Terrestrial Physics, P.O. Box 291, Irkutsk 664033, Russia

2 Pulkovo Astronomical Observatory, St. Petersburg 196140, Russia

3 Nordita, AlbaNova University Center, Roslagstullsbacken 23, SE-10691 Stockholm, Sweden

4 Department of Astronomy, Stockholm University, SE-10691 Stockholm, Sweden

The dates of receipt and acceptance should be inserted later

Key words stars: interiors – stars: rotation – stars: abundances – hydrodynamics – turbulence

Small levels of turbulence can be present in stellar radiative interiors due to, e.g., instability of rotational shear. In this paper we estimate turbulent transport coefficients for stably stratified rotating stellar radiation zones. Stable stratification induces strong anisotropy with a very small ratio of radial-to-horizontal turbulence intensities. Angular momentum is transported mainly due to the correlation between azimuthal and radial turbulent motions induced by the Coriolis force. This non-diffusive transport known as the Λ-effect has outward direction in radius and is much more efficient compared to the effect of radial eddy viscosity. Chemical species are transported by small radial diffusion only. This result is confirmed using direct numerical simulations combined with the test-scalar method. As a consequence of the non-diffusive transport of angular momentum, the estimated characteristic time of rotational coupling ( < ∼ 100 Myr) between radiative core and convective envelope in young solar-type stars is much shorter compared to the time-scale of Lithium depletion (∼ 1 Gyr).

\begin{matrix}c \\ ◯ \end{matrix} 2011 WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim

1 Introduction

Mixing of material in stably stratified stellar radiative interiors is opposed by buoyancy forces. Nevertheless, a certain level of turbulence resulting, e.g., from instabilities of rotational shear (Goldreich & Schubert 1967) or forcing from the upper convection zone (Blöcker et al. 1998) is possible and even necessary to explain rotational coupling between the radiative core and the convective envelope, and thus the depletion of light elements in solar-type stars.

Turbulent transport in stellar radiation zones is supposed to be strongly anisotropic (Spiegel & Zahn 1992; Denissenkov 2010). Buoyancy forces suppress radial mixing so that much higher horizontal compared with radial velocities might be expected. The degree of the anisotropy is not certain, however. It remains a free parameter of transport models. The anisotropy in rotating fluids is, however, not free (Rüdiger & Pipin 2001). The influence of the Coriolis force on horizontal motions produces radial displacements. The ratio of radial to azimuthal mixing intensities is controlled by the balance between Coriolis and buoyancy forces.

The aim of this paper is to estimate the effects of rotation and stable stratification on turbulence. We do not specify here the origin of the turbulence but just prescribe the so-called ‘original turbulence’ that would take place in nonrotating and neutrally stratified fluid. The combined effect of rotation and stable stratification on this original turbulence is then estimated.

We shall see that the resulting parameters of turbulent mixing depend only slightly on the properties of the prescribed original turbulence. In particular, the anisotropy,

\frac{⟨u_{ϕ}^{2}⟩}{⟨u_{r}^{2}⟩} ≈ \frac{τ^{2}N^{4}}{Ω^{2}} , (1)

varies only slightly between very different prescriptions (in this equation, u is the fluctuating velocity, Ω is the angular velocity, N is the buoyancy frequency, and τ is the eddy turnover time). Typically, Ω ≪ N in radiation zones and probably we also have τN ≫ 1. Therefore, the anisotropy of Eq. (1) is high. Our main finding is that the strongly anisotropic turbulence is much more efficient in transporting angular momentum than chemical species. The reason is that the Coriolis force makes the azimuthal and radial motions correlated. Fluid particles with positive (negative) azimuthal velocity are pushed by the Coriolis force in the positive (negative) radial direction. As a result, ⟨u_{r}u_{ϕ}⟩ > 0, and the turbulence transports angular momentum in the positive radial direction. This angular momentum flux cannot be interpreted as an effect of eddy viscosity, i.e., it is not diffusive by nature. Transport of chemical species, on the contrary, is only due to radial turbulent mixing.

The non-diffusive transport of angular momentum that we find is essentially the well-known Λ-effect of differential rotation theory (Lebedinskii 1941; Rüdiger 1989). The effect has been well studied for nearly adiabatically stratified stellar convection zones. Our paper suggests that it can be important for the rotational coupling between radiative interiors and convective envelopes as well.

\begin{matrix}c \\ ◯ \end{matrix} 2011 WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim

⋆ Corresponding author: e-mail: kit@iszf.irk.ru

2 Mathematical formulation

2.1 Relation between the original and background turbulence

To specify the effects of rotation and stable stratification on turbulence, we follow the standard approach of quasilinear theory by formulating the linear relation, u_{i} = D_{ij}u_{j}^{(0)} between the velocity fields of original turbulence u^{(0)}, which would take place in nonrotating neutrally stratified fluid, and actual or background turbulence u. The relation tensor D_{ij} includes the effects of rotation and stratification.

We apply a simple version of the τ-approximation to write the equations for fluctuating velocity and entropy as follows,

\frac{u}{τ} + 2Ω × u + \frac{1}{ρ}∇P − g = f,
\frac{s′}{τ} = −(u · ∇)⟨s⟩,  (2)

where f is a random force driving the turbulence, τ is the eddy turnover time, g is gravity and the τ-relaxation terms replace the nonlinear terms together with time-derivatives,

\frac{∂u}{∂t} + (u · ∇)u − ⟨(u · ∇)u⟩ → \frac{u}{τ},

\frac{∂s′}{∂t} + u · ∇s − ⟨u · ∇s⟩ → \frac{s′}{τ} .  (3)

Microscopic diffusion is neglected. The τ-approximation

(3) assumes that the nonlinear interaction of turbulent eddies leads to their effective dissipation (a turbulent fragmentation of scales) in a characteristic time τ ≈ ℓ/u, where ℓ is the characteristic spatial scale of the eddies. We assume that the motion is incompressible, divu = 0. This means that vertical displacements are small compared to the density scale height and horizontal velocities are small compared to the speed of sound.

Fourier transformation is applied,

u(r) = ∫ e^{ik·r}ũ(k) dk,
s′(r) = ∫ e^{ik·r}s̃′(k) dk,  (4)

to convert the partial differential equations (2) into algebraic equations. We neglect variations of mean fields on the spatial scale of the turbulence to obtain a closed equation for the fluctuating velocity:

ũ_{i} + N^{∗}^{2}(r̂_{i} − μk̂_{i})(r̂ · ũ) + σΩ^{∗}ε_{ipm}k̂_{p}ũ_{m} = ũ_{i}^{(0)} (5)

(cf. Kitchatinov & Rüdiger (2008) on how the pressure term is treated). In Eq. (5), r̂ is the radial unit vector, k̂ = k/k is the unit vector along the wave vector, μ = (k̂ · r̂) is the cosine of the angle between wave vector and radius, and σ = (k̂ · Ω)/Ω is cosine of the angle between the wave vector and rotation axis. The Coriolis number

Ω^{∗} = 2τΩ (6) and the normalized buoyancy frequency

N^{∗} = τN, N^{2} = \frac{g}{c_{p}} \frac{∂⟨s⟩}{∂r}, (7)

parameterize the effects of rotation and stratification. Solving Eq. (5) for ũ gives the relation tensor,

D_{ij} = \frac{1}{1 + N^{∗}^{2}(1 − μ^{2}) + σ^{2}Ω^{∗}^{2}} ×
( (1 + N^{∗}^{2}(1 − μ^{2})) (δ_{ij} + σΩ^{∗}ε_{ijp}k̂_{p}) −
− N^{∗}^{2}(r̂_{i} − μk̂_{i})r̂_{j} +

+ σΩ^{∗}N^{∗}^{2} ((r̂_{i} − μk̂_{i})ε_{jmp} − r̂_{j}ε_{imp}) r̂_{m}k̂_{p}),

ũ_{i} = D_{ij}ũ_{j}^{(0)}.  (8)

This equation describes the joint influence of rotation and stable stratification on the turbulence.

In the limit of neutral stratification, N^{∗} → 0, Eq. (8) reduces to the familiar expression (Kitchatinov 1986)

D_{ij} = \frac{δ_{ij} + σΩ^{∗}ε_{ijp}k̂_{p}}{1 + σ^{2}Ω^{∗}^{2}} , (9)

describing the influence of rotation in the τ-approximation. In the other limit of slow rotation, Ω^{∗} → 0, Eq. (8) reduces to

D_{ij} = δ_{ij} − \frac{N^{∗}^{2}}{1 + N^{∗}^{2}(1 − μ^{2})}(r̂_{i} − μk̂_{i})r̂_{j}. (10)

This equation accounts for the effect of stable stratification alone.

2.2 Original turbulence models

The original turbulence properties can be prescribed by specifying the spectral tensor Q̃_{ij}^{(0)}:

⟨ũ_{i}^{(0)}(k)ũ_{j}^{(0)}(k′)⟩ = Q̃_{ij}^{(0)}(k) δ(k + k′) (11)

(Rüdiger 1989). The relation tensor (8) can then be used to account for the effects of rotation and stratification,

Q̃_{ij} = D_{im}D_{jn}Q̃_{mn}^{(0)},
Q_{ij} = ⟨u_{i}u_{j}⟩ = ∫ Q̃_{ij}(k) dk.  (12)

For the simplest case of isotropic nonhelical turbulence the spectral tensor reads

Q̃_{ij}^{(0)} = \frac{E(k)}{8πk^{2}} (δ_{ij} − k̂_{i}k̂_{j}) . (13)

In this equation, E is the spectrum function

⟨u^{2}⟩^{(0)} = ∫_{0}^{∞}E(k)dk. (14)

A more realistic model is anisotropic turbulence with a preferred direction being the radial one. Horizontal velocities are expected to be much larger than the vertical ones.

For the extreme case of strictly horizontal random motions, we have

Q̃_{ij}^{(0)} = q(k, μ)×

((1 − μ^{2})(δ_{ij} − k̂_{i}k̂_{j}) − (r̂_{i} − μk̂_{i})(r̂_{j} − μk̂_{j})). (15)

The correlation lengths for vertical and horizontal directions may differ. If, however, the lengths are equal, the q-function in Eq. (15) does not depend on μ:

q = \frac{3E(k)}{8πk^{2}}. (16)

The case where the correlation length in radius is small compared to the horizontal correlation length can be modeled by the equation

q = \frac{3E(k_{⊥}) k_{⊥}}{8πk^{3}} , k_{⊥} = k\sqrt{1 − μ^{2}}. (17)

The spectrum functions E of Eqs. (16) and (17) are related to the turbulence intensity by the same equation (14) as in the case of the isotropic turbulence of Eq. (13).

Anisotropic turbulence that has finite radial velocities (but different from the horizontal velocities) can be modeled by a linear superpositions of the spectral tensors (13) and (15). We shall see that, if the buoyancy frequency is large compared to the rotation frequency, N ≫ Ω, and the normalized buoyancy frequency is large, N^{∗} ≫ 1, which are conditions typical of stellar radiation zones, then the turbulent transport parameters differ little between the representations (13), (16) and (17) for the original turbulence, i.e., the transport characteristics are not sensitive to a particular choice of the original turbulence model.

2.3 Direct numerical simulations

An independent verification of the effects of strong stratification on turbulent diffusion is provided by means of direct numerical simulations. In that case, we solve the full set of compressible hydrodynamic equations for ρ, u, and s:

\frac{Dρ}{Dt} = −ρ∇ · u,  (18)
ρ\frac{Du}{Dt} = −∇P + ∇ · (2ρνS) + ρ(f + g),  (19)
ρT \frac{Ds}{Dt} = ∇ · K∇T + 2ρνS^{2} − \frac{1}{τ_{s}} ρ(c_{s}^{2} − c_{s0}^{2}),  (20)

where D/Dt = ∂/∂t + u · ∇ is the advective derivative, S_{ij} = \frac{1}{2}(u_{i,j} + u_{j,i}) − \frac{1}{3}δ_{ij}∇ · u is the traceless rate of strain tensor, commas indicate partial differentiation, ν is the kinematic viscosity, the specific entropy is given by s = c_{v} ln P − c_{v} ln ρ, where c_{p} and c_{v} are the specific heats at constant pressure and constant volume, respectively, the temperature is related to P and ρ via (c_{p} − c_{v})T = P/ρ, which, in turn, is related to the sound speed c_{s} via c_{s}^{2} = γP/ρ, where γ = c_{p}/c_{v} is the ratio of specific heats, f is the external forcing function, and K is the thermal conductivity. The last term in the entropy equation (20) represents a cooling term that keeps the temperature (or sound speed) approximately constant with a given relaxation or cooling time τ_{s}. The flow is driven by a random forcing function consisting of nonhelical waves with wavenumbers whose modulus lie in a narrow band around an average wavenumber k_{f} (Haugen et al. 2004). We arrange the amplitude of the forcing function such that the Mach number based on the rms velocity remains below 0.1, so the effects of compressibility are negligible.

We consider a cubic domain of size L^{3} with periodic boundary conditions in the x and y directions and insulating impenetrable stress-free boundary conditions on z = ±L/2. The mass in the volume is therefore conserved and given by ∫ ρ d^{3}x = ρ_{0}L^{3}. The lowest wavenumber that fits into the domain is k_{1} = 2π/L. Our average forcing k_{f} wavenumber is chosen such that k_{f}/k_{1} = 5. We adopt isothermal initial conditions with c_{s} = c_{s0} and ρ = ρ_{0} exp(−z/H), where H = c_{s0}^{2}/γg is the scale height which is chosen such that Hk_{1} = 1. The turnover time based on the wavenumber k_{f} is τ = (u_{rms}k_{f})^{−1}, where u_{rms} is the rms velocity based on all three velocity components. We vary the normalized buoyancy frequency by varying the forcing amplitude. For each run the viscosity is adjusted such that the Reynolds number, Re = u_{rms}/νk_{f}, is around 40. The Prandtl number, Pr = νc_{p}ρ/K, is chosen to be equal to unity. The cooling time is chosen such that τ_{s}u_{rms}k_{f} is also around unity.

To quantify the suppression of vertical mixing in a numerical simulation of stratified turbulence we use the testscalar method by solving equations for the fluctuating passive scalar concentration in the presence of a prescribed mean passive scalar concentration,

\frac{∂c^{pq}}{∂t} = −∇ · (uC^{pq} + uc^{pq} − ⟨uc^{pq}⟩) + κ∇^{2}c^{pq}, (21)

where u is the velocity fluctuation obtained from Eq. (19), and the mean flow is zero. The superscript p (=1 or 3) stands for test scalars varying in the x or z directions while q stands for c and s that represent mean fields that are proportional to cos kx_{p} or sin kx_{p}, respectively,

C^{pc} = C_{0} cos kx_{p}, C^{ps} = C_{0} sin kx_{p}, (22)

where (x_{1}, x_{2}, x_{3}) = (x, y, z) are Cartesian coordinates and C_{0} is a normalization factor. Angle brackets denote planar averaging over the yz plane (if p = 1) or the xy plane (if p = 3). The passive scalar flux is then related to the gradient of C^{pc} via ⟨u_{i}c^{pq}⟩ = −κ_{ij}∇_{j}C^{pc}, where κ_{ij} are the components of the turbulent diffusivity tensor, obtained as

κ_{xx} = (sin kx⟨u_{x}c^{1c}⟩ − cos kx⟨u_{x}c^{1s}⟩)/kC_{0},  (23)
κ_{zz} = (sin kz⟨u_{z}c^{1c}⟩ − cos kz⟨u_{z}c^{1s}⟩)/kC_{0},  (24)

by using the solutions for four different test scalars. These coefficients depend on time and one spatial coordinate, but because of stationarity and approximate homogeneity of the turbulence intensity, we present in the following temporal and spatial averages of these coefficients.

The test scalar method has been used previously to quantify mixing in turbulence in the presence of rotation and magnetic fields (Brandenburg et al. 2009), shear (Madarassy & Brandenburg 2010), as well as isothermal density stratification (Brandenburg et al. 2012). However, unlike those earlier works, the entropy equation is here included and a non-isothermal equation of state for a perfect monatomic gas is used with γ = 5/3. Equations (18)–(21) are solved using the Pencil Code^{1} in three dimensions. The numerical resolution used for the simulations presented here is 128^{3} meshpoints. Like the Prandtl, the Schmidt number, Sc = ν/κ, is also chosen to be equal to unity.

3 The effect of stable stratification

We consider first the effect of stable stratification for a nonrotating fluid, Ω = 0. The only preferred direction in this case is the radial one. For any model of the original turbulence, we therefore have

⟨u_{i}u_{j}⟩ = \frac{⟨u_{h}^{2}⟩}{2} (δ_{ij} − r̂_{i}r̂_{j}) + ⟨u_{r}^{2}⟩r̂_{i}r̂_{j}, (25)

where ⟨u_{h}^{2}⟩ = ⟨u_{ϕ}^{2}⟩ + ⟨u_{θ}^{2}⟩ is the horizontal turbulence intensity. Standard spherical coordinates (r, θ, ϕ) are used.

On using the relation tensor (10) and Eq. (12) for the isotropic original turbulence of Eq. (13), we find

⟨u_{h}^{2}⟩ = ⟨u^{2}⟩^{(0)}(\frac{1}{2} + \frac{1}{4N^{∗}^{2}} −

−\frac{1}{8N^{∗}^{3}\sqrt{N^{∗}^{2} + 1}} ln (\frac{\sqrt{N^{∗}^{2} + 1} + N^{∗}}{\sqrt{N^{∗}^{2} + 1} − N^{∗}} )) , (26)

⟨u_{r}^{2}⟩ = \frac{⟨u^{2}⟩^{(0)}}{4N^{∗}^{2}(N^{∗}^{2} + 1)} ×

(\frac{2N^{∗}^{2} + 1}{2N^{∗}\sqrt{N^{∗}^{2} + 1}} ln (\frac{\sqrt{N^{∗}^{2} + 1} + N^{∗}}{\sqrt{N^{∗}^{2} + 1} − N^{∗}} ) − 1) . (27)

Figure 1 shows the ratio of vertical to horizontal turbulence intensities as a function of the normalized buoyancy frequency N^{∗}. Squares on the plot show the results of direct numerical simulations described in Sect. 2.3.

Horizontal turbulence intensity is influenced by the stable stratification only slightly. In the limit of very large N^{∗}, it is ⟨u_{h}^{2}⟩ = ⟨u^{2}⟩^{(0)}/2, i.e., ⟨u_{h}^{2}⟩ is reduced compared to the case of a non-stratified fluid by a factor of 3/4 only. The vertical motions, by contrast, are strongly suppressed:

⟨u_{r}^{2}⟩ ≃ \frac{⟨u^{2}⟩^{(0)}}{4N^{∗}^{4}} (ln(4N^{∗}^{2}) − 1) for N^{∗} ≫ 1 (28)

in the strong stratification limit. The originally isotropic turbulence is changed towards horizontal turbulence as N^{∗} increases.

As might be expected, the horizontal turbulence of Eq.

(15) is not influenced by the stratification at all, i.e., Q_{ij} = Q_{ij}^{(0)} in this case. Therefore, whatever model for the original turbulence in Section 2.2 is used, the resulting turbulence in the limit of large N^{∗} is almost the same. The only difference between isotropic and horizontal original turbulence models

Fig. 1 Ratio of the vertical to horizontal turbulence intensities in dependence on the normalized buoyancy frequency (7) for the case of a non-rotating fluid and isotropic original turbulence. Squares show the results of 3D direct numerical simulations for isotropically forced turbulence.

Fig. 2 Dependence of κ_{xx}/κ_{t0} (open symbols) and κ_{zz}/κ_{t0} (filled symbols) on the normalized buoyancy frequency. The dashed line shows that κ_{xx}/κ_{t0} ≈ 2.5 while the solid line gives κ_{zz}/κ_{t0} ≈ 0.09 (τN)^{−3/2}.

is that for the first case the radial velocities of Eq. (28) are finite, although very small.

Similar results are obtained numerically for the turbulent diffusivity using the test-scalar method; see Fig. 2, where we plot κ_{xx} and κ_{zz}, both normalized by the reference value κ_{t0} = u_{rms}/3k_{f}. Note that κ_{xx}/κ_{t0} ≈ 2.5, independent of τN, while κ_{zz}/κ_{t0} ≈ 0.09 (τN)^{−3/2}. Visualizations of the scalar quantity τ∇ · u on the periphery of the computational domain (Fig. 3) confirm our expectation that turbulent structures that become flatter as τN is increased. These structures become smoother in the xy plane, while their length scale in the z direction decreases.

1 http://www.pencil-code.googlecode.com \begin{matrix}c \\ ◯ \end{matrix} 2011 WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim www.an-journal.org

Fig. 3 Visualizations of τ∇ · u on the periphery of the computational domain for τN = 0.17, 0.9, and 3.7.

4 Mixing in rotating and stratified fluids

Similarity of the results for different original turbulence models is even more pronounced when rotation is included. Coriolis force produces radial motion even for original horizontal turbulence. The relative intensity of radial mixing is controlled by the balance between buoyancy and Coriolis forces.

Producing analytical results for arbitrary values of N^{∗} and Ω^{∗} is problematic. Derivations were performed for the practically interesting case of strong stratification, N^{∗} ≫ 1, and not too fast rotation, Ω^{∗} ≪ N^{∗}. The period of gravity waves for radiation zones of solar-type stars is of the order of one hour. The characteristic time of turbulence is probably much longer (it is about one month for convection and probably longer for radiation zone motions), so that the condition N^{∗} ≫ 1 is well satisfied. The condition Ω^{∗} ≪ N^{∗} should be satisfied as well in radiation zones of pressure supported stars (otherwise we are dealing with centrifugally supported disks).

We now use Eq. (8) for the relation tensor that includes the combined effects of rotation and stratification. For the case of the isotropic original turbulence of Eq. (13), the terms of lowest order in N^{∗}^{−1} in the most significant velocity correlations read

⟨u_{h}^{2}⟩ = \frac{1}{2}⟨u^{2}⟩^{(0)},
⟨u_{r}^{2}⟩ = \frac{⟨u^{2}⟩^{(0)}}{4N^{∗}^{4}} [ln (\frac{4N^{∗}^{2}}{1 + Ω^{∗}^{2} cos^{2} θ} ) − 1 + Ω^{∗}^{2}
+ Ω^{∗}^{2} cos^{2} θ (ln (\frac{4N^{∗}^{2}}{1 + Ω^{∗}^{2} cos^{2} θ} ) − 4)] ,
⟨u_{r}u_{θ}⟩ = ⟨u^{2}⟩^{(0)} \frac{Ω^{∗}^{2}}{4N^{∗}^{4}} (ln (\frac{4N^{∗}^{2}}{1 + Ω^{∗}^{2} cos^{2} θ} ) − 3)

× sin θ cos θ,

⟨u_{r}u_{ϕ}⟩ = ⟨u^{2}⟩^{(0)} \frac{Ω^{∗}}{4N^{∗}^{2}} sin θ.  (29)

The turbulence intensities ⟨u_{h}^{2}⟩ and ⟨u_{r}^{2}⟩ control the eddy diffusion in latitude and radius, respectively. The crosscorrelation ⟨u_{r}u_{ϕ}⟩ is important for transport of angular momentum. The correlation ⟨u_{r}u_{θ}⟩ may cause the temperature variation with latitude; positive ⟨u_{r}u_{θ}⟩ cos θ implies poleward eddy heat flux. It may be noted that rotation does not produce anisotropy in the horizontal plane in the strong stratification limit, ⟨u_{θ}^{2}⟩ = ⟨u_{ϕ}^{2}⟩ = ⟨u_{h}^{2}⟩/2 (more precisely, the anisotropy is small: (⟨u_{θ}^{2}⟩ − ⟨u_{ϕ}^{2}⟩)/⟨u_{h}^{2}⟩ ∼ Ω^{∗}^{2}/N^{∗}^{4}).

The results for horizontal turbulence of Eq. (15) with an isotropic correlation length of Eq. (16) or a short vertical correlation length of Eq. (17) are the same and read

⟨u_{h}^{2}⟩ = ⟨u^{2}⟩^{(0)},
⟨u_{r}^{2}⟩ = ⟨u^{2}⟩^{(0)} \frac{Ω^{∗}^{2}}{2N^{∗}^{4}} ,
⟨u_{r}u_{θ}⟩ = ⟨u^{2}⟩^{(0)} \frac{Ω^{∗}^{2}}{2N^{∗}^{4}} sin θ cos θ,
⟨u_{r}u_{ϕ}⟩ = ⟨u^{2}⟩^{(0)} \frac{Ω^{∗}}{2N^{∗}^{2}} sin θ.  (30)

If the Coriolis number (6) is not small (Ω^{∗} > 1), equations (29) and (30), though different in details, are of the same order of magnitude. As the τ-approximation, which was used to derive these equation, probably has the same accuracy, the results are practically the same. We focus in the following discussion on the results of Eq. (30) because the case of horizontal original turbulence is probably more adequate for stably stratified fluids.

5 Discussion

A remarkable feature of Eq. (30) is the finite crosscorrelation ⟨u_{r}u_{ϕ}⟩. This means that the turbulence transports angular momentum in radius. The cross-correlation is positive so that the angular momentum is transported outward. For the case of differential rotation caused by stellar spin-down, this angular momentum flux acts in same direction as turbulent viscosity, but it is not viscous by nature.

What we find here is the known Λ-effect of differential rotation theory (Lebedinskii 1941; Rüdiger 1989). Figure 4 illustrates its origin. The positive correlation ⟨u_{r}u_{ϕ}⟩ results from the Coriolis force action on the original horizontal motion. The fluid particles, which move in the direction of global rotation, are deflected outwards in radius by the Coriolis force. Fluid particles moving in the retrograde direction are deflected inwards. The product u_{r}u_{ϕ} is positive in both cases. Note that the influence of the Coriolis force

Fig. 4 Illustration of the origin of non-viscous angular momentum transport by anisotropic turbulence (see text).

on the original radial motion produces a negative correlation ⟨u_{r}u_{ϕ}⟩ < 0 (Fig. 4). The sense of radial angular momentum transport is controlled by the turbulence anisotropy. Angular momentum is transported inward if radial mixing prevails and outward for predominantly horizontal mixing (Rüdiger 1989). In the present case, the anisotropy of horizontal type is produced by the stable stratification. It may be noticed that we do not find a meridional Λ-effect. The meridional flux of angular momentum is of higher order in N^{∗}^{−2} compared to the radial flux. This is a consequence of small anisotropy in the horizontal plane. Therefore, mixing in stellar radiation zones does not produce latitudinal differential rotation.

In a star during its spinning-down, turbulent viscosity and Λ-effect both transport angular momentum outward. If the turbulent viscosity in the radial direction is estimated as ν_{T} ≈ τ⟨u_{r}^{2}⟩, the ratio of angular momentum fluxes produced by the viscosity and the Λ-effect can be estimated as

\frac{Q_{rϕ}^{ν}}{Q_{rϕ}^{Λ}} = \frac{−ν_{T}r \frac{∂Ω}{∂r} sin θ}{⟨u_{r}u_{ϕ}⟩} = q (\frac{Ω}{N} )^{2} ≪ 1,
q = −\frac{r}{Ω} \frac{dΩ}{dr},  (31)

where Eq. (30) was used. For the upper radiation zone of the Sun, Ω^{2}/N^{2} ∼ 10^{−5} (see Fig. 1 in Kitchatinov & Rüdiger 2008). For faster rotating young stars the ratio is larger, but still well below unity. Therefore, the angular momentum transport by the Λ-effect is much more efficient compared to the effect of the eddy viscosity.

Radial transport of chemical species is slow because it is only due to weak radial mixing and it is further reduced by relatively intensive horizontal mixing (Vincent, Michaud & Meneguzzi 1996; Michaud & Zahn 1998). This may explain why the characteristic times of rotational coupling between the core and the convective envelope ( < ∼ 100 Myr; Hartmann & Noyes 1987; Denissenkov et al. 2010) in solar-type stars is shorter than the characteristic time (∼ 1 Gyr; Skumanich 1972; Meléndez et al., 2010; Baumann et al. 2010) of Lithium depletion.

Papaloizou & Pringle (1978) found that differential rotation in radiative cores is unstable to r-modes or vortices that are global in latitude but of small scale in radius. For the estimates below we assume that the turbulent motions are global in the horizontal directions and that radial eddy diffusion of chemicals is estimated as

χ_{T} ≈ τ⟨u_{r}^{2}⟩ ≈ τ⟨u^{2}⟩^{(0)} \frac{Ω^{∗}^{2}}{2N^{∗}^{4}} ≈ \frac{2R^{2}Ω^{2}}{τ^{3}N^{4}} , (32)

where again Eq. (30) was used. Lithium in young solar-type stars is depleted about ten times in the first billion years of their main-sequence life (Meléndez et al. 2010). Lithium has to be transported over a short distance (∼ 0.1R) below the convection zone in order to be destroyed (see Fig. 1 of Rüdiger & Pipin 2001). This leads to an estimated fractional Lithium abundance (relative to the primordial abundance) at the age of 1 Gyr

exp (−\frac{10^{9} yr × χ_{T}}{10^{−2}R^{2}} ) ≈ 0.1. (33)

Taking the value of R ≈ 5 × 10^{10} cm for the radius of the base of (solar) convection zone, we get from Eq. (33)

χ_{T} ≈ 2 × 10^{3} cm^{2} s^{−1}. (34)

Equations (32) and (34) lead to an estimate of the characteristic eddy turnover time of

τ ≈ 10^{7} s (35)

for a young star rotating ten times faster than the Sun. With this value of τ, an estimate for the characteristic time of core-envelope rotational coupling due to the Λ-effect can be obtained,

T_{Ω} ≈ \frac{R^{2}Ω}{Q_{rϕ}^{Λ}} ≈ τ^{3}N^{2} ≈ 60 Myr, (36)

where again Eq. (30) for the non-diffusive flux of angular momentum Q_{rϕ}^{Λ} = ⟨u_{r}u_{ϕ}⟩ was used together with the estimate N ≈ 10^{−3} s^{−1} for radiation zones.

The estimates suggest the following scenario for angular momentum and chemical species transport in stellar radiation zones. Angular momentum loss by a stellar wind induces differential rotation with an angular velocity decreasing outwards. The differential rotation is hydrodynamically unstable and this instability produces turbulent mixing. The stable stratification makes the turbulence highly anisotropic with very small radial velocities. The eddy diffusion in radius lead to a slow decrease of the Lithium abundance with time. Angular momentum is transported not by the eddy viscosity but by the Λ-effect due to the correlation of azimuthal and radial turbulent motions induced by the Coriolis force. The non-viscous outward transport provides relatively fast rotational coupling between the radiative core and the convective envelope.

Other possibilities for rotational coupling discussed in the literature include the angular momentum transport by gravity waves (Charbonnel & Talon 2005) and by magnetic stress (Charbonneau & MacGregor 1993; Rüdiger & Kitchatinov 1996; Denissenkov 2010). Further studies are necessary to decide which of the possibilities are most relevant.

Acknowledgements. LLK is thankful to NORDITA for hospitality and support. This work was supported by the Russian Foundation for Basic Research (projects 10-02-00148, 10-02-00391) and the European Research Council under the AstroDyn Research Project 227952. We acknowledge the allocation of computing resources provided by the Swedish National Allocations Committee at the Center for Parallel Computers at the Royal Institute of Technology in Stockholm and the National Supercomputer Centers in Linköping.

References

Baumann, P., Ramirez, I., Meléndez, J., Asplund, M., Lind, K.:

2010, A&A 519, A87 Blöcker, T., Holweger, H., Freitag, B., Herwig, F., Ludwig, H.-G.,

Steffen, M.: 1998, SSRev 85, 105 Brandenburg, A., Rädler, K.-H., Kemel, K.: 2012, A&A (in press,

arXiv:1108.2264) Brandenburg, A., Svedin, A., Vasil, G.M.: 2009, MNRAS 395,

1599 Charbonneau, P., MacGregor, K.B.: 1993, ApJ 417, 762 Charbonnel, C., Talon, S.: 2005, Science 309, 2189 Denissenkov, P.A.: 2010, ApJ 719, 28 Denissenkov, P.A., Pinsonneault, M., Terndrup, D.M., New-

sham, G.: 2010, ApJ 716, 1269 Goldreich, P., Schubert, G.: 1967, Science 156, 1101 Hartmann, L.W., Noyes, R.W.: 1987, Ann. Rev. Astron. Astro-

phys. 25, 271 Haugen, N.E.L., Brandenburg, A., Dobler, W.: 2004, Phys. Rev. E

70, 016308 Kitchatinov, L.L.: 1986, GApFD 35, 93 Kitchatinov, L.L., Rüdiger, G.: 2008, A&A 478, 1 Lebedinskii, A.I.: 1941, Astron. Zh., 18, 10 Madarassy, E.J.M. and Brandenburg, A.: 2010, Phys. Rev. E 82,

016304 Meléndez, J., Ramirez, I., Casagrande, L. et al.: 2010, Ap&SS 328,

193 Michaud, G., Zahn, J.-P.: 1998, Theoret. Comput. Fluid Dynamics

11, 183 Papaloizou, J., Pringle, J.E.: 1978, MNRAS 182, 423 Rüdiger, G.: 1989, Differential Rotation and Stellar Convection,

Gordon & Breach, New York Rüdiger, G., Kitchatinov, L.L.: 1996, ApJ 466, 1078 Rüdiger, G., Pipin, V.V.: 2001, A&A 375, 149 Skumanich, A.: 1972, ApJ 171, 565 Spiegel, E.A., Zahn, J.-P.: 1992, A&A 265, 106 Vincent, A., Michaud, G., Meneguzzi, M.: 1996, Phys. Fluids 8,

1312