© ESO, 2011

1. Formulation of the problem

Outside the unstably stratified CZ (convective cones) of stellar interiors there are dynamically important stably stratified regions referred to as OV (overshooting regions), which over the years have attracted a great deal of attention, and yet the true extent of such regions is still uncertain. In the solar case (Roxburgh 1978, 1989; Anderson et al. 1990; Spiegel & Zahn 1992), helio-seismological data suggest an OV extent of 0.05–0.1 H_p (Basu & Antia 1994; Basu et al. 1994; Thompson et al. 2003; Christensen-Dalsgaard et al. 2010), while theoretical predictions from both models (Xiong 1986; Deng & Xiong 2008) and numerical simulations (Brummet et al. 2002) yield results that are an order of magnitude greater. The similarity of the latter results may actually be coincidental since it is known that numerical simulations are still too viscous; for example, the Prandtl number is O(1) while the real value of O(10^-8) and thus the numerical simulations do not yet resolve the full spectrum of eddies, while theoretical models are based on the assumption of infinite Reynolds numbers and thus account for all the eddies (for recent review of both methods, see Kupka 2008). This means not only that caution must be exercised in comparing the results of these different approaches but also that, should helio data be used for comparison, it must be kept in mind that they are derived under the assumption that the temperature gradient in the OV region is adiabatic. In our opinion, before a meaningful comparison is possible between; a) helio data; b) numerical simulations; and c) theoretical OV models, the last two must be improved internally, and in this paper we discuss method c) and suggest a reformulation of the procedure used thus far with the intent of overcoming the existing shortcomings.

