© ESO, 2011

1. Introduction

Helio-seismological data have presented us with an interesting new feature: solar rotation changes from differential in the convective zone CZ to uniform in the radiative, stably stratified region below the CZ. The problem has been studied for many years by many authors, and in what follows we present a brief summary of the present situation:

• 1.

  3D numerical simulations by Brummel et al.(2002) and Brun & Toomre (2002) have made “promisingcontacts” with the convective zone but have not yet explained theuniform rotation in the radiative regime. The reason identified bythe authors of those studies is the limited range of physicalparameters (compared to the true solar values) allowed bynumerical simulations, e.g., the simulations are still too viscous;

• 2.

  use of the angular momentum equation with Reynolds stresses with the “standard down gradient flux” for the Reynolds stresses (shear only) leads to an extraction of angular momentum from the interior that is too weak to explain the helio data. Such models yield large rotation gradients (in the stellar interiors) that are not consistent with helio data;

• 3.

  several authors suggested that a different process is needed, internal gravity waves (IGW) which are characterized by an up-gradient flux (Kumar & Quataert 1997; Zahn et al. 1997; Kumar et al. 1999; Charbonnel & Talon 2005, 2007; Talon & Charbonnel 2005). The time scale required by the IGW mechanism to arrive at a solar rotational curve compatible with the helio data was estimated to be  ~ 10⁷ yrs;

• 4.

  we show that, if in addition to shear generated down-gradient fluxes, one includes vorticity and stable stratification, the new Reynolds stresses contain an up-gradient term that in the radiative zone (RZ), where Ω(r) becomes close to rigid, dominates the shear down-gradient term. For small Ri corresponding to a strong shear, as expected in the CZ-RZ transition zone, the time scale can be of the same order as the one provided by the IGW mechanism.

In the next sections, we first show why the angular momentum equation employed in all previous studies is an approximation of the complete one. We then derive and show how to include vorticity and stratification, which contribute with opposite signs in the CZ and in the RZ zones where they produce an up-gradient flux in the angular momentum equation.

Finally, we show how the new Reynolds stress model includes unstable stratification, stable stratification, differential rotation, double-diffusion, arbitrary Peclet number (accounting for radiative losses), and meridional currents.

2. The solar angular momentum problem

Talon & Charbonnel (2005) have summarized the successes and failures of the “hydro-dynamically induced” processes in stellar structure studies. Here, we limit ourselves to discussing the “failures” or incompleteness of such models, the most notorious of which concerns the sun since, as Thompson et al. (2003, Sect. 7.3) have written, “calculations based on the angular momentum equation:

$\frac{\mathit{\partial }}{\mathit{\partial t}}\mathrm{\left(}{\mathit{r}}^{\mathrm{2}}\mathrm{\Omega }\mathrm{\right)}\mathrm{=}{\mathit{r}}^{-2}\frac{\mathit{\partial }}{\mathit{\partial r}}\left({\mathit{r}}^{\mathrm{4}}{\mathit{K}}_{\mathrm{m}}\frac{\mathit{\partial }\mathrm{\Omega }}{\mathit{\partial r}}\right)\mathrm{+}\mathrm{·}\mathrm{·}\mathrm{·}$(1a)predict rotation of the solar interior at a rate several times higher than the surface rate, in stark disagreement with helio data of nearly uniform rotation”. Equation (1a) is the standard equation used in all studies we have consulted (e.g., Talon & Zahn 1998; Maeder & Meynet 2001; Talon & Charbonnel 2003; Mathis et al. 2004; Palacios et al. 2003, 2006). Equation (1a) is a particular case of the general equation that can be derived from the Navier-Stokes equations for the mean velocity ${\overline{\mathit{u}}}_{\mathit{i}}$: $\frac{\mathit{\partial }}{\mathit{\partial t}}\mathrm{(}\overline{\mathit{\rho }}\hspace{0.17em}{\overline{\mathit{u}}}_{\mathit{i}}\mathrm{)}\mathrm{+}\frac{\mathit{\partial }}{\mathit{\partial }{\mathit{x}}_{\mathit{j}}}\overline{\mathit{\rho }}\mathrm{(}{\mathit{R}}_{\mathit{ij}}\mathrm{+}{\overline{\mathit{u}}}_{\mathit{i}}{\overline{\mathit{u}}}_{\mathit{j}}\mathrm{)}\mathrm{=}\mathrm{-}{\mathit{g}}_{\mathit{i}}\overline{\mathit{\rho }}\mathrm{-}\frac{\mathit{\partial }\overline{\mathit{p}}}{\mathit{\partial }{\mathit{x}}_{\mathit{i}}}$(1b)where ${\mathit{R}}_{\mathit{ij}}\mathrm{=}\overline{{\mathit{u}}_{\mathit{i}}^{\mathrm{\prime }}{\mathit{u}}_{\mathit{j}}^{\mathrm{\prime }}}$ are the Reynolds stresses representing the effect of turbulence on the mean flow (the total velocity field has a mean component denoted by an overbar and a fluctuating component denoted by a prime). The ${\overline{\mathit{u}}}_{\mathit{\varphi }}\mathrm{=}\mathrm{\Gamma }\mathit{r}\mathrm{\Omega }\mathit{,}\mathrm{\Gamma }\mathrm{=}\sin \mathit{\theta }$ component of Eq. (1b) reads, with J = r²Ω, as $\frac{\mathit{\partial J}}{\mathit{\partial t}}\mathrm{=}\mathrm{-}\hspace{0.17em}{\mathrm{\Gamma }}^{-1}{\mathit{r}}^{-2}\frac{\mathit{\partial }}{\mathit{\partial r}}{\mathit{r}}^{\mathrm{3}}\mathrm{(}{\mathit{R}}_{\mathit{r\varphi }}\mathrm{+}{\overline{\mathit{u}}}_{\mathit{r}}{\overline{\mathit{u}}}_{\mathit{\varphi }}\mathrm{)}\mathrm{-}{\mathrm{\Gamma }}^{-3}\frac{\mathit{\partial }}{\mathit{\partial \theta }}{\mathrm{\Gamma }}^{\mathrm{2}}\mathrm{(}{\mathit{R}}_{\mathit{\theta \varphi }}\mathrm{+}{\overline{\mathit{u}}}_{\mathit{\theta }}{\overline{\mathit{u}}}_{\mathit{\varphi }}\mathrm{)}\mathit{.}$(1c)Since the equations for the meridional currents ${\overline{\mathit{u}}}_{\mathit{r}}$, ${\overline{\mathit{u}}}_{\mathit{\theta }}$, given in Sect. 6, also depend on the Reynolds stresses, we concentrate on the latter. Thus far, all models have assumed that R_ij is contributed only by shear S_ij, that is, they adopted the following model: Down-gradient type: ${\mathit{R}}_{\mathit{ij}}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\mathrm{-}\hspace{0.17em}\mathrm{2}{\mathit{K}}_{\mathrm{m}}{\mathit{S}}_{\mathit{ij}}$(1d)where $\mathrm{2}{\mathit{S}}_{\mathit{ij}}\mathrm{=}\hspace{0.17em}{\overline{\mathit{u}}}_{\mathit{i,j}}\mathrm{+}{\overline{\mathit{u}}}_{\mathit{j,i}}\mathit{,} {\mathit{S}}_{\mathit{r\varphi }}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}\mathit{r}\mathrm{\Gamma }\frac{\mathit{\partial }\mathrm{\Omega }}{\mathit{\partial r}}$(1e)where K_m > 0 is a momentum diffusivity. Use of Eqs. (1d), (1e) in (1c) and integration over angles yield (1a). There are several problems with (1d). The first is the need to justify why only shear enters, and second is how one determines the momentum diffusivity. The first item is almost never discussed and the determination of K_m is handled with heuristic models discussed in Sect. 5 of Paper II, that is, most of the attention was devoted to the determination of K_m and very little, if any, to the completeness of the “shear alone” model (1d). This leads to (1a), which is not a diffusion equation, in spite of being generally referred to as such. In fact, we show that inclusion of vorticity and stratification adds to (1a) a new term that has a truly diffusive character: $\frac{\mathit{\partial J}}{\mathit{\partial t}}\mathrm{=}\mathrm{-}{\mathit{r}}^{-2}\frac{\mathit{\partial }}{\mathit{\partial r}}\left({\mathit{r}}^{\mathrm{2}}\mathit{x}{\mathit{K}}_{\mathrm{m}}\frac{\mathit{\partial J}}{\mathit{\partial r}}\right)\mathrm{+}\mathrm{·}\mathrm{·}\mathrm{·}\hspace{0.17em}\mathit{,} \mathit{x}\mathrm{=}\mathrm{(}\mathit{\tau N}{\mathrm{)}}^{\mathrm{2}}\mathit{.}$(1f)Here, N is the Brunt-Väisäla frequency ${\mathit{N}}^{\mathrm{2}}\mathrm{=}\mathrm{-}\mathit{g}{\mathit{\rho }}_{\mathrm{0}}^{-1}\mathit{\partial }\overline{\mathit{\rho }}\mathit{/}\mathit{\partial z}$ and τ = 2K / ε (K is the eddy kinetic energy and ε its rate of dissipation) is the dynamical time scale that will be determined later. The key point is that in the CZ where x < 0, Eq. (1f) is still of the down-gradient type like (1a), but in the RZ where x > 0, Eq. (1f) becomes of the up-gradient type¹: $\frac{\mathit{\partial J}}{\mathit{\partial t}}\mathrm{=}\mathrm{-}{\mathit{r}}^{-2}\frac{\mathit{\partial }}{\mathit{\partial r}}\left({\mathit{r}}^{\mathrm{2}}\left|\mathit{x}\right|{\mathit{K}}_{\mathrm{m}}\frac{\mathit{\partial J}}{\mathit{\partial r}}\right)\mathrm{+}\mathrm{·}\mathrm{·}\mathrm{·}$(1g) This means that the generally assumed down-gradient nature of shear instabilities is incorrect at least in the context of stellar interiors where the transition from unstable to stable stratification changes the nature of the instabilities making it up-gradient. The new feature is entirely due to the combination of vorticity+buoyancy, as Eq. (1g) shows. An additional, interesting feature of (1g) is that it does not vanish for Ω = const. characterizing the radiative zone, whereas in the same regime, Eq. (1a) has a zero rhs. In conclusion, as one approaches the radiative zone, the dominant term is the non-vanishing up-gradient term (1g).

