Mean-field theory of differential rotation in density stratified turbulent convection
I. Rogachevskii N. Kleeorin
physics.flu-dyn astro-ph.SR physics.plasm-ph physics.space-ph


[
    Tulasi Ram Reddy
    December 30, 2023
=====================






A mean-field theory of differential rotation in a density stratified turbulent convection
has been developed.
This theory is based on a combined effect of the turbulent heat flux and
anisotropy of turbulent convection on the Reynolds stress.
A coupled system of dynamical budget equations consisting in
the equations for the Reynolds stress, the entropy fluctuations
and the turbulent heat flux has been solved.
To close the system of these equations, the spectral tau approach
which is valid for large Reynolds and Peclet numbers, has been applied.
The adopted model of the background turbulent convection takes into account an
increase of the turbulence anisotropy and a decrease of the turbulent
correlation time with the rotation rate.
This theory yields the radial profile of the differential
rotation which is in agreement
with that for the solar differential rotation.




§ INTRODUCTION


Origin of the solar and stellar magnetic fields is associated with
a mean-field dynamo (refereed as αΩ
or α^2 Ω dynamos) that is based on the combined effect of
helical turbulent motions and a differential rotation
<cit.>.
A non-zero mean kinetic helicity produced by a rotating density stratified turbulent
convection, causes the α effect in the solar convective zone.
An origin of the solar differential rotation is
related to an anisotropic eddy viscosity <cit.>.
This idea has been applied in developing a theory of the
differential rotation <cit.>.
The turbulent heat flux in these theories has been introduced
phenomenologically using the mixing-length theory relation:
⟨ u'^2 ⟩∝ g τ_0 ⟨ u'_z s' ⟩,
where ⟨ u'_z s' ⟩ is the vertical turbulent heat flux,
u' and s' are fluctuations of fluid velocity and entropy,
g is the gravity acceleration  and τ_0 is the characteristic
turbulent time. Also a quasi-linear approach that is valid for small
fluid Reynolds numbers has been applied in these studies.

Additional possibility for the production of the solar differential rotation
is associated with an effect of the turbulent heat flux on the Reynolds
stress in a rotating density stratified turbulent convection.
Based on this idea, <cit.> develop a mean field theory
of the differential rotation, where a coupled system of dynamical equations
for the Reynolds stress, the entropy fluctuations and the
turbulent heat flux has been solved adopting a spectral τ approach.
It was demonstrated  <cit.> that the ratio of the contributions
to the Reynolds stress caused by the turbulent heat flux and the anisotropic
eddy viscosity is of the order of ∼ 10   (H_ρ / ℓ_0)^2,
where ℓ_0 is the maximum scale of turbulent motions and
H_ρ is the fluid density variation scale.
This theory allows to determine the profiles of the differential rotation
in the upper part of the solar convection zone where the rotation is slow
in comparison with the turbulent time.

In the low part of the solar convective zone, the rotation is fast
in comparison with the turbulent time.
This causes a strong anisotropy of the turbulent convection
that is an additional source of the solar differential rotation.
One of the key theoretical questions is how can turbulent convection
be modified by the fast rotation, and how can it affect the production
of the differential rotation.
This issue remains to be an open unresolved problem in the solar physics
and astrophysics.
Note that different theories of the solar differential rotation can be validated
using data from the surface measurements of the solar angular velocity
<cit.> and helioseismology based on measurements
of the frequency of p-mode oscillations
<cit.>.

In the present study a combined effect of the
turbulent heat flux and the turbulence anisotropy
increasing with the rotation rate on
the Reynolds stress has been studied for a rotating
density stratified turbulent convection.
The spectral tau approach which is valid for large
Reynolds and Peclet numbers, has been used
in this study.
This allows us to advance the mean-field theory
of the solar differential rotation and obtain
the profiles of the differential rotation versus
radius which are in agreement with the
measured profiles of the solar differential rotation.



§ EFFECT OF ROTATION ON THE REYNOLDS STRESS, ENTROPY FLUCTUATIONS
AND TURBULENT HEAT FLUX


To develop the theory of differential
rotation in a small-scale density stratified turbulent convection,
we use a mean-field approach whereby the
velocity, pressure and entropy are decomposed into
mean and fluctuating parts.
This approach implies that there is a separation of temporal
and spatial scales, so that the mean fields are varied in much larger scales
in comparison with those for fluctuations.

Let us determine the dependencies of the Reynolds stresses
⟨ u'_i(t, x)   u'_j(t, x)
⟩ and the turbulent heat flux ⟨
s'(t, x)   u'_i(t, x) ⟩ on
the mean fields, where angular brackets denote the ensemble
averaging. To this end we use equations for fluctuations of
velocity and entropy in a rotating turbulent
convection, which are obtained by subtracting
equations for the mean fields from the
corresponding equations for the total
fields.
The equations for fluctuations of velocity u' and entropy s' are given by

    ∂ u' ∂ t   =    - ( U·∇)  u' - ( u' ·∇)  U -
    ∇(p' ρ_0)
    -  g  s' + 2  u' ×Ω +  U^N,
    
    ∂ s' ∂ t   =    - Ω_b^2 g (
    u' · e) - ( U·∇) s' + S^N .

Equations (<ref>) and (<ref>) are written in the reference
frame rotating with the angular velocity Ω. Here p' are fluctuations of fluid
pressure, the entropy fluctuations are
determined by s' = (γ P_0)^-1 p' -
ρ_0^-1ρ', the mean fields U and S
are the mean velocity and entropy, e is the unit vector directed opposite to g and Ω_b^2 = - 
g·∇ S. The fluid velocity for a low Mach number
flows satisfies the continuity equation written
in the anelastic approximation, div 
(ρ_0   U) = 0 and div 
(ρ_0   u') = 0. The variables with the
subscript "0" correspond to the hydrostatic
nearly isentropic basic reference state, i.e.,
∇ P_0 = ρ_0 g and g· [(γ P_0)^-1∇ P_0
- ρ_0^-1∇ρ_0] ≈
0, where γ is the ratio of specific
heats. The turbulent convection is regarded as a
small deviation from a well-mixed adiabatic
reference state.
The nonlinear terms U^N and S^N in
Eqs. (<ref>) and (<ref>) which include the molecular
dissipative terms, are given by

    U^N   =   ⟨ ( u' ·∇)  u'
    ⟩ - ( u' ·∇)  u' + 
    f_ν( u'),
    
    
    S^N   =   ⟨ ( u' ·∇) s' ⟩ - ( u' ·∇) s' - (1/
    T_0)  ∇· F_κ( u', s'),

where ρ_0   f_ν( U) is the
mean molecular viscous force, F_κ( U, S) is the mean heat flux
associated with the molecular thermal conductivity.

To study the rotating turbulent convection we
perform the derivations which include the
following steps: (i) adopting new variables for
fluctuations of velocity v = √(ρ_0)  u' and entropy s = √(ρ_0) 
s'; (ii) derivation of the equations for the
second moments of the velocity
fluctuations ⟨ v_i   v_j ⟩, the
entropy fluctuations ⟨ s^2 ⟩ and
the turbulent heat flux ⟨ v_i   s
⟩ in the k space, where we apply a multi-scale approach
<cit.>, which separates the mean fields
varied in large scales from fluctuations varied in
small scales; (iii) application of the spectral τ approximation
and solution of the derived second-moment equations in the
k space; (iv) returning to the physical
space to obtain formulae for the Reynolds
stress and the turbulent heat flux as the
functions of the rotation rate.

Using Eqs. (<ref>)-(<ref>) for the fluctuations of velocity and entropy
in k space derived in Appendix A, we obtain equations
for the following correlation functions:
f_ij( k, K) = ⟨ v_i(t, k_1) v_j(t,
k_2) ⟩,
F_i( k, K) = ⟨ s(t, k_1) v_i(t,
k_2) ⟩, and
Θ( k, K) = ⟨ s(t, k_1) s(t,
k_2) ⟩,
where k_1 =  k +  K / 2 and k_2 = - k +  K / 2.
Here the wave vectors K and k are related to the large and small scales, respectively.
Hereafter we omit the argument t in the correlation functions to simplify notations.
The equations for these second moments are given by

    ∂ f_ij( k, K) ∂ t   =   (I_ijmn^U +
    L_ijmn^Ω)   f_mn + M_ij^F + N̂f̃_ij ,
    
    ∂ F_i( k, K) ∂ t   =   (J_im^U +
    D_im^Ω)   F_m + g e_m P_im( k_1)  Θ
    + N̂F̃_i ,
    
    ∂Θ( k, K) ∂ t   =     -  div ( UΘ) + N̂Θ ,

where

    I_ijmn^U    =    J^U_im( k_1)  δ_jn + J^U_jn(
    k_2)  δ_im
    = [2 k_iqδ_mpδ_jn + 2 k_jqδ_imδ_pn - δ_imδ_jqδ_np
            - δ_iqδ_jnδ_mp + δ_imδ_jn k_q∂∂ k_p] ∇_p
    U_q - δ_imδ_jn ( div  U +  U·∇),
    
    
    M_ij^F    =    g e_m [P_im( k_1) F_j(
    k, K)  + P_jm( k_2) F_i(-
    k, K)],

and L_ijmn^Ω = D_im^Ω( k_1)  δ_jn +
D_jn^Ω( k_2)  δ_im,
J_ij^U( k) = 2 k_in∇_j U_n - ∇_j U_i -
(1/ 2)   div  U δ_ij
and D_ij^Ω( k) = 2 ε_ijmΩ_n k_mn.
Here δ_ij is the Kronecker tensor, k_ij = k_i  k_j /
k^2, ε_ijk is the Levi-Civita tensor,
and F_i(- k, K) = ⟨ s( k_2) v_i(
k_1) ⟩. The correlation functions f_ij, F_i and Θ are proportional to the non-uniform fluid density
ρ_0. Here N̂f̃_ij, N̂F̃_i and N̂Θ are the terms which are
related to the third-order moments appearing due to the nonlinear
terms. In particular,

    N̂f̃_ij   =   ⟨ P_im( k_1)
    v^N_m( k_1) v_j( k_2) ⟩
    + ⟨ v_i( k_1) P_jm( k_2) v^N_m( k_2)
    ⟩ ,
    
    N̂F̃_i   =   ⟨ s^N( k_1) u_j( k_2)
    ⟩ + ⟨ s( k_1) P_im( k_2) v^N_m(
    k_2) ⟩ ,
    
    N̂Θ   =   ⟨ s^N( k_1) s( k_2) ⟩
    + ⟨ s( k_1) s^N( k_2) ⟩ .


The equations for the second-order moments contain high-order
moments and a closure problem arises <cit.>.
We apply the spectral τ approximation
that is a sort of third-order closure procedure
<cit.>.
The spectral τ approximation postulates that the deviations of the
third-order-moment terms, N̂f_ij( k), from the
contributions to these terms afforded by the background turbulent
convection, N̂f_ij^(0)( k), are expressed
through the similar deviations of the second moments, f_ij(
k) - f_ij^(0)( k), i.e.,

    N̂f_ij( k) - N̂f_ij^(0)( k) = -
    f_ij( k) - f_ij^(0)( k) τ_r(k),

and similarly for other tensors, where N̂f_ij = N̂f̃_ij +
M_ij^F(F^Ω=0) and N̂ F_i =
N̂F̃_i + g e_n P_in(k)
Θ^Ω=0, the superscript (0)
corresponds to the background turbulent
convection (i.e., a turbulent convection with ∇_i U_j = 0), τ_r (k) is the
characteristic relaxation time of the statistical
moments, which  can be identified with the
correlation time τ(k) of the turbulent
velocity field for large Reynolds numbers. The
quantities F^Ω=0 and Θ^Ω=0
are for a nonrotating turbulent convection with
nonzero spatial derivatives of the mean velocity.
Validation of the τ approximation has been done in various numerical
simulations and analytical studies
<cit.>.
Note that we apply the
τ-approximation (<ref>) only to study the
deviations from the background turbulent
convection which are caused by the spatial
derivatives of the mean velocity.
The background turbulent convection is assumed to be known (see below).

We use the following model of the background
turbulent convection which  takes into account an
increase of the anisotropy of turbulence with
increase of the rate of rotation:

    f_ij^(0)≡⟨ v_i( k_1)  
    v_j( k_2)  ⟩^(0)   =   E(k)
      [1 + 2 k  ε_u  δ(k_z)]
    8 π  k^2   (k^2 + λ̃^2)   (1 +
    ε_u)[δ_ij  (k^2 +
    λ̃^2)  - k_i   k_j - λ̃_i
     λ̃_j
    
           + i  (λ̃_i  
    k_j - λ̃_j   k_i) ] ρ_0  ⟨ u'^2 ⟩^(0),
    
    
    F_i^(0)≡⟨ v_i( k_1)  
    s( k_2)  ⟩^(0)   =   3   E(k)  8 π k^4[k^2   e_j   P_ij(
    k) - i λ̃  k_j   P_ij( e)
    ] ρ_0  ⟨ s'   u'_z ⟩^(0) ,

and Θ^(0)≡⟨ s( k_1)   s( k_2) ⟩^(0)  =  ρ_0  ⟨ (s')^2 ⟩^(0)  E(k) / 4π k^2, where λ̃ =
(λ - ∇) / 2, λ = -(∇ρ_0) / ρ_0.
We assume that the background turbulent convection is the Kolmogorov type turbulence with
a constant flux of energy over the spectrum,
i.e., the kinetic energy spectrum
E(k) = - d τ̅(k) / dk, τ̅(k) =
(k / k_0)^1-q with the exponent of the kinetic 
energy spectrum 1 < q < 3, e.g., q =5/3 is for Kolmogorov spectrum.
The turbulent correlation time
τ(k) = 2 τ__Ωτ̅(k), where
τ__Ω = ℓ_0 / u_0, and
ℓ_0 is the energy containing scale of
turbulent motions, u_0=√(⟨ u'^2 ⟩^(0)) is the
characteristic turbulent velocity in the scale
ℓ_0 and k_0 = 1 /ℓ_0. 
We consider an anisotropic turbulent convection as a combination of a three-dimensional isotropic turbulence and two-dimensional turbulence in the plane perpendicular to the rotational axis. The degree of anisotropy ε_u is defined as  the ratio of turbulent kinetic energies of two-dimensional to three-dimensional motions.
In this model we neglect effects which are
O(λ^3, ∇^3 ⟨ v^2 ⟩^(0)).
 
The effect of rotation on the turbulent
correlation time is described just by an heuristic argument, i.e.,
we assume that
τ_Ω^-2 = τ_0^-2 + Ω^2 /  C_Ω^2,
that yields:

    τ__Ω = τ_0 [1 + (C__Ω^-1 Ω τ_0)^2]^1/2 .

This implies that for fast rotation, Ω τ_0 ≫ 1, the
parameter ω = 8 Ω τ__Ω
tends to be limiting value ω_m = 8 C__Ω, where the dimensionless constant C__Ω∼ 1.

The solution of Eqs. (<ref>)–(<ref>) after application of the spectral τ
approximation, and the integration over the k space
(see Appendix <ref>) allow
us to determine the Reynolds stress and the effective
force versus angular velocity.
The latter yields the mean-field equation for the differential rotation
(see next section), which takes into account the effects of rotating
density stratified turbulent convection.



§ MEAN-FIELD EQUATION FOR DIFFERENTIAL ROTATION


The differential rotation in the axisymmetric fluid flow is
determined by linearized Navier-Stokes equation for the toroidal component
U_φ(r, θ) ≡ r  sinθ δΩ of
the mean velocity:

    ρ_0  ∂ U_φ∂ t   =   1  r^3∂∂ r (r^3 σ_r φ) + 1  r
    sin^2 θ∂∂θ (sin^2 θ σ_θφ)
    + 2  ρ_0   ( U×Ω)_φ ,

where the tensor σ_ij = - ⟨ v_i   v_j ⟩ is
determined by the Reynolds stress:

    σ_r φ   ≡    - e_j^φ  e_i^r  ⟨ v_i  
    v_j ⟩ = σ_r φ^ν_T + σ_r φ^F + σ_r φ^u ,
    
    σ_θφ   ≡    - e_j^φ  e_i^θ ⟨ v_i   v_j ⟩ = σ_θφ^ν_T
    + σ_θφ^F + σ_θφ^u .

and e^r,  e^θ and e^φ are
the unit vectors along the radial, meridional and toroidal
directions of the spherical coordinates r, θ, φ.
There are three contributions to the tensor σ_ij = - ⟨ v_i
  v_j ⟩ in Eqs. (<ref>) and  (<ref>). The first term in the right hand
side of Eqs. (<ref>) and  (<ref>) describes the contribution
σ_ij^ν_T to the Reynolds stress caused by turbulent viscosity ν__T:

    σ_r φ^ν_T = ρ_0  ν__T  r  ∂∂ r(U_φ r),
    
    σ_θφ^ν_T = ρ_0  ν__T sinθ r ∂∂θ(U_φsinθ) .

The second term in Eqs. (<ref>) and  (<ref>) determines the contribution
σ^F to the Reynolds stress caused by the turbulent
heat flux:

    σ_r φ^F    =   1  6 ρ_0  τ__Ω^2   g  ⟨ s'   u'_z ⟩^(0) Ω sinθ  [Φ_1(ω) + cos^2 θ Φ_2(ω)]  ,
    
    σ_θφ^F    =   1  3 ρ_0  τ__Ω^2
      g  ⟨ s'   u'_z ⟩^(0) Ω sin^2 θ cosθ Φ_2(ω) ,

where the parameter ω = 8  Ω τ__Ω.
The functions Φ_1(ω) and Φ_2(ω)
are given by Eqs. (<ref>)–(<ref>) in Appendix <ref> and are shown in Fig. <ref>.
When the turbulent correlation time is independent of the rotation rate,
Eqs. (<ref>) and (<ref>) coincide with those obtained by
<cit.>.




The third term in Eqs. (<ref>) and (<ref>)
determines the contribution σ^u to the Reynolds
stress caused by the anisotropy of turbulence due to the
nonuniform fluid density and fast uniform rotation (see
Eq. (<ref>) in Appendix <ref>):

    σ_r φ^u    =    - λ^2  ℓ_0^2  20 ρ_0  ⟨ u'^2 ⟩^(0) τ__Ω Ω sinθ  (1 + cos^2 θ) ,
    
    σ_θφ^u    =   λ^2  ℓ_0^2  20 ρ_0  ⟨ u'^2 ⟩^(0) τ__Ω Ω sin^2 θ cosθ .

Equation (<ref>) in a steady-state that determines the profiles of the differential rotation, reads:

    Ŵ(r)  ρ_0 ν__T{∂∂ rδΩΩ   +   1  r [a_F  (Φ_1 (ω) + Φ_2 (ω) X^2 )
    - 2 a_u  λ^2 ℓ_0^2 (1 + X^2 ) ] }
       -   ρ_0 ν__T r^2 M̂(X)  {δΩΩ - [a_F  Φ_2 (ω) + a_u λ^2 ℓ_0^2 ]   X^2 } = 0,

where the operators Ŵ(r) and M̂(X) are defined as

    Ŵ(r) f(r)= 1  r^4 ∂∂ r [r^4 f(r)],     M̂(X) ϕ(X)= [(X^2 - 1) ∂^2 ∂ X^2 + 4 X ∂∂ X] ϕ(X),

X = cosθ, and the parameters a_F and a_u are given by a_F = τ_Ω^2   g  ⟨ s'   u'_z ⟩^(0) / 6  ν__T and
a_u = τ_Ω ⟨ u'^2 ⟩^(0) / 40  ν__T.
We seek a solution of Eq. (<ref>) in the form:

    δΩΩ = ∑_n=0^∞  C_2n^3/2(X)  Ω̃_2n(r) ,

where the radius r is measured in units of the solar radius R_⊙,
and the function C_n^3/2(X) satisfies the equation for the
ultra-spherical polynomials:

    [M̂(X) - n(n+3)]   C_n^3/2(X) = 0  .

The function C_n^3/2(X) has the following properties:

    ∫_-1^1   (1-X^2)   C_n^3/2(X)   C_m^3/2(X)  dX
    = (n+1)(n+2)  n+3/2 δ_nm ,

C_0^3/2(X) = 1 and C_2^3/2(X) = (3/2) (5X^2 - 1).
Substituting Eq. (<ref>) into Eq. (<ref>), we obtain equations for
the functions Ω̃_0(r):

    Ω̃_0(r) = Ω̃_∗ - 1  5∫^1_r/R_⊙{12 a_u λ^2 ℓ_0^2 - a_F  [5Φ_1 (ω) + Φ_2 (ω)] } dr  r,

and Ω̃_2(r):

    Ŵ(r)  ρ_0 ν__T [∂Ω̃_2(r) ∂ r + 2  15   r (a_F  Φ_2 (ω) - 2 a_u λ^2 ℓ_0^2 )]
    
                   - 10  ρ_0 ν__T r^2[Ω̃_2(r) - 2  15(a_F  Φ_2 (ω) + a_u λ^2 ℓ_0^2)]
    = 0 ,

where Ω̃_∗ is the free constant determined by the surface boundary condition.



In Fig. <ref> we show the total angular velocity Ω̃_ tot = Ω̃_0 + 1 that includes the uniform rotation Ω versus the radius r/R_⊙. This theoretical profile is compared with the radial profile of the solar angular velocity
obtained from the helioseismology observational data <cit.> specified for the latitude ϕ=30^∘ and normalized by the solar angular velocity Ω_⊙(ϕ=0) at the equator. Note that at ϕ=30^∘ the contribution from the term C_2^3/2(X)  Ω̃_2(r) to the differential rotation vanishes, because the function C_2^3/2(X) = (3/2) (5X^2 - 1) at the angle around ϕ=30^∘ vanishes.
To determine Ω̃_ tot we use the rotation rate dependence of the
turbulent viscosity ν__T(ω)=ν__T^∗Φ_ν(ω), where
ν__T^∗ = τ_0 ⟨ u'^2 ⟩^(0) / 6, the functions Φ_ν(ω) is given by Eq. (<ref>) in Appendix <ref> and is shown in Fig. <ref>. Strong change of the turbulent viscosity is caused by the fast rotation during the transition from isotropic three-dimensional turbulence to strongly anisotropic quasi two-dimensional turbulence.



For the comparison of the theoretical profiles of the differential rotation and observational data, we use the radial profiles of Ω τ_0(r) (see Fig. <ref>) and the ratio  ℓ_M(r) / H_ρ(r) (see Fig. <ref>) of the mixing length ℓ_M to the density stratification length H_ρ based on the model of the solar convective zone by <cit.>.
Inspection of Fig. <ref> demonstrates that the theoretical profile of the differential
rotation is in agreement with the profile of the solar differential rotation when
ℓ_M /ℓ_0 = 5. The latter is justified by the results of analytical study <cit.> and laboratory experiments <cit.>, which show that the integral scale ℓ_0 of the turbulent convection is smaller in 5 times in comparison with the size of the coherent structures (the large-scale circulations).
We compare the theoretical and observation profiles of the differential rotation for the latitude ϕ=30^∘ because for the latitudes which are far from ϕ=30^∘, the contribution of the term ∝Ω̃_2 (determined by Eq. (<ref>)) to the differential rotation cannot be ignored. More detail comparison of the theoretical and observation profiles of the differential rotation for different latitudes requires the mean-field numerical modelling that is a subject of the separate study.








§ CONCLUSIONS


We discuss a new theory of differential rotation
based on a combined effect of the turbulent heat flux and
the turbulence anisotropy increasing with the rate of rotation on the Reynolds stress
in a density stratified turbulent convection.
We solve a coupled system of dynamical budget equations which includes
the equations for the Reynolds stress, the entropy fluctuations
and the turbulent heat flux, applying
a spectral τ approach to close the system of these equations.
The model of the background turbulent convection takes into account an
increase of the turbulence anisotropy and a decrease of the turbulent
correlation time with the rotation rate.
This theory allows to obtain the profile of the differential
rotation versus radius which is in agreement with the
profile of the solar differential rotation.

The mechanism of the differential rotation that is related to
the effect of the turbulent heat flux on Reynolds stress in a
rotating turbulent convection is as follows. The total angular velocity
includes the uniform rotation Ω and the
differential rotation δΩ.
The uniform rotation results in the counter-rotation turbulent
heat flux ⟨ s'   u'_φ⟩ that is
directed opposite to the uniform rotation Ω.
The counter-rotation  turbulent heat flux is similar to the
counter-wind turbulent heat flux that is
directed opposite to the mean wind known in the atmospheric physics  <cit.>.
In turbulent convection an ascending fluid element obeys
larger temperature than the temperature of the surrounding fluid and smaller
toroidal fluid velocity, while a descending fluid element obeys
smaller temperature and larger toroidal fluid velocity. This results in
the turbulent heat flux in the direction opposite to the uniform rotation.
The entropy fluctuations produce fluctuations of the buoyancy force,
that increases fluctuations of the vertical and
meridional components of the velocity which are correlated with the
fluctuations of the toroidal component of the velocity. This implies that
the off-diagonal components of the Reynolds stress,
⟨ u'_r u'_φ⟩ and ⟨ u'_θ u'_φ⟩ are non-zero, producing the toroidal component of the effective force.
The latter results in the formation of the differential rotation
δΩ in turbulent convection.




This work was supported in part by the Research Council of Norway
under the FRINATEK (grant No. 231444).
The authors acknowledge the hospitality of NORDITA and Ural Federal University.







§ DERIVATION OF EQUATIONS FOR THE SECOND MOMENTS


Equations (<ref>) and (<ref>) in the new variables for
fluctuations of velocity v= √(ρ_0)  u' and entropy s = √(ρ_0) 
s' are given by

    1 √(ρ_0)∂ v( x,t) ∂ t   =    - ∇(p' ρ_0)
    + 1 √(ρ_0)[2  v×Ω -
    ( v·∇)  U
    - G^U   v -  g  s ] +  F_M +  v^N,
    
    
    ∂ s( x,t) ∂ t   =    - Ω_b^2
    g ( v· e) - G^U   s + s^N,

where G^U = (1 / 2)   div  U +  U·∇, v^N and s^N are the nonlinear
terms which include the molecular viscous and dissipative terms. The fluid velocity
fluctuations v satisfy the equation ∇·
v =  v·λ/2.

Let us derive equations for the second-order moments. For this
purpose we rewrite the momentum equation and the entropy equation in
a Fourier space. In particular,

    dv_i( k)  dt   =    [D_im^Ω( k) +
    J̃_im^U( k)] v_m( k) + g   e_m   P_im( k)  
    s( k) + v_i^N( k),
    
    ds( k)  dt   =    - G^U( k)   s( k) + s^N ,

where

    J̃_ij^U( k) = 2 k_in∇_j U_n - ∇_j U_i -
    G^U( k) δ_ij,
       G^U( k) = 1  2  div  U + i ( U· k),

D_ij^Ω( k) = 2 ε_ijmΩ_n k_mn,
δ_ij is the Kronecker tensor, k_ij = k_i  k_j /
k^2 and ε_ijk is the Levi-Civita tensor.
To derive Eq. (<ref>) we multiply the momentum equation written
in k-space by P_ij( k) = δ_ij - k_ij
to exclude the pressure term.
We also use the following identities:

    √(ρ_0) [∇× [∇× ( u' ×Ω)] ]
    = (Ω×∇^(λ))
      (λ· v)
    + (Ω·∇^(λ))  (∇^(λ)× v),
    
          √(ρ_0) [∇× [∇× ( g s')] ] = - g_j  [ δ_ij (∇^(λ))^2 -
    ∇_i^(λ)∇_j^(λ)] s,
    
          √(ρ_0) [∇× [∇× ( u)]
    ]_ k = - [Λ^2  δ_ij -
    Λ_i λ_j ] v_j( k),

where ∇^(λ) = ∇  +
λ /2, λ =
-(∇ρ_0) / ρ_0, Λ = i  k + λ /2.
Using Eqs. (<ref>) and (<ref>) we derive equations for the second moments
which are given by Eqs. (<ref>)–(<ref>).



§ SOLUTIONS FOR THE SECOND MOMENTS



Equations (<ref>)-(<ref>) in a
steady state and after  applying the spectral
τ approximation (<ref>), read

    f_ij( k)    =    L_ijmn^-1[f^(0)_mn + τ M̃_mn^F +
    τ  (I_mnpq^U + L_mnpq^∇
    + L_mnpq^λ  + L_mnpq^∇^2 +
    L_mnpq^λ^2)    f_pq],
    
    
    
    F_i( k)    =    D_im^-1[F^(0)_m( k)  + τ  (J_mn^U +
    D_mn^∇ + D_mn^λ
    + D_mn^∇^2 + D_mn^λ^2) F_n],

and Θ( k) = [1 - τ ( U·∇)]
Θ^(0)( k),
where

    M̃_ij^F    =    g e_m {[P_im(
    k) + k_im^∇  + k_im^λ -
    k_im^∇^2 + k_im^λ^2]
    F̃_j( k)
    + [P_jm( k) - k_jm^∇
           -k_jm^λ - k_jm^∇^2  + k_jm^λ^2] F̃_i(- k)
    },

and F̃_i=F_i-F_i^Ω=0 and we
neglected  terms ∼ O(∇^3, λ^3).
Here the operator D_ij^-1= χ(ψ)   (δ_ij +
ψ ε_ijm k̂_m + ψ^2   k_ij) is the inverse of
δ_ij - τD̃_ij and the
operator L_ijmn^-1(Ω) is the
inverse of δ_imδ_jn - τ L̃_ijmn <cit.>,
where

    L_ijmn^-1(Ω)    =   1  2 [B_1
     δ_imδ_jn + B_2   k_ijmn +
    B_3   (ε_impδ_jn
    + ε_jnpδ_im) k̂_p +B_4   (δ_im k_jn
           + δ_jn k_im) + B_5  ε_ipmε_jqn
    k_pq + B_6   (ε_imp k_jpn +
    ε_jnp k_ipm) ] ,

and k̂_i = k_i / k, χ(ψ) = 1 / (1
+ ψ^2), ψ = 2 τ(k)   ( k·Ω) / k, B_1 = 1 + χ(2
ψ) , B_2 = B_1 + 2 - 4 χ(ψ) , B_3 = 2 ψ χ(2 ψ) , B_4 = 2
χ(ψ) - B_1 , B_5 = 2 - B_1 and B_6
= 2 ψ  [χ(ψ) - χ(2 ψ)].

To obtain solutions for the second moments,
we extract in tensors D_ij^Ω and
L_ijmn^Ω  the parts which
depend on large-scale spatial derivatives
and on the density stratification effects:

    D_ij^Ω   =   D̃_ij + D_ij^∇
    +  D_ij^∇^2 + D_ij^λ +
    D_ij^λ^2 + O(∇^3) ,
    
    
    L_ijmn^Ω   =   L̃_ijmn+
    L_ijmn^∇  + L_ijmn^∇^2 +
    L_ijmn^λ + L_ijmn^λ^2
    + O(∇^3) ,

where

    L̃_ijmn   =    2  Ω_q  
    (ε_imp δ_jn +
    ε_jnp δ_im)   k_pq,
      
    L_ijmn^∇ = - 2 Ω_q  
    (ε_imp δ_jn -
    ε_jnp δ_im)  
    k_pq^∇ ,
    
    
    L_ijmn^λ   =    - 2 Ω_q  [(ε_imp δ_jn -
    ε_jnp δ_im)  
    k_pq^λ
    +i  k^2 (ε_ilq δ_jn λ_m- ε_jlq δ_im λ_n)   k_l],
    
    
    L_ijmn^∇^2   =    2 Ω_q  
    (ε_imp δ_jn +
    ε_jnp δ_im)  
    k_pq^∇^2,
      
    L_ijmn^λ^2 = 2 Ω_q  
    (ε_imp δ_jn +
    ε_jnp δ_im)  
    k_pq^λ^2,

and D̃_ij = 2 ε_ijpΩ_q k_pq,
D_ij^∇ = 2 ε_ijpΩ_q k_pq^∇,
D_ij^λ = 2 ε_ijpΩ_q k_pq^λ,
D_ij^∇^2 = 2 ε_ijpΩ_q k_pq^∇^2,
and D_ij^λ^2 = 2 ε_ijpΩ_q k_pq^λ^2.
Here

    k_ij^∇   =   i  2 k^2  [k_i ∇_j
    + k_j ∇_i - 2 k_ij ( k·∇)],
      
    k_ij^λ = i  2 k^2  [k_i λ_j
    + k_j λ_i - 2 k_ij ( k·λ)],
    
    
    k_ij^∇^2   =   1  4 k^2 [k_ij∇^2  + 2 (k_ip∇_j
    + k_jp∇_i)  ∇_p
    - 4 k_ijpq∇_p ∇_q - ∇_i ∇_j ],
    
    
    k_ij^λ^2   =   1  4 k^2[λ_i ∇̃_j  + λ_j
    ∇̃_i - 2 λ_m (k_im∇̃_j + k_jm∇̃_i
    + k_ij∇̃_m) + λ_i λ_j -
    k_ijλ^2  + 4 k_ijpqλ_p
    λ_q],

and ∇̃_i = ∇_i - 4 k_il∇_l.

After integration in k space we obtain
contributions to the Reynolds stress caused by
turbulence anisotropy due to the rapid rotation:

    f_ij^(u)   =   λ 20{[(ω̂× e)_i e_j +
    (ω̂× e)_j e_i ]
      (λ- ∇_z)
    + (ω̂· e)
    [(ω̂× e)_i
    ω̂_j
    
           + (ω̂× e)_j
    ω̂_i ]   (λ+ ∇_z)
    }ε_u  1+ ε_u ρ_0  ⟨ u'^2 ⟩^(0) Ωτ__Ω ℓ_0^2 .

To derive Eq. (<ref>), we use the following integrals:

    ∫ k_ij^⊥ dφ = πδ_ij^(2),
      ∫ k_ijmn^⊥ dφ = π 4Δ_ijmn^(2),

where δ_ij^(2)≡ P_ij(Ω) = δ_ij - Ω_i Ω_j /Ω^2
and Δ_ijmn^(2) = δ_ij^(2)δ_mn^(2)
+ δ_im^(2)δ_jn^(2)+ δ_in^(2)δ_jm^(2).

The contributions to the Reynolds stress caused by
the turbulent heat flux are given by Eqs. (<ref>) and (<ref>),
where the functions Φ_1(ω) and Φ_2(ω) are given by

    Φ_1(ω)    =    2 Ψ_1(ω) + Ψ_2(ω/2),   Φ_2(ω) = 2 Ψ_2(ω) + Ψ_2(ω/2),
    
    Ψ_1(ω)    =    - 6 ω^4[arctan ωω (1 + ω^2) - 8
    ω^2 3 -1
    + 2  ω  Y(ω) ],
    
    Ψ_2(ω)    =   6 ω^4[5  arctan ωω (1 + ω^2) + 8 ω^2 3 - 5
    - 6  ω  Y(ω) ] ,

ω = 8 Ωτ__Ω and Y(ω) = ∫_0^ω
[arctan  y / y]  d y.
When the turbulent correlation time is independent of the rotation rate,
Eqs. (<ref>) and (<ref>) coincide with those obtained by <cit.>.

To determine the profile of the differential rotation, we use
the rotation rate dependence of the turbulent viscosity ν__T(ω)
=ν__T^∗Φ_ν(ω), where ν__T^∗ = τ_0 ⟨
u'^2 ⟩^(0) / 6 and the functions Φ_ν(ω) is given by

    Φ_ν(ω)    =   1  8(1+ ε_u) {(q+3)  ε_u + 2  [A_1^(1)(ω) - A_1^(1)(0) + (q+ 2) C_1^(1)(0)
    + C_1^(1)(ω)]
    
           + A_2^(1)(ω) + C_3^(1)(ω)} .

Here

    A_1^(1)(ω)    =    12 [arctan (ω) ω(1 - 1 ω^2)
    + 1 ω^2[1 - ln(1 + ω^2)] ] ,
    
    
    A_2^(1)(ω)    =    - 12 [arctan (ω) ω(1 - 3 ω^2)
    + 1 ω^2[3 - 2 ln(1 + ω^2)] ]  ,
    
    
    C_1^(1)(ω)    =   arctan (ω) ω(3 -
    6 ω^2 - 1 ω^4)
    + 1 ω^2(17  3 + 1 ω^2 - 4 ln(1 + ω^2) ) ,

where C_3^(1)(ω) = A_1^(1)(ω) - 5 C_1^(1)(ω),
and we use equations derived by <cit.>, which are adopted for the spherical
geometry.



jpp