The procedure used to compute the OV in method c), within the context of stellar structure-evolutionary codes, is to solve the stellar structure equations, together with a turbulence model (of variable complexity), where a key component is the equation for the eddy kinetic energy K. That is the formulation that we discuss. In general, the equation for K has the following form, see Paper I, Eqs. (7c)–(7f): $\frac{\mathit{\partial K}}{\mathit{\partial t}}\mathrm{+}\underset{\mathrm{large}\hspace{0.17em}\mathrm{scale}\hspace{0.17em}\mathrm{velocities}}{\underset{\mathrm{􏽼}0\mathrm{􏽻􏽺}0\mathrm{􏽽}}{{\overline{\mathit{u}}}_{\mathit{j}}\frac{\mathit{\partial K}}{\mathit{\partial }{\mathit{x}}_{\mathit{j}}}}}\mathrm{+}\underset{\mathrm{non}\mathrm{-}\mathrm{local}\hspace{0.17em}\mathrm{term}}{\underset{\mathrm{􏽼}0\mathrm{􏽻􏽺}0\mathrm{􏽽}}{\frac{\mathit{\partial }{\mathit{F}}_{\mathit{i}}^{\mathrm{ke}}}{\mathit{\partial }{\mathit{x}}_{\mathit{i}}}}}\mathrm{=}\mathrm{-}\begin{array}{}\underset{\mathrm{􏽼}0\mathrm{􏽻􏽺}0\mathrm{􏽽}}{{\mathit{R}}_{\mathit{ij}}{\mathit{S}}_{\mathit{ij}}}\\ \mathrm{shear}\end{array}\mathrm{+}\underset{\mathrm{buoyancy}\hspace{0.17em}\mathrm{flux}}{\underset{\mathrm{􏽼}0\mathrm{􏽻􏽺}0\mathrm{􏽽}}{\mathit{g}{\mathit{\lambda }}_{\mathit{i}}{\mathit{J}}_{\mathit{i}}^{\mathit{\rho }}}}\mathrm{-}\underset{\mathrm{dissipation}}{\underset{\mathrm{􏽼}0\mathrm{􏽻􏽺}0\mathrm{􏽽}}{\mathit{\epsilon }}}$(1a)where shear S_ij, Reynolds stresses R_ij, buoyancy flux ${\mathit{J}}_{\mathit{i}}^{\mathit{\rho }}$, and flux of kinetic energy ${\mathit{F}}_{\mathit{i}}^{\mathrm{ke}}$ are defined as $\begin{array}{ccc}& & \mathrm{2}{\mathit{S}}_{\mathit{ij}}\mathrm{=}{\overline{\mathit{u}}}_{\mathit{i,j}}\mathrm{+}{\overline{\mathit{u}}}_{\mathit{j,i}}\mathit{,} {\overline{\mathit{u}}}_{\mathit{i,j}}\mathrm{=}\mathit{\partial }{\overline{\mathit{u}}}_{\mathit{i}}\mathit{/}\mathit{\partial }{\mathit{x}}_{\mathit{j}}\mathit{,} {\mathit{R}}_{\mathit{ij}}\mathrm{=}\overline{{\mathit{u}}_{\mathit{i}}^{\mathrm{\prime }}{\mathit{u}}_{\mathit{j}}^{\mathrm{\prime }}}\mathit{,}\\ & & {\mathit{J}}_{\mathit{i}}^{\mathit{\rho }}\mathrm{=}\mathrm{-}{\overline{\mathit{\rho }}}^{-1}{\overline{{\mathit{\rho }}^{\mathrm{\prime }}{\mathit{u}}^{\mathrm{\prime }}}}_{\mathit{i}}\mathrm{=}{\mathit{\alpha }}_{\mathit{T}}{\mathit{J}}_{\mathit{i}}^{\mathit{h}}\mathrm{-}{\mathit{\alpha }}_{\mathit{c}}{\mathit{J}}_{\mathit{i}}^{\mathit{c}}\\ & & {\mathit{F}}_{\mathit{i}}^{\mathrm{ke}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}\overline{{\mathit{q}}^{\mathrm{2}}{\mathit{u}}_{\mathit{i}}^{\mathrm{\prime }}}\mathit{,} {\mathit{u}}_{\mathit{j}}^{\mathrm{\prime }}{\mathit{u}}_{\mathit{j}}^{\mathrm{\prime }}\mathrm{=}{\mathit{q}}^{\mathrm{2}}\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{·}\end{array}$(1b)Here, ${\mathit{J}}_{\mathit{i}}^{\mathit{h}}\mathrm{=}\overline{{\mathit{u}}_{\mathit{i}}^{\mathrm{\prime }}{\mathit{T}}^{\mathrm{\prime }}}\mathrm{and}\hspace{0.17em}\hspace{0.17em}{\mathit{J}}_{\mathit{i}}^{\mathit{c}}\mathrm{=}\overline{{\mathit{u}}_{\mathit{i}}^{\mathrm{\prime }}{\mathit{c}}^{\mathrm{\prime }}}$ represent the heat and concentration fluxes (c is an arbitrary concentration field and α_T =  − ρ^-1(∂ρ / ∂T)_p,c,   α_c =  + ρ^-1(∂ρ / ∂c)_p,T are the expansion coefficients). Several comments are in order concerning Eq. (1a). First, the second term on the lhs represents the advection of K by azimuthal and meridional currents ${\overline{\mathit{u}}}_{\mathit{j}}$ whose structure and topology are only recently being revealed by helio-seismological data (Thompson et al. 2003; Mitra-Kraev & Thompson 2007; Gough & Hindman 2009). At present, such currents are not included in any OV model. Second, since  − R_ijS_ij > 0, the first term on the rhs of (1a) represents the production of K by shear. Differential rotation is included in this term whose evaluation requires the knowledge of R_ij, the Reynolds stresses whose general form was constructed in Papers I–III. No OV model has thus far included shear production. Third, the second and third terms on the rhs of (1a) represent the buoyancy flux that contains the possibility of DD processes (double diffusion: salt-fingers and semi-convection),whose effect we now study (Vauclair 2004, 2008; Charbonnel & Zahn 2007, and references cited there of previous work). Using the second of (1b), we have in the 1D case $\begin{array}{ccc}& & {\mathit{J}}_{\mathit{\rho }}\mathrm{=}{\mathit{\alpha }}_{\mathit{T}}{\mathit{J}}_{\mathit{h}}\mathrm{-}{\mathit{\alpha }}_{\mathit{\mu }}{\mathit{J}}_{\mathit{\mu }}\mathrm{=}{\mathit{H}}_{\mathrm{p}}^{-1}{\mathit{K}}_{\mathit{h}}\mathrm{(}\mathrm{\nabla }\mathrm{-}{\mathrm{\nabla }}_{\mathrm{ad}}\mathrm{)}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{r}}^{-1}\mathrm{)}\mathit{,}\\ & & \mathit{r}\mathrm{\equiv }\frac{{\mathit{K}}_{\mathit{h}}}{{\mathit{R}}_{\mathit{\mu }}{\mathit{K}}_{\mathit{\mu }}}\mathit{,}{\mathit{R}}_{\mathit{\mu }}\mathrm{=}\frac{{\mathrm{\nabla }}_{\mathit{\mu }}}{\mathrm{\nabla }\mathrm{-}{\mathrm{\nabla }}_{\mathrm{ad}}}\mathrm{·}\end{array}$(1c)In the case of salt-fingers, the T-field is stable and the μ-field is unstable and the instability of the μ-field generates mixing. Since the process is characterized by the relations N² > 0:R_μ < 1,   ∇ − ∇_ad < 0, we have $\frac{{\mathit{K}}_{\mathit{h}}}{{\mathit{K}}_{\mathit{\mu }}}\mathrm{(}\mathrm{\nabla }\mathrm{-}{\mathrm{\nabla }}_{\mathrm{ad}}\mathrm{)}\mathit{<}{\mathrm{\nabla }}_{\mathit{\mu }}\mathit{,} \mathit{r}\mathit{<}\mathrm{1}\mathit{,} {\mathit{K}}_{\mathit{h}}\mathit{<}{\mathit{K}}_{\mathit{\mu }}\mathit{,}$(1d)in which case ${\mathit{J}}_{\mathit{\rho }}\mathit{>}\mathrm{0.}$(1e)For semi-convection, the μ-field is stable and T-field is unstable and the instability of the T-field generates mixing. Since in this case the process is characterized by the relations N² > 0:   R_μ > 1,   ∇ − ∇_ad > 0, we have $\frac{{\mathit{K}}_{\mathit{h}}}{{\mathit{K}}_{\mathit{\mu }}}\mathrm{(}\mathrm{\nabla }\mathrm{-}{\mathrm{\nabla }}_{\mathrm{ad}}\mathrm{)}\mathit{>}{\mathrm{\nabla }}_{\mathit{\mu }}\mathit{,} \mathit{r}\mathit{>}\mathrm{1} \mathit{,} {\mathit{K}}_{\mathit{h}}\mathit{>}{\mathit{K}}_{\mathit{\mu }}\mathit{,}$(1f)in which case ${\mathit{J}}_{\mathit{\rho }}\mathit{>}\mathrm{0.}$(1g)As a result, both salt fingers and semi-convection generate a positive J_ρ which acts like a source of K that increases the extent of the OV. Within this background, one must consider that all present formulations of the OV extent within method c) have adopted the following equation for K: $\frac{\mathit{\partial K}}{\mathit{\partial t}}\mathrm{+}\underset{\mathrm{non}\mathrm{-}\mathrm{local}\hspace{0.17em}\mathrm{term}}{\underset{\mathrm{􏽼}0\mathrm{􏽻􏽺}0\mathrm{􏽽}}{\frac{\mathit{\partial }{\mathit{F}}_{\mathit{i}}^{\mathrm{ke}}}{\mathit{\partial }{\mathit{x}}_{\mathit{i}}}}}\mathrm{=}\mathit{g}{\mathit{\lambda }}_{\mathit{i}}{\mathit{\alpha }}_{\mathit{T}}{\mathit{J}}_{\mathit{i}}^{\mathit{h}}\mathrm{-}\mathit{\epsilon ,}$(1h)which is the limit of (1a) in the case of no shear (no differential rotation), no meridional currents and no DD. The effect of DD processes is particularly interesting, since even without numerical estimates, one can understand qualitatively what the effect is on the extent of the OV. Without DD, in Eq. (1h) J_h < 0 and thus the buoyancy flux term in (1h) is negative and acts like a sink of K. As such it tends to decrease the extent of the OV. Since the whole rhs of (1h) is now negative, the only way to avoid (1h) from representing a case of decaying turbulence is to have a negative gradient of the non-local term flux of K denoted by ${\mathit{F}}_{\mathit{i}}^{\mathrm{ke}}$, which represents the fourth problem represented by the need to model the third-order moment (TOM) ${\mathit{F}}_{\mathit{i}}^{\mathrm{ke}}$ in terms of some second-order moment. In their analysis, Xiong (1986) and Deng & Xiong (2008) adopted a down-gradient approximation, DGA, which, in the 1D case has the form ${\mathit{F}}^{\mathrm{ke}}\mathrm{=}\mathrm{-}{\mathit{K}}_{\mathrm{m}}\frac{\mathit{\partial K}}{\mathit{\partial z}}$(1i) where K_m is a momentum diffusivity. The problem of modeling TOMs was discussed in detail in Sect. 11 of Canuto (2009) where the relevant work by previous authors was discussed together with most recent advances. In spite of considerable progress in modeling TOMs documented in the work just cited, the present author is unable to suggest an expression for ${\mathit{F}}_{\mathit{i}}^{\mathrm{ke}}$ representative of the physical conditions characterizing the OV. This is because the assessments of the TOMs models were made using primarily planetary boundary layer processes that lack many key stellar features, such as meridional currents, differential rotation, radiative losses, and double-diffusion.