In the next sections, we show that shear alone giving rise to (1a) is not a justifiable approximation, and then we proceed to extend (1d) to include vorticity, buoyancy, radiative losses, double diffusion, and meridional currents.

3. Why shear alone?

Consider the first relation in Eq. (1d). It represents the first term in a Taylor expansion in the parameter τ / T where τ is the dynamical time scale while T is a typical time scale characterizing the mean flow. This can be seen by rewriting the first of (1d) as follows: ${\mathit{R}}_{\mathit{ij}}\mathit{/}\mathit{K}\mathrm{~}\mathrm{(}{\mathit{K}}_{\mathrm{m}}\mathit{/}\mathit{TK}\mathrm{)}\mathrm{\left(}\mathit{T}{\mathit{S}}_{\mathit{ij}}\mathrm{\right)}\mathrm{~}\frac{\mathit{\tau }}{\mathit{T}}\mathrm{\left(}\mathit{T}{\mathit{S}}_{\mathit{ij}}\mathrm{\right)}$(2a)where we have employed the following relations: $\frac{{\mathit{K}}_{\mathrm{m}}}{\mathit{TK}}\mathrm{~}\frac{\mathit{\tau }}{\mathit{T}}\mathit{,} {\mathit{K}}_{\mathrm{m}}\mathrm{~}\mathit{\tau K}\mathit{.}$(2b)There is no physical reason why τ / T should be small, which would justify stopping at the first term (2a), and non-perturbative derivations of the Reynolds stresses (Taulbee 1992; Gatski & Speziale 1993) showed there are higher order terms in powers of τ / T. A second consideration is that (1d) requires that the principal axes of the tensor τ_ij representing turbulence be aligned with those of S_ij representing the mean flow. This is true for the case of pure strain but not for flows with a mean vorticity, which is defined as ${\mathit{V}}_{\mathit{ij}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}\mathrm{(}{\overline{\mathit{u}}}_{\mathit{i,j}}\mathrm{-}{\overline{\mathit{u}}}_{\mathit{j,i}}\mathrm{)}\mathit{,} {\mathit{V}}_{\mathit{r\varphi }}\mathrm{=}\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}{\mathit{r}}^{-1}\mathrm{\Gamma }\frac{\mathit{\partial }{\mathit{r}}^{\mathrm{2}}\mathrm{\Omega }}{\mathit{\partial r}}$(2c)which is not zero even when shear is zero since shear and vorticity are two independent, orthogonal tensors, that is, when Ω →  rigid body, S_rϕ → 0,   V_rϕ ≠ 0. For a 3D flow, in general, the measured flow distribution can only be predicted by choosing different viscosities for each stress component (Markatos 1987). In fact, a complete derivation of (1d) shows the presence of non-isotropic terms that break the “alignment assumption” and that are ultimately responsible for the extra terms discussed above. In other words, in lieu of the first of (1d), one has (Taulbee 1992; Gatski & Speziale 1993) an expression of the type ${\mathit{b}}_{\mathit{ij}}\mathrm{=}{\mathit{R}}_{\mathit{ij}}\mathrm{-}\mathrm{2}\mathit{K}\mathit{/}\mathrm{3}{\mathit{\delta }}_{\mathit{ij}}\mathrm{:} {\mathit{b}}_{\mathit{ij}}\mathrm{~}{\mathit{a}}_{\mathrm{1}}{\mathit{S}}_{\mathit{ij}}\mathrm{+}{\mathit{a}}_{\mathrm{2}}{\mathit{b}}_{\mathit{ik}}{\mathit{S}}_{\mathit{kj}}\mathrm{+}{\mathit{a}}_{\mathrm{3}}{\mathit{b}}_{\mathit{ik}}{\mathit{V}}_{\mathit{kj}}\mathrm{+}\mathrm{·}\mathrm{·}\mathrm{·}$(2d)of which the form (8h) of Paper I is the specific form after the closure constants have been determined (a ≈ 0,b = c = 1 / 10), a relation that, for completeness, we repeat here: ${\mathit{b}}_{\mathit{ij}}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\mathrm{-}\frac{\mathrm{8}}{\mathrm{75}}\mathit{K\tau }{\mathit{S}}_{\mathit{ij}}\mathrm{-}\mathit{a\tau }{\mathrm{\Sigma }}_{\mathit{ij}}\mathrm{-}\mathit{b\tau }{\mathit{Z}}_{\mathit{ij}}\mathrm{+}\mathit{c\tau }{\mathit{B}}_{\mathit{ij}}\mathit{.}$(2e)The traceless tensors Σ_ij,Z_ij,B_ij are defined as$\begin{array}{ccc}& & \mathrm{Shear:} {\mathrm{\Sigma }}_{\mathit{ij}}\mathrm{=}{\mathit{b}}_{\mathit{ik}}{\mathit{S}}_{\mathit{kj}}\mathrm{+}{\mathit{b}}_{\mathit{jk}}{\mathit{S}}_{\mathit{ik}}\mathrm{-}\frac{\mathrm{2}}{\mathrm{3}}{\mathit{\delta }}_{\mathit{ij}}{\mathit{b}}_{\mathit{km}}{\mathit{S}}_{\mathit{km}}\\ & & \mathrm{Vorticity:} {\mathit{Z}}_{\mathit{ij}}\mathrm{=}{\mathit{b}}_{\mathit{ik}}{\mathit{V}}_{\mathit{jk}}\mathrm{+}{\mathit{b}}_{\mathit{jk}}{\mathit{V}}_{\mathit{ik}}\\ & & \mathrm{Buoyancy:} {\mathit{B}}_{\mathit{ij}}\mathrm{=}\mathit{g}\left({\mathit{\lambda }}_{\mathit{i}}{\mathit{J}}_{\mathit{j}}^{\mathit{\rho }}\mathrm{+}{\mathit{\lambda }}_{\mathit{j}}{\mathit{J}}_{\mathit{i}}^{\mathit{\rho }}\mathrm{-}\frac{\mathrm{2}}{\mathrm{3}}{\mathit{\delta }}_{\mathit{ij}}{\mathit{\lambda }}_{\mathit{k}}{\mathit{J}}_{\mathit{k}}^{\mathit{\rho }}\right)\end{array}$where ${\mathit{J}}_{\mathit{i}}^{\mathit{\rho }}$ is the buoyancy flux, where Eqs. (7f), in Paper I comprises both heat and μ fluxes and where ${\mathit{\lambda }}_{\mathit{i}}\mathrm{\equiv }\mathrm{-}\hspace{0.17em}\mathrm{(}\mathit{g}\overline{\mathit{\rho }}{\mathrm{)}}^{-1}\mathit{\partial }\overline{\mathit{p}}\mathit{/}\mathit{\partial }{\mathit{x}}_{\mathit{i}}$. Thus, in going from (1d) to (2e), we have gone from ${\mathit{R}}_{\mathit{ij}}\mathrm{\left(}\mathit{S}\mathrm{\right)}\mathrm{\to }\hspace{0.17em}{\mathit{R}}_{\mathit{ij}}\mathrm{\left(}\mathit{S,V,B}\mathrm{\right)}$(2i)where S, V, B stand for shear, vorticity and buoyancy. Given the incompleteness of the lhs of (2i), the failures of (1a) referred to by Thompson et al. (2003) may look less surprising, since in reality (1a) should contain shear, vorticity, and buoyancy. Other considerations are also in order. First, as we have already alluded to, one would expect the rhs of (1a) to exhibit a diffusion form while (1a) does not, though it is generally, but improperly, referred to as a “diffusion equation”. Since diffusion involves small scales, it is clear that shear, which is governed by large scales, cannot represent diffusion. Since small scales have large vorticity, it is only natural that the latter be present in (2i). Indeed, one may notice that the second of (2c) contains the angular momentum r²Ω, whereas shear, Eq. (1e), does not. Therefore, both the physical interpretation of vorticity dominated by small-scale processes (of diffusive nature) and its mathematical structure, indicate that its presence in the angular momentum equation will give rise to “diffusion”. Second, to properly describe the qualitative difference in the rotation curves in the convective and radiative regimes, it is only natural to have the buoyancy flux, which is positive in the first and negative in the second regime. The presence of B in (2i) is therefore not only natural but physically required. Third, what about the Peclet number representing radiative losses? If we recall Eq. (13j) of Paper I, it is clear that Pe is large in the CZ and small in the radiative region. This is because both the eddy velocity and the length scale are large in the CZ but small in the stably stratified radiative regime. Thus, Eq. (2i) must be extended to read as ${\mathit{R}}_{\mathit{ij}}\mathrm{\left(}\mathit{S,V,B}|\mathrm{Pe}\mathrm{\right)}\mathit{.}$(2j)Finally, what about meridional currents? As shown above, they enter directly into the angular momentum Eqs. (1b), (1c) and thus we further generalize (2j) to the form ${\mathit{R}}_{\mathit{ij}}\mathrm{\left(}\mathit{S,V,B}|\mathrm{Pe}\mathit{,M}\mathrm{\right)}$(2k) where M stands for meridional currents. If we succeed in constructing (2k), we would have included shear, vorticity, different regimes of both unstable (convective zone, CZ) and stable stratification (radiative regime, RZ), radiative losses, and meridional currents. While there is no guarantee that the resulting rotational curve will explain the helio data, we would have at least made sure that we have included the key processes that characterize the two regimes of interest, CZ and RZ. The RSM (Reynolds stress model) that we have presented in Paper I, is the model we employ next to compute the complete Reynolds stress R_ij.

4. Previous results of the RSM model

It is fair to ask whether the RSM has been previously employed in a stellar context. In answer to this question, we point to the work of Kupka and collaborators (reviewed in details in Kupka & Muthsam 2007; Kupka 2009), who studied turbulent convection using an RSM model in which the stresses and heat flux depend on time, but without rotation. Use of (2k) with rotation was made under the following conditions and with the following results. Under unstably stratified conditions, without meridional currents and using the observed rotational curve (Ulrich et al. 1988), $\mathrm{\Omega }\mathrm{\left(}\mathit{\theta }\mathrm{\right)}\mathrm{=}{\mathrm{\Omega }}_{\mathrm{0}}\mathrm{(}\mathrm{1}\mathrm{-}\mathit{a}{\cos }^{\mathrm{2}}\mathit{\theta }\hspace{0.17em}\mathrm{-}\hspace{0.17em}\mathit{b}{\cos }^{\mathrm{4}}{\mathit{\theta }}^{\mathrm{)}}$(3a)the model results for the Reynolds stresses ${\mathit{R}}_{\mathit{\theta \varphi }}\mathrm{=}\overline{{\mathit{u}}_{\mathit{\theta }}^{\mathrm{\prime }}{\mathit{u}}_{\mathit{\varphi }}^{\mathrm{\prime }}} \mathrm{vs}\mathit{.} \mathit{\theta ,}$(3b)were compared with existing solar surface measurements (Virtanen 1989; Pulkkinen et al. 1993). None of the models with $\mathrm{shear} \mathrm{alone}\mathit{,} \mathrm{shear}\hspace{0.17em}\mathrm{+}\hspace{0.17em}\mathrm{buoyancy}\mathit{,} \mathrm{shear}\hspace{0.17em}\mathrm{+}\hspace{0.17em}\mathrm{vorticity}$(3c)were able to reproduce the measured surface data of (3b). In particular, the first two combinations in (3c) gave the wrong sign in both hemispheres. Only the combination $\mathrm{vorticity} \mathrm{+} \mathrm{buoyancy}$(3d)was able to reproduce the data (Canuto et al. 1994), a conclusion that should be viewed as a useful hint when describing regimes removed from the surface. However, as already discussed, thus far all the solutions of the angular momentum equation, which yielded results in disagreement with helio data, were based on the same assumption, the first of (3c), which failed to reproduce the surface values of (3b).

5. Reynolds stresses, heat, and concentration fluxes

The specific form of the equation for the Reynolds stresses has the structure (2e), and its complete form is (see Paper I, Eq. (8h)) ${\mathit{R}}_{\mathit{ij}}\mathrm{=}\frac{\mathrm{2}\mathit{K}}{\mathrm{3}}{\mathit{\delta }}_{\mathit{ij}}\mathrm{-}\frac{\mathrm{8}}{\mathrm{75}}\mathit{K\tau }{\mathit{S}}_{\mathit{ij}}\mathrm{-}\frac{\mathrm{1}}{\mathrm{10}}\mathit{\tau }{\mathit{Z}}_{\mathit{ij}}\mathrm{+}\frac{\mathrm{1}}{\mathrm{10}}\mathit{\tau }{\mathit{B}}_{\mathit{ij}}\mathit{.}$(4a)In (4a), the vorticity and buoyancy tensors were defined in Eqs. (2f) − (2h), while shear and vorticity S_ij and V_ij were defined in Eqs. (1e) and (2c). The density fluxes (heat and concentration fluxes) appearing in (2h) are defined as ${\mathit{J}}_{\mathit{i}}^{\mathit{\rho }}\mathrm{=}{\mathit{\alpha }}_{\mathit{T}}{\mathit{J}}_{\mathit{i}}^{\mathit{h}}\mathrm{-}{\mathit{\alpha }}_{\mathit{c}}{\mathit{J}}_{\mathit{i}}^{\mathit{c}}\mathit{.}$(4b)The heat ${\mathit{J}}_{\mathit{i}}^{\mathit{h}}$ and concentration fluxes ${\mathit{J}}_{\mathit{i}}^{\mathit{c}}$ are given by the following algebraic equations (Paper I, Eqs. (10a), (10b) and (12d), (12e)):