The reliability of the DGA closure (1i) in a stellar context is uncertain, to say the least, and yet the extent of the OV depends crucially on the structure of ${\mathit{F}}_{\mathit{i}}^{\mathrm{ke}}$. Several studies have shed light on the more general problem of the TOMs. The most extensive studies are those by Kupka & Muthsam (2007, 2008) and Kupka (2008), the general conclusion of which is that the DGAs rarely ever work, and if they do so, it is only near the boundary of convection zones. Due to these shortcomings, the key suggestion here is to abandon attempts at modeling ${\mathit{F}}_{\mathit{i}}^{\mathrm{ke}}$ and instead construct a differential equation for ${\mathit{F}}_{\mathit{i}}^{\mathrm{ke}}$ vs. radius. The point at which such flux has decreased to say 1% of its value at the bottom of the CZ can be considered a good measure of the OV extent. In the new procedure, we further include

• a)

  the effect of meridional currents,

• b)

  the flux of the total energy E;

• c)

  the flux of the Reynolds stresses; and

• d)

  the effect of unstable-stable stratification, differential rotation, shear, radiative losses (Peclet number), and double diffusion.

2. New approach

We consider the energy conservation law (I, Eq. (4a) with $\mathrm{D}\mathit{/}\mathrm{D}\mathit{t}\mathrm{=}{\mathit{\partial }}_{\mathit{t}}\mathrm{+}{\overline{\mathit{u}}}_{\mathit{i}}{\mathit{\partial }}_{\mathit{i}}\mathrm{)}$$\overline{\mathit{\rho }}\frac{\mathrm{D}\mathit{E}}{\mathrm{D}\mathit{t}}\mathrm{=}\mathrm{-}\frac{\mathit{\partial }{\mathrm{\Phi }}_{\mathit{i}}}{\mathit{\partial }{\mathit{x}}_{\mathit{i}}}\mathrm{+}\frac{\mathit{\partial }\overline{\mathit{p}}}{\mathit{\partial t}}$(2a)where $\mathit{E}\mathrm{=}{\mathit{c}}_{\mathrm{p}}\mathit{T}\mathrm{+}\mathit{K}\mathrm{+}{\mathit{K}}_{\mathrm{M}}\mathrm{+}\mathit{G,} {\mathit{K}}_{\mathrm{M}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}{\overline{{u}}}^{\mathrm{2}}\mathit{,} {\mathit{g}}_{\mathit{i}}{\overline{\mathit{u}}}_{\mathit{i}}\mathrm{=}\frac{\mathrm{D}\mathit{G}}{\mathrm{D}\mathit{t}}\mathit{,}$(2b)which represents the total energy of the system given by the sum of enthalpy c_pT, turbulent kinetic energy K, mean flow (meridional and azimuthal) kinetic energy K_M, and gravitational energy G. Furthermore, we have the total flux given by ${\mathrm{\Phi }}_{\mathit{i}}\mathrm{=}\overline{\mathit{\rho }}\mathrm{(}{\mathit{F}}_{\mathit{i}}^{\mathrm{rad}}\mathrm{+}{\mathit{F}}_{\mathit{i}}^{\mathrm{h}}\mathrm{+}{\mathit{F}}_{\mathit{i}}^{\mathrm{ke}}\mathrm{+}{\mathit{R}}_{\mathit{ij}}{\overline{\mathit{u}}}_{\mathit{j}}\mathrm{)}\mathit{,}$(2c)which represents the sum of the radiative, heat, and turbulent kinetic energy fluxes as well as the flux of the Reynolds stresses by the mean flow, where $\begin{array}{c}\\ \mathrm{radiative}\hspace{0.17em}\mathrm{flux}\mathrm{:} {\mathit{F}}_{\mathit{i}}^{\mathrm{r}}\mathrm{=}\mathrm{-}{\mathit{c}}_{\mathrm{p}}{\mathit{\chi }}_{\mathit{r}}\frac{\mathit{\partial T}}{\mathit{\partial }{\mathit{x}}_{\mathit{i}}}\\ \mathrm{heat}\hspace{0.17em}\hspace{0.17em}\mathrm{flux}\mathrm{:} {\mathit{F}}_{\mathit{i}}^{\mathrm{h}}\mathrm{=}{\mathit{c}}_{\mathrm{p}}{\mathit{J}}_{\mathit{i}}^{\mathrm{h}}\\ \mathrm{Reynolds}\hspace{0.17em}\hspace{0.17em}\mathrm{stresses}\mathrm{:} {\mathit{R}}_{\mathit{ij}}\mathrm{=}\overline{{\mathit{u}}_{\mathit{i}}^{\mathrm{\prime }}{\mathit{u}}_{\mathit{j}}^{\mathrm{\prime }}}\mathit{,}\end{array}$(2d)where χ_r is the thermometric conductivity with units of cm² s^-1. Next, we consider Eqs. (I, 5) for the mean concentration: $\overline{\mathit{\rho }}\frac{\mathrm{D}\overline{\mathit{C}}}{\mathrm{D}\mathit{t}}\mathrm{=}\mathrm{-}\frac{\mathit{\partial }}{\mathit{\partial }{\mathit{x}}_{\mathit{i}}}\mathrm{\left(}\overline{\mathit{\rho }}{\mathit{J}}_{\mathit{i}}^{\mathit{c}}\mathrm{\right)}\mathrm{-}\frac{\mathit{\partial }}{\mathit{\partial }{\mathit{x}}_{\mathit{i}}}\mathrm{\left(}\overline{\mathit{\rho }}{\mathit{F}}_{\mathit{i}}^{\mathit{c}}\mathrm{\right)}\mathit{,} {\mathit{F}}_{\mathit{i}}^{\mathit{c}}\mathrm{=}\mathrm{-}{\mathit{\chi }}_{\mathit{c}}\frac{\mathit{\partial }\overline{\mathit{C}}}{\mathit{\partial }{\mathit{x}}_{\mathit{i}}}$(3a)where ${\mathit{J}}_{\mathit{i}}^{\mathit{c}}\mathrm{=}\overline{{\mathit{u}}_{\mathit{i}}^{\mathrm{\prime }}{\mathit{c}}^{\mathrm{\prime }}}$. The variables ${\mathit{F}}_{\mathit{i}}^{\mathit{c}}$ and χ_c are the physical equivalents of ${\mathit{F}}_{\mathit{i}}^{\mathit{r}}$ and χ_r. We must note that χ_c is often called D in the literature, e.g., Michaud (1970, Eq. (1)) and Vauclair & Vauclair (1983, Eq. (1)). Multiplying Eqs. (2a) and (3a) by α_T,α_c, respectively, and subtracting the two, one obtains an equation of the same structure of (2a), namely, $\overline{\mathit{\rho }}\frac{\mathrm{D}{\mathit{E}}_{\mathrm{\ast }}}{\mathrm{D}\mathit{t}}\mathrm{=}\mathrm{-}\frac{\mathit{\partial }{\mathrm{\Phi }}_{\mathit{i}}^{\mathrm{\ast }}}{\mathit{\partial }{\mathit{x}}_{\mathit{i}}}\mathrm{+}\frac{\mathit{\partial }\overline{{\mathit{p}}_{\mathrm{\ast }}}}{\mathit{\partial t}}$(3b)where ${\mathit{E}}_{\mathrm{\ast }}\mathrm{=}{\mathit{c}}_{\mathrm{p}}\mathrm{\left[}{\mathit{\alpha }}_{\mathit{T}}\mathit{T}\mathrm{\right(}\mathrm{1}\mathrm{+}\mathit{\xi }\mathrm{)}\mathrm{-}{\mathit{\alpha }}_{\mathit{c}}\mathit{C}\mathrm{]}\mathit{,}\hspace{0.17em}\mathit{\xi }\mathrm{=}\frac{\mathit{K}\mathrm{+}{\mathit{K}}_{\mathrm{M}}\mathrm{+}\mathit{G}}{{\mathit{c}}_{\mathrm{p}}\mathit{T}}\mathit{,}\hspace{0.17em}\overline{{\mathit{p}}_{\mathrm{\ast }}}\mathrm{=}{\mathit{\alpha }}_{\mathit{T}}\overline{\mathit{p}}$(3c)${\overline{\mathit{\rho }}}^{-1}{\mathrm{\Phi }}_{\mathit{i}}^{\mathrm{\ast }}\mathrm{=}{\mathit{c}}_{\mathrm{p}}{\mathit{J}}_{\mathit{i}}^{\mathit{\rho }}\mathrm{+}{\mathit{\alpha }}_{\mathit{T}}{\mathit{F}}_{\mathit{i}}^{\mathrm{ke}}\mathrm{-}{\mathit{c}}_{\mathrm{p}}{\mathit{J}}_{\mathit{i}}^{\mathrm{r}\mathrm{-}\mathrm{m}}\mathrm{+}{\mathit{\alpha }}_{\mathit{T}}{\overline{\mathit{u}}}_{\mathit{j}}{\mathit{R}}_{\mathit{ij}}\mathit{,}$(3d)where for brevity we call ${\mathit{F}}_{\mathit{i}}^{\mathrm{r}\mathrm{-}\mathrm{m}}$ the radiative–molecular terms: ${\mathit{J}}_{\mathit{i}}^{\mathrm{r}\mathrm{-}\mathrm{m}}\mathrm{\equiv }{\mathit{\alpha }}_{\mathit{T}}{\mathit{\chi }}_{\mathit{r}}\frac{\mathit{\partial T}}{\mathit{\partial }{\mathit{x}}_{\mathit{i}}}\mathrm{-}{\mathit{\alpha }}_{\mathit{c}}{\mathit{\chi }}_{\mathit{c}}\frac{\mathit{\partial }\overline{\mathit{C}}}{\mathit{\partial }{\mathit{x}}_{\mathit{i}}}\mathrm{·}$(3e)Clearly, the variable ξ represents the ratio of the sum of the turbulent energy K, the mean flow energy K_M, and the gravitational energy to the thermal energy. The first four terms in (3d) represent the turbulent fluxes of buoyancy, eddy kinetic energy, radiative-molecular fluxes. The last term represents the flux of the Reynolds stresses R_ij by the large-scale fields ${\overline{\mathit{u}}}_{\mathit{j}}$. In the stationary limit, Eqs. (3b)–(3d) reduce to the following generalized flux conservation law: $\begin{array}{}\underset{\mathrm{􏽼}0\mathrm{􏽻􏽺}0\mathrm{􏽽}}{{\mathit{c}}_{\mathrm{p}}{\mathit{J}}_{\mathit{i}}^{\mathit{\rho }}\mathrm{+}{\mathit{\alpha }}_{\mathit{T}}{\mathit{F}}_{\mathit{i}}^{\mathrm{ke}}}\\ \mathrm{turbulentfluxes}\end{array}\mathrm{-}\underset{\mathrm{radiative}\mathrm{-}\mathrm{molecular}\hspace{0.17em}\mathrm{fluxes}\hspace{0.17em}}{\underset{\mathrm{􏽼}0\mathrm{􏽻􏽺}0\mathrm{􏽽}}{{\mathit{c}}_{\mathrm{p}}{\mathit{J}}_{\mathit{i}}^{\mathrm{r}\mathrm{-}\mathrm{m}}}}\mathrm{+}\hspace{0.17em}\begin{array}{}\underset{\mathrm{􏽼}0\mathrm{􏽻􏽺}0\mathrm{􏽽}}{{\overline{\mathit{u}}}_{\mathit{j}}\mathrm{(}{\mathit{E}}_{\mathrm{\ast }}{\mathit{\delta }}_{\mathit{ij}}\mathrm{+}{\mathit{\alpha }}_{\mathit{T}}{\mathit{R}}_{\mathit{ij}}\mathrm{)}}\\ \mathrm{large}\hspace{0.17em}\mathrm{scale}\hspace{0.17em}\mathrm{fluxes}\end{array}\mathrm{=}\hspace{0.17em}\mathrm{const}\mathit{.}$(4a)We note that each term in (4a) has the units of [α_Tu³]. The “traditional” flux conservation law used in stellar models contains only the first and the radiative parts, in which case Eq. (4a) reduces to ${\mathit{J}}_{\mathit{i}}^{\mathit{h}}\mathrm{-}{\mathit{\chi }}_{\mathit{r}}\frac{\mathit{\partial T}}{\mathit{\partial }{\mathit{x}}_{\mathit{i}}}\mathrm{=}\hspace{0.17em}\mathrm{const}\mathit{.}\mathit{,}$(4b)which expresses the constancy of the convective plus radiative fluxes. Eliminating the buoyancy flux ${\mathit{J}}_{\mathit{i}}^{\mathit{\rho }}$ between Eqs. (1a) and (4a), we obtain a differential equation for ${\mathit{F}}_{\mathit{i}}^{\mathrm{ke}}$: $\frac{\mathit{\partial }{\mathit{F}}_{\mathit{i}}^{\mathrm{ke}}}{\mathit{\partial }{\mathit{x}}_{\mathit{i}}}\mathrm{+}\left(\frac{\mathit{g}{\mathit{\alpha }}_{\mathit{T}}}{{\mathit{c}}_{\mathrm{p}}}\right){\mathit{\lambda }}_{\mathit{i}}{\mathit{F}}_{\mathit{i}}^{\mathrm{ke}}\mathrm{=}{\mathrm{\Phi }}_{\mathrm{old}}\mathrm{+}{\mathrm{\Phi }}_{\mathrm{new}}\mathit{,}$(5a)where we have used the notation $\begin{array}{ccc}& & {\mathrm{\Phi }}_{\mathrm{old}}\mathrm{=}\mathit{C}\mathrm{-}\mathit{g}{\mathit{\lambda }}_{\mathit{i}}{\mathit{J}}_{\mathit{i}}^{\mathrm{r}\mathrm{-}\mathrm{m}}\mathrm{-}\mathit{\epsilon ,}\\ & & {\mathrm{\Phi }}_{\mathrm{new}}\mathrm{=}\mathrm{-}{\mathit{R}}_{\mathit{ij}}{\mathit{S}}_{\mathit{ij}}\mathrm{-}\frac{\mathit{g}}{{\mathit{c}}_{\mathrm{p}}}{\mathit{\lambda }}_{\mathit{i}}{\overline{\mathit{u}}}_{\mathit{j}}\mathrm{(}{\mathit{E}}_{\mathrm{\ast }}{\mathit{\delta }}_{\mathit{ij}}\mathrm{+}{\mathit{\alpha }}_{\mathit{T}}{\mathit{R}}_{\mathit{ij}}\mathrm{)}\mathrm{-}{\overline{\mathit{u}}}_{\mathit{j}}\frac{\mathit{\partial K}}{\mathit{\partial }{\mathit{x}}_{\mathit{j}}}\mathrm{·}\end{array}$(5b)Each term of Eq. (5a) has the units of a dissipation [L² / t³]. If one keeps only Φ_old, Eq. (5a) coincides with the equation originally suggested by Roxburgh (1978). Since the constant C is related to the total stellar luminosity, the final criterion is expressed in terms of the integral of the difference L(N) − L(r) between the nuclear and radiative luminosities. Since such a difference is positive in the convective zone CZ and negative in the OV, the integral changes sign at some distance r_∗ that was identified with the OV extent. The presence of the new term Φ_new in (5a) introduces the following new features $\begin{array}{ccc}& & \mathrm{differentialrotation,largescalecurrents,concentrationfield,}\\ & & \mathrm{stable}\mathrm{-}\mathrm{unstablestratificationandradiativelosses}\mathrm{.}\end{array}$(5c)The solution of the new equation for the flux of K, Eq. (5a) can only be carried out in conjunction with a model for the Reynolds stresses, which we have discussed in Papers II and III. One also needs a stellar code to provide the radiative flux ${\mathit{F}}_{\mathit{i}}^{\mathrm{rad}}$. Finally, one must model the dissipation ε, a topic that we discussed in Papers I and II and which we slightly amplify here by listing some new options.