Heat fluxes: $\begin{array}{ccc}& & \mathrm{(}{\mathit{\delta }}_{\mathit{ij}}\mathrm{+}{\mathit{\mu }}_{\mathit{ij}}\mathrm{)}{\mathit{J}}_{\mathit{j}}^{\mathit{h}}\mathrm{=}{\mathit{\gamma }}_{\mathit{ij}}{\mathit{\beta }}_{\mathit{j}}\\ & & {\mathit{\gamma }}_{\mathit{ij}}\mathrm{=}{\mathit{\pi }}_{\mathrm{4}}\mathit{\tau }\mathrm{(}{\mathit{R}}_{\mathit{ij}}\mathrm{-}\mathit{g}{\mathit{\alpha }}_{\mathit{c}}{\mathit{\pi }}_{\mathrm{2}}\mathit{\tau }{\mathit{\lambda }}_{\mathit{i}}{\mathit{J}}_{\mathit{j}}^{\mathit{c}}\mathrm{)}\\ & & \end{array}$(5a)Concentration fluxes: $\begin{array}{ccc}& & \mathrm{(}{\mathit{\delta }}_{\mathit{ij}}\mathrm{+}{\mathit{\eta }}_{\mathit{ij}}\mathrm{)}{\mathit{J}}_{\mathit{j}}^{\mathit{c}}\mathrm{=}\mathrm{-}\hspace{0.17em}{\mathit{d}}_{\mathit{ij}}{\overline{\mathit{C}}}_{\mathit{,j}}\\ & & {\mathit{d}}_{\mathit{ij}}\mathrm{=}{\mathit{\pi }}_{\mathrm{1}}\mathit{\tau }\mathrm{(}{\mathit{R}}_{\mathit{ij}}\mathrm{+}\mathit{g}{\mathit{\alpha }}_{\mathit{T}}{\mathit{\pi }}_{\mathrm{2}}\mathit{\tau }{\mathit{\lambda }}_{\mathit{i}}{\mathit{J}}_{\mathit{j}}^{\mathit{h}}\mathrm{)}\\ & & \end{array}$(5b)\arraycolsep1.75ptwhere ${\mathit{\beta }}_{\mathit{i}}\mathrm{=}\mathrm{-}\frac{\mathit{\partial T}}{\mathit{\partial }{\mathit{x}}_{\mathit{i}}}\mathrm{-}{\mathit{\lambda }}_{\mathit{i}}\mathit{g}{\mathit{c}}_{\mathit{p}}^{-1}\mathit{,} {\mathit{\lambda }}_{\mathit{i}}\mathrm{=}\mathrm{-}\mathrm{(}\mathit{g}\overline{\mathit{\rho }}{\mathrm{)}}^{-1}\frac{\mathit{\partial }\overline{\mathit{p}}}{\mathit{\partial }{\mathit{x}}_{\mathit{i}}}\mathrm{·}$(5c)The dimensionless time scales represented by the functions π′s defined in Paper I Eq. (10c), are given by: $\begin{array}{ccc}{\mathit{\pi }}_{\mathrm{1}}& \mathrm{=}& \hspace{0.17em}{\mathit{\pi }}_{\mathrm{1}}^{\mathrm{0}}{\left(\mathrm{1}\mathrm{+}\frac{\mathit{Ri}{\mathit{R}}_{\mathit{\mu }}}{\mathit{a}\mathrm{+}{\mathit{R}}_{\mathit{\mu }}}\right)}^{-1}\mathit{,} {\mathit{\pi }}_{\mathrm{4}}\mathrm{=}\hspace{0.17em}{\mathit{\pi }}_{\mathrm{4}}^{\mathrm{0}}\mathit{f}\mathrm{\left(}\mathrm{Pe}\mathrm{\right)}{\left(\mathrm{1}\mathrm{+}\frac{\mathit{Ri}}{\mathrm{1}\mathrm{+}\mathit{a}{\mathit{R}}_{\mathit{\mu }}}\right)}^{-1}\\ {\mathit{\pi }}_{\mathrm{2}}& \mathrm{=}& {\mathit{\pi }}_{\mathrm{2}}^{\mathrm{0}}\mathrm{(}\mathrm{1}\mathrm{+}\mathit{Ri}{\mathrm{)}}^{-1}\mathrm{[}\mathrm{1}\mathrm{+}\mathrm{2}\mathit{Ri}{\mathit{R}}_{\mathit{\mu }}\mathrm{(}\mathrm{1}\mathrm{+}{\mathit{R}}_{\mathit{\mu }}^{\mathrm{2}}{\mathrm{)}}^{-1}\mathrm{]}\mathit{,}\hspace{0.17em}{\mathit{\pi }}_{\mathrm{5}}\mathrm{=}{\mathit{\pi }}_{\mathrm{5}}^{\mathrm{0}}\mathit{g}\mathrm{\left(}\mathrm{Pe}\mathrm{\right)}\mathit{,}\\ {\mathit{\pi }}_{\mathrm{1}}^{\mathrm{0}}& \mathrm{=}& {\mathit{\pi }}_{\mathrm{4}}^{\mathrm{0}}\mathrm{=}\mathrm{(}\mathrm{27}\mathit{K}{\mathit{o}}^{\mathrm{3}}\mathit{/}\mathrm{5}{\mathrm{)}}^{\mathrm{-}\mathrm{1}\mathit{/}\mathrm{2}}\mathrm{(}\mathrm{1}\mathrm{+}{\mathit{\sigma }}_{\mathit{t}}^{-1}{\mathrm{)}}^{-1}\mathit{,}\\ {\mathit{\pi }}_{\mathrm{2}}^{\mathrm{0}}& \mathrm{=}& \mathrm{1}\mathit{/}\mathrm{3}\mathit{,}{\mathit{\pi }}_{\mathrm{3}}\mathrm{=}{\mathit{\pi }}_{\mathrm{3}}^{\mathrm{0}}\mathrm{=}{\mathit{\pi }}_{\mathrm{5}}^{\mathrm{0}}\mathrm{=}{\mathit{\sigma }}_{\mathit{t}}\\ \mathit{f}\mathrm{\left(}\mathrm{Pe}\mathrm{\right)}& \mathrm{=}& \mathit{b}\mathrm{Pe}\mathrm{(}\mathrm{1}\mathrm{+}\mathit{b}\mathrm{Pe}{\mathrm{)}}^{-1}\mathit{,}\\ \mathit{g}\mathrm{\left(}\mathrm{Pe}\mathrm{\right)}& \mathrm{=}& \end{array}$(5d)In addition, the Peclet number Pe that quantifies radiative losses is given by (Canuto & Dubovikov 1998; Paper II, Eq. (5b), with τ = 2K / ε) $\mathrm{Pe}\mathrm{=}\frac{{\mathit{\pi }}^{\mathrm{2}}}{\mathrm{125}}{\mathrm{Re}}_{\mathit{\chi }}\mathrm{(}\mathit{\tau N}{\mathrm{)}}^{\mathrm{2}}\mathit{,} {\mathrm{Re}}_{\mathit{\chi }}\mathrm{=}\frac{\mathit{\epsilon }}{\mathit{\chi }{\mathit{N}}^{\mathrm{2}}}$(5e)where Re_χ is the Reynolds number based on the thermometric diffusivity (χis in cm² s^-1; Pr = ν / χ ≈ 10^-8 is the Prandtl number). The Richardson number Ri, the density ratio R_μ, the Brunt-Väisäla frequency, and the mean shear are defined as $\mathit{Ri}\mathrm{=}\frac{{\mathit{N}}^{\mathrm{2}}}{{\mathrm{\Sigma }}^{\mathrm{2}}}\mathit{,} {\mathit{R}}_{\mathit{\mu }}\mathrm{\equiv }\frac{\mathit{g}{\mathit{\alpha }}_{\mathit{c}}{\overline{\mathit{C}}}_{\mathit{,z}}}{\mathit{g}{\mathit{\alpha }}_{\mathit{T}}{\overline{\mathit{T}}}_{\mathit{,z}}}\mathrm{=}\frac{{\mathrm{\nabla }}_{\mathit{\mu }}}{\mathrm{\nabla }\mathrm{-}{\mathrm{\nabla }}_{\mathrm{ad}}}$(6a)${\mathit{N}}^{\mathrm{2}}\mathrm{=}\mathrm{-}\mathit{g}{\mathit{H}}_{\mathit{p}}^{-1}\mathrm{(}\mathrm{\nabla }\mathrm{-}{\mathrm{\nabla }}_{\mathrm{ad}}\mathrm{)}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{R}}_{\mathit{\mu }}\mathrm{)}\mathit{,} \mathrm{\Sigma }\mathrm{=}\mathrm{(}\mathrm{2}{\mathit{S}}_{\mathit{ij}}{\mathit{S}}_{\mathit{ij}}{\mathrm{)}}^{\mathrm{1}\mathit{/}\mathrm{2}}\mathit{.}$(6b)In summary, the above formalism, which is entirely algebraic, contains the following physical ingredients: $\begin{array}{ccc}& & \mathrm{Differential} \mathrm{rotation}\\ & & \mathrm{Unstable}\mathrm{/}\mathrm{Stable} \mathrm{stratification}\\ & & \mathrm{DoubleDiffusion}\\ & & \mathrm{RadiativelossesarbitraryPecletnumber}\\ & & \mathrm{Meridionalcurrents}\mathrm{.}\end{array}$(6c)Once Eqs. (4a) and (5a), (5b) are solved in spherical coordinates, the resulting Reynolds stresses have the general form (2k) and can then be substituted in (1c) whose solution may then be compared with the helio data. Clearly, such a computation must be done in conjunction with a solar structure code to provide the mean variables. One can actually perform two computations, with and without double-diffusion, and compare the results, a process that would be quite instructive.

5.1. the variables K, τ

The above equations still depend on two turbulence variables represented by τ = 2Kε^-1 and K, which often enter together, a product that can be written as $\mathit{\tau K}\mathrm{=}\frac{\mathrm{2}{\mathit{K}}^{\mathrm{2}}}{\mathit{\epsilon }}\mathrm{=}\mathrm{4}{\mathit{\ell }}^{\mathrm{2}}{\mathit{\tau }}^{-1}\mathrm{=}\mathrm{2}\mathit{\ell }{\mathit{K}}^{\mathrm{1}\mathit{/}\mathrm{2}}\mathrm{=}\mathrm{2}{\mathit{\epsilon }}^{\mathrm{1}\mathit{/}\mathrm{3}}{\mathit{\ell }}^{\mathrm{4}\mathit{/}\mathrm{3}}\mathit{.}$(7a)These are all equivalent expressions highlighting different variables; for example, the last expression is the well known 1926 Richardson law discovered fifteen years before the advent of the Kolmogorov law. In (7a) use was made of the relation ε = ℓ^-1K^3 / 2, where ℓ represents a typical eddy size. No matter which of the representations (7a) one chooses, the basic fact is that one must specify two independent variables, K − ε or, alternatively, K − τ. As discussed in Sect. 3 of Paper II, one must in principle solve two differential equations for the two variables K − ε given by Eqs. (II, (8a), (8b)). However, this is seldom done because of the added complexity. Instead of the equation for K, one uses its local limit corresponding to assuming production = dissipation, that is, ${\mathit{P}}_{\mathit{s}}\mathrm{+}{\mathit{P}}_{\mathit{b}}\mathrm{=}\mathit{\epsilon ,}$(7b)where the shear and buoyancy components are given by ${\mathit{P}}_{\mathit{s}}\mathrm{=}\mathrm{-}{\mathit{b}}_{\mathit{ij}}{\mathit{S}}_{\mathit{ij}}\mathrm{=}\mathrm{-}{\mathit{R}}_{\mathit{ij}}{\mathit{S}}_{\mathit{ij}}\mathit{>}\mathrm{0}\mathit{,} {\mathit{P}}_{\mathit{b}}\mathrm{=}\hspace{0.17em}\mathit{g}{\mathit{\lambda }}_{\mathit{i}}{\mathit{J}}_{\mathit{i}}^{\mathit{\rho }}\mathit{,}$(7c) which provides the time scale τ as a function of the large-scale variables such as shear, vorticity and temperature gradients. See for example Eqs. (9d), (9e) of Paper II. The other variable ε is discussed in Sect. 3 Paper II.

6. Meridional currents

In addition to Eq. (1c) for ${\overline{\mathit{u}}}_{\mathit{\varphi }}\mathrm{=}\mathit{r}\mathrm{\Omega \Gamma }$, for completeness we also write down the model independent general dynamic equations for the meridional currents ${\overline{\mathit{u}}}_{\mathit{r}}\mathit{,}{\overline{\mathit{u}}}_{\mathit{\theta }}$: $\begin{array}{ccc}\frac{\mathit{\partial }}{\mathit{\partial t}}\mathrm{\left(}\mathit{\rho }{\overline{\mathit{u}}}_{\mathit{r}}\mathrm{\right)}& \mathrm{=}& {\mathit{A}}_{\mathrm{1}}\mathrm{-}\frac{\mathit{\partial P}}{\mathit{\partial r}}\mathrm{+}\mathit{\rho }\frac{\mathit{\partial }\mathrm{\Phi }}{\mathit{\partial r}}\\ \frac{\mathit{\partial }}{\mathit{\partial t}}\mathrm{\left(}\mathit{\rho }{\overline{\mathit{u}}}_{\mathit{\theta }}\mathrm{\right)}& \mathrm{=}& {\mathit{A}}_{\mathrm{2}}\mathrm{-}{\mathit{r}}^{-1}\frac{\mathit{\partial P}}{\mathit{\partial \theta }}\mathrm{+}\mathit{\rho }{\mathit{r}}^{-1}\frac{\mathit{\partial }\mathrm{\Phi }}{\mathit{\partial \theta }}\mathit{,}\end{array}$where $\begin{array}{ccc}{\mathit{A}}_{\mathrm{1}}& \mathrm{=}& \mathrm{-}{\mathit{r}}^{-2}\frac{\mathit{\partial }}{\mathit{\partial r}}\mathrm{\left(}{\mathit{r}}^{\mathrm{3}}{\mathit{\psi }}_{\mathit{rr}}\mathrm{\right)}\mathrm{-}\mathrm{\left(}\mathit{r}\mathrm{\Gamma }{\mathrm{)}}^{-1}\frac{\mathit{\partial }}{\mathit{\partial \theta }}\mathrm{\right(}\mathrm{\Gamma }{\mathit{\psi }}_{\mathit{r\theta }}\mathrm{)}\mathrm{+}{\mathit{r}}^{-1}\mathrm{(}{\mathit{\psi }}_{\mathit{\theta \theta }}\mathrm{+}{\mathit{\psi }}_{\mathit{\varphi \varphi }}\mathrm{)}\\ {\mathit{A}}_{\mathrm{2}}& \mathrm{=}& \mathrm{-}{\mathit{r}}^{-3}\frac{\mathit{\partial }}{\mathit{\partial r}}\mathrm{\left(}{\mathit{r}}^{\mathrm{3}}{\mathit{\psi }}_{\mathit{rr}}\mathrm{\right)}\hspace{0.17em}\mathrm{-}\mathrm{\left(}\mathit{r}\mathrm{\Gamma }{\mathrm{)}}^{-1}\frac{\mathit{\partial }}{\mathit{\partial \theta }}\mathrm{\right(}\mathrm{\Gamma }{\mathit{\psi }}_{\mathit{\theta \theta }}\mathrm{)}\mathrm{+}\mathrm{(}\mathit{rtg\theta }{\mathrm{)}}^{-1}{\mathit{\psi }}_{\mathit{\varphi \varphi }}\\ {\mathit{\psi }}_{\mathit{ij}}& \mathrm{=}& \mathit{\rho }\mathrm{(}{\mathit{R}}_{\mathit{ij}}\mathrm{+}{\overline{\mathit{u}}}_{\mathit{i}}{\overline{\mathit{u}}}_{\mathit{j}}\mathrm{)}\end{array}$where P is the mean pressure and Φ is gravitational potential. Equations (8) show that the meridional currents depend on the Reynolds stresses as well.

7. Example: no meridional currents – no double diffusion