• a)

  The most obvious choice consists of adopting a Kolmogorovtype law, $\mathit{\epsilon }\mathrm{=}{\mathit{K}}^{\mathrm{3}\mathit{/}\mathrm{2}}{\mathit{\ell }}^{-1}\mathit{,}$(5d)where ℓ is a dissipation scale. In a stably stratified regime, such as the one encountered in a stellar radiative zone, the validity of the Kolmogorov spectrum is questionable and though one can formally write (5d), it is of little help since ℓ turns out to be a rather complex function of the stability functions such as the flux Richardson number (Cheng & Canuto 1994). This option is not advisable because it introduces a variable that is difficult to control.

• b)

  Following Paper II, Eq. (10a), one can identify ε with the estimate provided by the internal gravity waves model (IGW) $\mathit{\epsilon }\mathrm{=}\mathit{\epsilon }\mathrm{\left(}\mathrm{IGW}\mathrm{\right)}\mathit{,}$(5e)which in the solar case applies to the radiative zone below the convective zone, but which may be nontrivial to extend to other stars.

• c)

  The most reliable procedure is probably the one discussed in Sect. 6 of Paper I, which calls for the solution of the following differential equation for ε (c_1,2 = 2.88, 3.84) $\mathit{\tau }\frac{\mathrm{D}\mathit{\epsilon }}{\mathrm{D}\mathit{t}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}\frac{\mathit{\partial }{\mathit{F}}_{\mathit{i}}^{\mathrm{ke}}}{\mathit{\partial }{\mathit{x}}_{\mathit{i}}}\mathrm{+}{\mathit{c}}_{\mathrm{1}}{\mathit{P}}_{\mathit{s}}\mathrm{+}{\mathit{c}}_{\mathrm{3}}{\mathit{P}}_{\mathit{b}}\hspace{0.17em}\mathrm{-}\hspace{0.17em}{\mathit{c}}_{\mathrm{2}}\mathit{\epsilon ,}$(5f)