In Canuto & Minotti (2001), we presented the explicit solution of Eqs. (4a) and (5a) including meridional currents. The solutions represented by nested algebraic relations, were suited to a numerical treatment but not to physical considerations. Here, we present a solution of the same equations with no meridional currents and no double diffusion, a simplification that allows highlighting a new feature of the model, the existence of a counter-gradient angular momentum flux within the hydrodynamic instability framework. We begin by introducing the following dimensionless variables: $\begin{array}{ccc}{\mathit{r}}_{\mathit{ij}}& \mathrm{=}& \frac{{\mathit{R}}_{\mathit{ij}}}{\mathit{K}}\mathit{,} {\mathit{s}}_{\mathit{ij}}\mathrm{=}\mathit{\tau }{\mathit{S}}_{\mathit{ij}}\mathit{,} {\mathit{v}}_{\mathit{ij}}\mathrm{=}\mathit{\tau }{\mathit{V}}_{\mathit{ij}}\mathit{,} {\mathit{z}}_{\mathit{ij}}\mathrm{=}\frac{\mathit{\tau }{\mathit{Z}}_{\mathit{ij}}}{\mathit{K}}\mathit{,}\\ {\mathit{b}}_{\mathit{ij}}& \mathrm{=}& \end{array}$(9a)The dimensionless form of Eq. (4a) then reads ${\mathit{r}}_{\mathit{ij}}\mathrm{=}\frac{\mathrm{2}}{\mathrm{3}}{\mathit{\delta }}_{\mathit{ij}}\mathrm{-}\frac{\mathrm{8}}{\mathrm{75}}{\mathit{s}}_{\mathit{ij}}\mathrm{-}\frac{\mathrm{1}}{\mathrm{10}}{\mathit{z}}_{\mathit{ij}}\mathrm{+}\frac{\mathrm{1}}{\mathrm{10}}{\mathit{b}}_{\mathit{ij}}\mathit{,}$(9b)where $\begin{array}{ccc}{\mathit{z}}_{\mathit{ij}}& \mathrm{=}& \left({\mathit{r}}_{\mathit{ik}}\mathrm{-}\frac{\mathrm{2}}{\mathrm{3}}{\mathit{\delta }}_{\mathit{ik}}\right){\mathit{v}}_{\mathit{jk}}\mathrm{+}\left({\mathit{r}}_{\mathit{jk}}\mathrm{-}\frac{\mathrm{2}}{\mathrm{3}}{\mathit{\delta }}_{\mathit{jk}}\right){\mathit{v}}_{\mathit{ik}}\\ {\mathit{b}}_{\mathit{ij}}& \mathrm{=}& {\mathit{\lambda }}_{\mathit{i}}{\mathit{j}}_{\mathit{j}}\mathrm{+}{\mathit{\lambda }}_{\mathit{j}}{\mathit{j}}_{\mathit{i}}\mathrm{-}\frac{\mathrm{2}}{\mathrm{3}}{\mathit{\delta }}_{\mathit{ij}}{\mathit{\lambda }}_{\mathit{k}}{\mathit{j}}_{\mathit{k}}\mathit{.}\end{array}$The dimensionless form of Eq. (5a) reads as$\begin{array}{ccc}& & \mathrm{(}{\mathit{\delta }}_{\mathit{ik}}\mathrm{+}{\mathit{\mu }}_{\mathit{ik}}\mathrm{)}{\mathit{j}}_{\mathit{k}}\mathrm{=}{\mathit{\pi }}_{\mathrm{4}}{\mathit{r}}_{\mathit{ik}}{\mathit{n}}_{\mathit{k}}\\ & & \end{array}$(9e)Equations (9b) − (9e) form a system of linear, coupled, algebraic equations that may be solved using a method of symbolic algebra; additionally, we must solve the relation P = ε which, using (9a), reads as: $\mathrm{-}\mathrm{(}{\mathit{r}}_{\mathit{\varphi r}}{\mathit{s}}_{\mathit{\varphi r}}\mathrm{+}{\mathit{r}}_{\mathit{\theta \varphi }}{\mathit{s}}_{\mathit{\theta \varphi }}\mathrm{)}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}{\mathit{j}}_{\mathit{r}}\mathrm{=}\mathrm{1.}$(9f)For illustrative purposes, we take ${\mathit{\lambda }}_{\mathit{r}}\mathrm{=}\mathrm{1}\mathit{,}\hspace{0.17em}{\mathit{\lambda }}_{\mathit{\theta ,j}}\mathrm{=}\mathrm{0}\mathit{,}\hspace{0.17em}{\mathit{n}}_{\mathit{i}}\mathrm{=}\mathrm{\left(}{\mathit{n}}_{\mathit{r}}\mathit{,}\mathrm{0}\mathit{,}\mathrm{0}\mathrm{\right)}\mathit{,}\hspace{0.17em}\mathit{x}\mathrm{=}{\mathit{\tau }}^{\mathrm{2}}{\mathit{N}}^{\mathrm{2}}\mathit{,}{\mathit{N}}^{\mathrm{2}}\mathrm{=}\mathrm{-}\mathit{g}{\mathit{\rho }}_{\mathrm{0}}^{-1}{\overline{\mathit{\rho }}}_{\mathit{z}}\mathit{.}$(10a)Solving (9b) − (9e), the resulting Reynolds stress and heat fluxes have the forms: ${\mathit{R}}_{{\mathit{r}}_{\mathit{\varphi }}}\mathrm{=}\hspace{0.17em}\mathrm{-}\hspace{0.17em}\mathrm{2}{\mathit{\nu }}_{\mathit{t}}\mathrm{(}{\mathit{S}}_{\mathit{m}}^{\mathrm{\left(}\mathrm{1}\mathrm{\right)}}{\mathit{S}}_{\mathit{r\varphi }}\hspace{0.17em}\mathrm{+}\hspace{0.17em}\mathit{x}{\mathit{S}}_{\mathit{m}}^{\mathrm{\left(}\mathrm{2}\mathrm{\right)}}{\mathit{V}}_{\mathit{r\varphi }}\mathrm{)}\mathit{,} {\mathit{J}}_{\mathit{r}}\mathrm{=}\mathrm{-}{\mathit{\nu }}_{\mathit{t}}{\mathit{S}}_{\mathit{h}}\frac{\mathit{\partial T}}{\mathit{\partial r}}\mathit{,}$(10b)where ${\mathit{\nu }}_{\mathit{t}}\mathrm{=}\hspace{0.17em}\mathit{\tau K}\mathrm{=}\frac{\mathrm{2}{\mathit{K}}^{\mathrm{2}}}{\mathit{\epsilon }}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}{\mathit{\tau }}^{\mathrm{2}}\mathit{\epsilon }\mathit{.}$(10c)Here, S_m,h are dimensionless structure functions that depend on the parameters of the problem and that are obtained once Eq. (9b) is solved, as plotted in Fig. 2. Using the relations u_ϕ = rΓΩ(r,θ), u_r = u_θ = 0, Γ = sinθ, we further have ${\mathit{S}}_{\mathit{\varphi r}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}\mathit{r}\mathrm{\Gamma }\frac{\mathit{\partial }\mathrm{\Omega }}{\mathit{\partial r}}\mathrm{=}\mathrm{-}\frac{\mathit{p}}{\mathrm{2}}\frac{\mathit{\omega }}{\mathit{\tau }}\mathrm{\Gamma }\mathit{,} {\mathit{V}}_{\mathit{r\varphi }}\mathrm{=}\mathrm{-}\frac{\mathrm{\Gamma }}{\mathrm{2}\mathit{r}}\frac{\mathit{\partial }\mathrm{\left(}{\mathit{r}}^{\mathrm{2}}\mathrm{\Omega }\mathrm{\right)}}{\mathit{\partial r}}\mathrm{=}\mathrm{-}\left(\mathrm{1}\mathrm{-}\frac{\mathit{p}}{\mathrm{2}}\right)\frac{\mathit{\omega }}{\mathit{\tau }}\mathrm{\Gamma }$(10d)with S_ϕr + V_ϕr =  − ΓΩ, and we introduced the dimensionless variables: $\mathit{p}\mathrm{=}\mathrm{-}\frac{\mathit{r}}{\mathrm{\Omega }}\frac{\mathit{\partial }\mathrm{\Omega }}{\mathit{\partial r}}\mathit{,} \mathit{\omega }\mathrm{=}\mathit{\tau }\hspace{0.17em}\mathrm{\Omega }\mathit{,} \mathit{Ri}\mathrm{=}\frac{\mathrm{1}}{{\mathit{p}}^{\mathrm{2}}{\mathrm{\Gamma }}^{\mathrm{2}}}\frac{{\mathit{N}}^{\mathrm{2}}}{\hspace{0.17em}{\mathrm{\Omega }}^{\mathrm{2}}}\mathrm{·}$(10e)Solving (9f), one obtains the functions $\mathit{\omega }\mathrm{\left(}\mathit{Ri,p}\mathrm{\right)}\mathit{,} {\mathit{S}}_{\mathit{m,h}}\mathrm{\left(}\mathit{Ri,p}\mathrm{\right)}\mathit{.}$(10f)The angular momentum Eq. (1c) then becomes: $\frac{\mathit{\partial }}{\mathit{\partial t}}\mathrm{\left(}{\mathit{r}}^{\mathrm{2}}\mathrm{\Omega }\mathrm{\right)}\mathrm{=}{\mathit{r}}^{-2}\frac{\mathit{\partial }}{\mathit{\partial r}}\left({\mathit{\nu }}_{\mathit{t}}{\mathit{r}}^{\mathrm{4}}{\mathit{S}}_{\mathit{m}}^{\mathrm{\left(}\mathrm{1}\mathrm{\right)}}\frac{\mathit{\partial }\mathrm{\Omega }}{\mathit{\partial r}}\right)\mathrm{-}{\mathit{r}}^{-2}\frac{\mathit{\partial }}{\mathit{\partial r}}\left[{\mathit{\nu }}_{\mathit{t}}\mathit{x}{\mathit{S}}_{\mathit{m}}^{\mathrm{\left(}\mathrm{2}\mathrm{\right)}}{\mathit{r}}^{\mathrm{2}}\frac{\mathit{\partial }{\mathit{r}}^{\mathrm{2}}\mathrm{\Omega }}{\mathit{\partial r}}\right]\mathrm{·}$(11)The first term has the same form as in (1d), while the second term includes the contribution of vorticity and stratification x. Several comments are in order:

• 1.

  Vorticity and buoyancy appear together in (10b), which issomething of a surprise since such a combination was not obviousin the starting Eqs. (9b) − (9e).

• 2.

  Since x represents stratification, we have x < 0 in the CZ (convective zone) and x > 0 in the RZ (radiative zone).