where we have used the form of the flux of ε suggested by Kupka & Mutsham (2007): ${\mathit{F}}_{\mathit{i}}^{\mathit{\epsilon }}\mathrm{=} \mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}{\mathit{\tau }}^{-1}{\mathit{F}}_{\mathit{i}}^{\mathrm{ke}}\mathit{.}$(5g)The presence of ${\mathit{F}}_{\mathit{i}}^{\mathrm{ke}}$ in (5f) means that the latter must be solved with Eqs. (5a), (5b). Even this procedure still requires the knowledge of one variable, the dissipation time scale τ. Since we no longer have an independent differential equation for the latter (since we have used both the K and ε equations), we suggest a parametric form consisting of assuming that in stably stratified flows: $\mathit{\tau N}\mathrm{=}\mathit{O}\mathrm{\left(}\mathrm{1}\mathrm{\right)}\mathit{.}$(5h)Admittedly, this introduces a proportionality constant of order unity in (5h), but we consider this a far more physically justified and general procedure than the previous choices.

Once Eqs. (5a), (5b) are solved with (5f), the OV extent in the radiative zone can be obtained as the distance from the bottom of the CZ where the flux of the kinetic energy ${\mathit{F}}_{\mathit{i}}^{\mathrm{ke}}$ has decreased, to say, 1% of its value at the bottom of the CZ.

Since we have eliminated the buoyancy flux that no longer appears explicitly in (5a), one may wonder about the effect of DD on the extent of the OV.

3. Physical interpretation of the new term Φ_new

The novel feature in (5a) is represented by Φ_new, which brings about features that were not accounted for before and which we now analyze. The first term in Φ_new represents the production of kinetic energy by shear: ${\mathit{P}}_{\mathit{s}}\mathrm{=}\mathrm{-}{\mathit{R}}_{\mathit{ij}}{\mathit{S}}_{\mathit{ij}}\mathit{>}\mathrm{0.}$(6a)The second term represents the flux by the mean currents of the total energy of the system $\mathrm{-}\mathit{g}{\mathit{c}}_{\mathrm{p}}^{-1}{\mathit{\lambda }}_{\mathit{i}}{\overline{\mathit{u}}}_{\mathit{i}}{\mathit{E}}_{\mathrm{\ast }}\mathit{,}$(6b)while the third term represents the flux of the Reynolds stresses by the mean currents $\mathrm{-}\mathit{g}{\mathit{\alpha }}_{\mathit{T}}{\mathit{c}}_{\mathrm{p}}^{-1}{\mathit{\lambda }}_{\mathit{i}}{\overline{\mathit{u}}}_{\mathit{j}}{\mathit{R}}_{\mathit{ij}}\mathit{,}$(6c)which we now analyze separately. Using the results presented in Paper IV, we have ${\mathit{P}}_{\mathit{s}}\mathrm{=}\mathrm{-}{\mathit{R}}_{\mathit{ij}}{\mathit{S}}_{\mathit{ij}}\mathrm{=}\mathrm{-}\mathrm{2}{\mathit{R}}_{\mathit{r\phi }}{\mathit{S}}_{\mathit{r\phi }}\mathrm{=}\mathit{A\epsilon ,}$(7a)where (here r is the radial distance): $\mathit{A}\mathrm{\equiv }\frac{{\mathrm{\Gamma }}^{\mathrm{2}}}{\mathrm{4}}\mathrm{(}\mathit{\tau }\mathrm{\Omega }{\mathrm{)}}^{\mathrm{2}}{\mathit{S}}_{\mathit{m}}{\left(\frac{\mathit{r}}{\mathrm{\Omega }}\frac{\mathit{\partial }\mathrm{\Omega }}{\mathit{\partial r}}\right)}^{\mathrm{2}}\mathit{,} \mathrm{\Gamma }\mathrm{=}\sin \mathit{\theta }$(7b)where the structure function S_m was plotted in Fig. 2 of Paper IV. Next, we consider the second term in Φ_new, which is given by ${\mathrm{\Phi }}_{\mathrm{new}}\mathrm{\left(}\mathrm{secondterm}\mathrm{\right)}\mathrm{=}\mathrm{-}\mathit{g}{\mathit{c}}_{\mathrm{p}}^{-1}{\overline{\mathit{u}}}_{\mathit{r}}{\mathit{E}}_{\mathrm{\ast }}\mathit{.}$(8)Since the available data on the meridional currents by Mitra-Kraev & Thompson (2007) and Gough & Hindman (2009) do not extend below the convective zone, we are unable at present to determine the size of this term. The third term is given by $\begin{array}{ccc}{\mathrm{\Phi }}_{\mathrm{new}}\mathrm{\left(}\mathrm{thirdterm}\mathrm{\right)}& \mathrm{=}& \mathrm{-}\mathit{g}{\mathit{c}}_{\mathrm{p}}^{-1}{\mathit{\alpha }}_{\mathit{T}}{\mathit{\lambda }}_{\mathit{r}}{\mathit{R}}_{\mathit{rj}}{\overline{\mathit{u}}}_{\mathit{j}}\\ & \mathrm{=}& \mathrm{-}\mathit{g}{\mathit{c}}_{\mathrm{p}}^{-1}{\mathit{\alpha }}_{\mathit{T}}\mathrm{(}{\mathit{R}}_{\mathit{rr}}{\overline{\mathit{u}}}_{\mathit{r}}\mathrm{+}{\mathit{R}}_{\mathit{r\phi }}{\overline{\mathit{u}}}_{\mathit{\phi }}\mathrm{+}{\mathit{R}}_{\mathit{r\theta }}{\overline{\mathit{u}}}_{\mathit{\theta }}\mathrm{)}\mathit{,}\end{array}$(9a)where the dynamic equations for the meridional velocities ${\overline{\mathit{u}}}_{\mathit{r}}\mathit{,}{\overline{\mathit{u}}}_{\mathit{\theta }}$ are given by Eqs. (8) of Paper IV. If we assume that the second term in the last expression on the rhs of (9) is the largest, we then have, with ${\overline{\mathit{u}}}_{\mathit{\phi }}\mathrm{=}\mathit{r}\mathrm{\Gamma \Omega }$: ${\mathrm{\Phi }}_{\mathrm{new}}\mathrm{\left(}\mathrm{thirdterm}\mathrm{\right)}\mathrm{\approx }\mathrm{-}\mathit{g}{\mathit{c}}_{\mathrm{p}}^{-1}{\mathit{\alpha }}_{\mathit{T}}{\mathit{R}}_{\mathit{r\phi }}{\overline{\mathit{u}}}_{\mathit{\phi }}\mathrm{\approx }\mathrm{-}\mathit{B\epsilon ,}$(9b)where $\mathit{B}\mathrm{=}\frac{{\mathrm{\Gamma }}^{\mathrm{2}}}{\mathrm{2}}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{\gamma }}^{-1}\mathrm{)}\mathrm{(}\mathit{\tau }\mathrm{\Omega }{\mathrm{)}}^{\mathrm{2}}{\mathit{S}}_{\mathit{m}}\left(\mathrm{-}\frac{\mathit{r}}{\mathrm{\Omega }}\frac{\mathit{\partial }\mathrm{\Omega }}{\mathit{\partial r}}\right)\frac{\mathit{r}}{{\mathit{H}}_{\mathrm{p}}}\mathit{>}\mathrm{0.}$(9c)

4. Conclusions

In this paper we have put forward the suggestion of abandoning the canonical attempt to solve Eq. (1a) which requires a model for the flux of kinetic energy ${\mathit{F}}_{\mathit{i}}^{\mathrm{ke}}$; see Eq. (1i) and discussion thereafter, for which the present models are neither unique nor fully trustworthy. Rather, we suggest viewing the latter function as an unknown to be determined by solving the differential equation, Eq. (5a). The previous form of (5a) suggested by Roxburgh (1978) did not include the new term Φ_new, which brings the following new features:

• a)

  the effect of meridional currents,

• b)

  the flux of the total energy E,

• c)

  the flux of the Reynolds stresses and

• d)

  a full model of the Reynolds stresses that accounts for stable stratification, unstable stratification, differential rotation, shear, radiative losses (Peclet number) and double diffusion.

References