• 3.

  The first term in the rhs of (11) is independent of stratification and does not have the form of a diffusion of J = r²Ω, as we have already pointed out in the discussion after Eq. (1f). The second term in (11) is of the diffusion type, and its sign depends on whether one is in the CZ x < 0 or in the RZ zone, x > 0. There is an alternative representation of (11) with the form $\begin{array}{ccc}\frac{\mathit{\partial }}{\mathit{\partial t}}\mathrm{\left(}{\mathit{r}}^{\mathrm{2}}\mathrm{\Omega }\mathrm{\right)}& \mathrm{=}& {\mathit{r}}^{-2}\frac{\mathit{\partial }}{\mathit{\partial r}}\left({\mathrm{\Gamma }}_{\mathit{m}}^{\mathrm{\left(}\mathrm{1}\mathrm{\right)}}\frac{\mathit{\epsilon }}{{\mathrm{\Omega }}^{\mathrm{2}}}{\mathit{r}}^{\mathrm{4}}\frac{\mathit{\partial }\mathrm{\Omega }}{\mathit{\partial r}}\right)\\ & & \mathrm{-}{\mathit{r}}^{-2}\frac{\mathit{\partial }}{\mathit{\partial r}}\left[{\mathrm{\Gamma }}_{\mathit{m}}^{\mathrm{\left(}\mathrm{2}\mathrm{\right)}}\frac{\mathit{\epsilon }}{{\mathrm{\Omega }}^{\mathrm{2}}}{\mathit{r}}^{\mathrm{2}}\frac{\mathit{\partial }{\mathit{r}}^{\mathrm{2}}\mathrm{\Omega }}{\mathit{\partial r}}\right]\mathit{,}\end{array}$(12a)where we have defined the two dimensionless variables: ${\mathrm{\Gamma }}_{\mathit{m}}^{\mathrm{\left(}\mathrm{1}\mathrm{\right)}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}{\mathit{\omega }}^{\mathrm{2}}{\mathit{S}}_{\mathit{m}}^{\mathrm{\left(}\mathrm{1}\mathrm{\right)}}\mathit{,} {\mathrm{\Gamma }}_{\mathit{m}}^{\mathrm{\left(}\mathrm{2}\mathrm{\right)}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}\mathit{x}{\mathit{\omega }}^{\mathrm{2}}{\mathit{S}}_{\mathit{m}}^{\mathrm{\left(}\mathrm{2}\mathrm{\right)}}\mathit{.}$(12b)

$\frac{\mathit{\partial }\mathrm{\Omega }}{\mathit{\partial r}}$

$\frac{\mathit{\partial }}{\mathit{\partial t}}\mathrm{\left(}{\mathit{r}}^{\mathrm{2}}\mathrm{\Omega }\mathrm{\right)}\mathrm{=}{\mathit{r}}^{-2}\frac{\mathit{\partial }}{\mathit{\partial r}}\left({\mathrm{\Gamma }}_{\mathit{m}}\frac{\mathit{\epsilon }}{{\mathrm{\Omega }}^{\mathrm{2}}}{\mathit{r}}^{\mathrm{4}}\frac{\mathit{\partial }\mathrm{\Omega }}{\mathit{\partial r}}\right)\mathrm{-}\mathrm{2}{\mathit{r}}^{-2}\frac{\mathit{\partial }}{\mathit{\partial r}}\left({\mathrm{\Gamma }}_{\mathit{m}}^{\mathrm{\left(}\mathrm{2}\mathrm{\right)}}\frac{\mathit{\epsilon }}{{\mathrm{\Omega }}^{\mathrm{2}}}{\mathit{r}}^{\mathrm{3}}\mathrm{\Omega }\right)\mathit{,}$(13)

${\mathrm{\Gamma }}_{\mathit{m}}\mathrm{=} {\mathrm{\Gamma }}_{\mathit{m}}^{\mathrm{\left(}\mathrm{1}\mathrm{\right)}}\mathrm{-}{\mathrm{\Gamma }}_{\mathit{m}}^{\mathrm{\left(}\mathrm{2}\mathrm{\right)}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}{\mathit{\omega }}^{\mathrm{2}}\mathrm{(}{\mathit{S}}_{\mathit{m}}^{\mathrm{\left(}\mathrm{1}\mathrm{\right)}}\mathrm{-}\mathit{x}{{\mathit{S}}_{\mathit{m}}^{\mathrm{\left(}\mathrm{2}\mathrm{\right)}}}^{\mathrm{)}}\mathit{>}\mathrm{0.}$(14)

Fig. 1 The dimensionless structure functions Γm defined in Eqs. (12b) and (14) vs. Ri. Ri < 0 corresponds to the CZ, while Ri > 0 corresponds to the RZ. For comments, see the text.

Fig. 2 The dimensionless functions Sm in Eq. (10b).

7.1. unstable-stable stratification

In the CZ, x = τ²N² < 0, Ri < 0 and since ${\mathrm{\Gamma }}_{\mathit{m}}^{\mathrm{\left(}\mathrm{2}\mathrm{\right)}}\mathit{<}\mathrm{0}$, the last term in (13) has a positive sign. The first term also has a positive sign, and thus the two terms on the rhs in (13) both have a positive sign, which means that they are both of the counter-gradient type. The value of ε is large since the convective flux ${\mathit{F}}_{\mathit{c}}\mathrm{=}\mathit{\rho }{\mathit{c}}_{\mathit{p}}\overline{\mathit{w\theta }}\mathrm{=}\mathrm{(}\mathit{\rho }{\mathit{c}}_{\mathit{p}}\mathit{/}\mathit{g\alpha }\mathrm{)}\mathit{g\alpha }\overline{\mathit{w\theta }}\mathrm{=}\mathrm{(}\mathit{\rho }{\mathit{c}}_{\mathit{p}}\mathit{/}\mathit{g\alpha }\mathrm{)}\mathit{\epsilon }$ carries almost all the solar flux. Owing to the p dependence exhibited in relation (10d), the value p = 2 corresponds to the case of pure shear (no vorticity).

In the RZ, x = τ²N² > 0, Ri > 0, we have ${\mathrm{\Gamma }}_{\mathit{m}}^{\mathrm{\left(}\mathrm{2}\mathrm{\right)}}\mathit{>}\mathrm{0}\mathit{,}{\mathrm{\Gamma }}_{\mathit{m}}\mathit{>}\mathrm{0}$. The first term in (13) is of the counter-gradient type, while the second term is up-gradient. At the bottom of the CZ, ε can be related to the power generated by internal gravity waves (Kumar et al. 1999), as already discussed in Eq. (10a) of Paper II. The term ${\mathrm{\Gamma }}_{\mathit{m}}^{\mathrm{\left(}\mathrm{2}\mathrm{\right)}}$ increases as the rotation curve becomes increasingly flatter (from p = 2 to p = 1), and so does Γ_m in the first term but its increase is largely cancelled by the decrease in the term ∂Ω / ∂r, which in principle becomes negligible as one approaches the RZ where a rigid body rotation sets in.

7.2. comparison with previous models

It is instructive to compare the momentum diffusivity K_m obtained in this work with the one used in the literature. It is usually denoted by D_v (Charbonnel & Talon 2005; Talon & Charbonnel 2005, Eq. (5); Zahn 2008, Eq. (3.5)) $\mathit{Ri}\mathit{<}\mathit{Ri}\mathrm{\left(}\mathit{cr}\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}}{\mathrm{4}}\mathrm{:} \frac{{\mathit{K}}_{\mathrm{m}}}{\mathit{\chi }}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}\mathit{Ri}}\mathit{,} \mathit{Ri}\mathrm{=}\frac{\mathrm{1}}{{\mathit{p}}^{\mathrm{2}}{\mathrm{\Gamma }}^{\mathrm{2}}}\frac{{\mathit{N}}^{\mathrm{2}}}{\hspace{0.17em}{\mathrm{\Omega }}^{\mathrm{2}}}\mathit{,}$(15)where χ(cm²   s^-1) is thermometric conductivity that enters the radiative diffusivity K_r = c_pρχ (see Eq. (4d) of Paper I). Several comments about (15) are in order.

The first is about the absence of a factor representing the amount of energy (or power) that creates the mixing. In fact, turbulence is a process that does not generate or destroy energy, rather, it distributes whatever energy is put into the system among a wide variety of scales. Without such an energy, there would be no turbulent motion or, in the presence of turbulence, turning such source of energy off would lead to a decaying turbulent mixing. The well known Kolmogorov law E(k) = Koε^2 / 3k ^− 5 / 3 gives the spectrum of the eddies generated by the nonlinear interaction, but it contains the rate of dissipation or energy input ε which is considered an outside, given variable.

The second comment is about the lack of universality of ε since different stellar interiors have different ε, which then yield different rates of mixing. Thus, the presence of such a factor would ensure that different stars give rise to different states of mixing, as ought to be the case. In the present formalism, one sees from (13) that ε is present and that the dimensionless factor, ${\mathrm{\Gamma }}_{\mathit{m}}\frac{\mathit{\epsilon }}{\mathit{\chi }{\mathrm{\Omega }}^{\mathrm{2}}}\mathit{,}$(16a)is a measure how much stronger the turbulent diffusivity is than the radiative one represented by χ. The factor (16a) is clearly different for different stars. Third, in (15), the momentum diffusivity is assumed to be a decreasing function of Ri. On the other hand, in the present formalism, the momentum diffusivity can be written in a variety of forms beginning with Eq. (10b) and using (11): $\begin{array}{ccc}\frac{{\mathit{K}}_{\mathrm{m}}^{\mathrm{\left(}\mathrm{1}\mathrm{\right)}}}{\mathit{\chi }}& \mathrm{=}& \frac{{\mathit{\nu }}_{\mathit{t}}}{\mathit{\chi }}{\mathit{S}}_{\mathit{m}}^{\mathrm{\left(}\mathrm{1}\mathrm{\right)}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}{\mathit{\tau }}^{\mathrm{2}}\frac{\mathit{\epsilon }}{\mathit{\chi }}{\mathit{S}}_{\mathit{m}}^{\mathrm{\left(}\mathrm{1}\mathrm{\right)}}\mathrm{=}\frac{\mathit{\epsilon }}{\mathit{\chi }{\mathrm{\Omega }}^{\mathrm{2}}}{\mathrm{\Gamma }}_{\mathit{m}}^{\mathrm{\left(}\mathrm{1}\mathrm{\right)}}\\ \frac{{\mathit{K}}_{\mathrm{m}}^{\mathrm{\left(}\mathrm{2}\mathrm{\right)}}}{\mathit{\chi }}& \mathrm{=}& \frac{{\mathit{\nu }}_{\mathit{t}}}{\mathit{\chi }}{\mathit{S}}_{\mathit{m}}^{\mathrm{\left(}\mathrm{2}\mathrm{\right)}}\mathrm{=}\frac{\mathit{\epsilon }}{\mathit{\chi }{\mathrm{\Omega }}^{\mathrm{2}}}{\mathrm{\Gamma }}_{\mathit{m}}^{\mathrm{\left(}\mathrm{2}\mathrm{\right)}}\mathit{.}\end{array}$(16b) Inspection of Figs. 1 shows the basic differences between (15) and (16): a) the function Γ_m does not decrease with Ri, whereas (15) does and b) (15) is limited to values of Ri < 1 / 4, whereas (16b) embraces any value of Ri. As already discussed in Sect. 7 of Paper I and in Sect. 5 of Paper II, relation (15) assumes the existence of a critical Ri above which there is no mixing, while the most recent data show that such a limit does not exist.

7.3. internal gravity waves

As discussed by the authors cited previously (Kumar & Quataert 1997; Zahn et al. 1997; Kumar et al. 1999; Charbonnel & Talon 2005, 2007), the most recently employed form of the angular momentum equation (Talon & Charbonnel 2005, Eq. (9)), has an additional term due to IGW in the righthand side of Eq. (1f) $\mathrm{-}\frac{\mathrm{3}}{\mathrm{8}\mathit{\pi \rho }}\frac{\mathrm{1}}{{\mathit{r}}^{\mathrm{2}}}\frac{\mathit{\partial }}{\mathit{\partial r}}{\mathit{L}}_{\mathrm{IGW}}\mathrm{\left(}\mathit{r}\mathrm{\right)}$(17a)where the explicit form of the IGW luminosity L_IGW can be found in the reference above. Its magnitude is approximately 10²⁹ erg s^-1 corresponding to 0.004% of the luminosity at the base of the CZ of 2.5 × 10³³ erg s^-1 (Talon & Charbonnel 2005, Sect. 3.3). From the new equation for the angular momentum, one obtains a time scale of $\mathit{\tau }\mathrm{\left(}\mathrm{IGW}\mathrm{\right)}\mathrm{=}\frac{\mathrm{8}\mathit{\pi }}{\mathrm{3}}\frac{\mathit{\rho }{\mathit{r}}^{\mathrm{5}}\mathrm{\Omega }}{{\mathit{L}}_{\mathrm{igw}}}\mathrm{\approx }\mathrm{a}\hspace{0.17em}\mathrm{few}\hspace{0.17em}{\mathrm{10}}^{\mathrm{7}} \mathrm{yrs}\mathit{,}$(17b)where we used r = 3.5 × 10⁹   cm,ρ = 125   gr   cm^-3,Ω = 2 × 10^-6   s^-1,L_IGW = 10²⁹   erg   s^-1. The value (17b) agrees with the results obtained from more detailed computations discussed in the cited literature. By comparison, if we use the standard Eq. (1a) corresponding to a down-gradient model (DG), we obtain $\mathit{\tau }\mathrm{\left(}\mathrm{DG}\mathrm{\right)}\mathrm{=}\frac{{\mathit{r}}^{\mathrm{2}}}{{\mathit{K}}_{\mathrm{m}}}\mathrm{\approx }\frac{\mathrm{4}\mathrm{\times }{\mathrm{10}}^{\mathrm{4}}}{{\mathit{K}}_{\mathrm{m}}}{\mathrm{10}}^{\mathrm{7}} \mathrm{yrs}\mathit{.}$(17c)Since the values of K_m given in Fig. 14 of Talon & Charbonnel (2005), where K_m is called D_v at say the bottom of the CZ, are of the order of 50 cm² s^-1, the τ resulting from (17c) is two-to-three orders of magnitude larger than (17b). Finally, let us take the last term in (13) corresponding to an up-gradient model (UG). We obtain $\mathit{\tau }\mathrm{\left(}\mathrm{UP}\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}\frac{{\mathit{r}}^{\mathrm{2}}{\mathrm{\Omega }}^{\mathrm{2}}}{\mathit{\epsilon }{\mathrm{\Gamma }}_{\mathit{m}}^{\mathrm{\left(}\mathrm{2}\mathrm{\right)}}}\mathrm{\approx }\frac{\mathrm{1}}{\mathit{\epsilon }{\mathrm{\Gamma }}_{\mathit{m}}^{\mathrm{\left(}\mathrm{2}\mathrm{\right)}}} \mathrm{yrs}\mathrm{\approx }\frac{\mathrm{1}}{\mathrm{300}{\mathrm{\Gamma }}_{\mathrm{m}}^{\mathrm{\left(}\mathrm{2}\mathrm{\right)}}}\mathrm{\times }{\mathrm{10}}^{\mathrm{7}} \mathrm{yrs}\mathit{,}$(17d)where we have taken ε ≈ 3 × 10^-5   cm²   s^-3 that corresponds to L_IGW = 10²⁹   erg   s^-1. Relation (17d) is compatible with the IGW result (17b) if $\mathrm{300}{\mathrm{\Gamma }}_{\mathit{m}}^{\mathrm{\left(}\mathrm{2}\mathrm{\right)}}⩽\mathrm{1.}$(17e) Figures 1c − d show that this is possible.

8. Conclusions

Given the well known challenge of trying to reproduce the helio data on the solar rotation curve, we have examined the ingredients of the angular momentum equation. It contains two main terms, the Reynolds stresses and the meridional currents which, in a large Re regime such as that characterizing a stellar interior and in a steady state, must balance each other. It is worth noticing that this balance has not yet been exhibited by the numerical simulations published thus far (see Figs. 11 of Brun & Toomre 2002), most probably because they are still too viscous. Since the equations for the meridional currents also depend on the Reynolds stresses, the latter constitute a key ingredient and we have therefore concentrated on how they have been modeled thus far and what the missing terms are that must be included. The final new formula for the Reynolds stress is quite simple, Eq. (4a) and the hope is that it will be tested and assessed to ascertain which angular momentum profiles it produces and what improvement it brings with respect to the expression used thus far, Eq. (1d).

A key feature of the new RSM is that all the relevant equations governing the Reynolds stresses (4a) and the heat and concentration fluxes (5a,b) are obtained by solving linear algebraic equations. This is a welcome feature if one considers the large amount of information that the new model contains: stable stratification, unstable stratification, double-diffusion, differential rotation, shear, radiative losses (arbitrary Peclet number) and meridional currents.

As an illustrative example we have worked out the case of no double diffusion and no meridional currents so as to highlight a key feature of the model. The standard RSM model, based on shear alone, is of the down-gradient type, and it fails to reproduce the helio data that point to a rigid body rotation below the solar convective zone. It was then suggested that IGW, which operate on an up-gradient flux, may be responsible for such a rotational state, and quantitative computations by several authors have confirmed that on a time scale of the order of 10⁷ yrs, the rigid body rotation can be achieved.

Here, we have introduced an alternative that we hasten to stress is not ad hoc but is part and parcel of the RSM model. It is unavoidable and it turns out that the combination vorticity+stable stratification gives rise to an up-gradient term in the Reynolds stresses, which, in the angular momentum equation, produces a characteristic time scale comparable to the one by the IGW model. As Talon & Charbonnel (2005) point out, a mixing model must do more than just reproduce the solar rotation curve and the model presented here must thus await those tests to be performed before a final judgment can be made.

References

All Figures

	Fig. 1 The dimensionless structure functions Γm defined in Eqs. (12b) and (14) vs. Ri. Ri < 0 corresponds to the CZ, while Ri > 0 corresponds to the RZ. For comments, see the text.
In the text

	Fig. 2 The dimensionless functions Sm in Eq. (10b).
In the text