BEYOND MIXING-LENGTH THEORY: A STEP TOWARD 321D

Erika Böhm-Vitense developed the version of MLT used in stellar evolution in the 1950s (Vitense 1953; Böhm-Vitense 1958), prior to the publication in the west of Andrey Kolmogorov's theory of the turbulent cascade (Kolmogorov 1962). MLT might have been different had she been aware of the original work (Kolmogorov 1941). Edward Lorenz showed that a simple convective roll had chaotic behavior (a strange attractor, Lorenz 1963). Ludwig Prandtl developed the theory of boundary layers (Prandtl & Tietjens 1934), as well as the original version of MLT (Prandtl 1925). All these ideas will be relevant to our discussion, which is based, as far as possible, upon experimentally verified turbulence theory and 3D simulations, and free of astronomical calibration.

The 3D turbulent energy cascade is illustrated in Figure 1. The turbulent motion is driven at the largest scale (the "integral" scale), which contains most of the kinetic energy. These motions are unstable and break up into smaller-scale flow patterns dominated by inertial forces (the "inertial subrange"). This continues to scales small enough for microscopic effects (viscosity) to finally provide damping of the flow at the Kolmogorov scale. Both the inertial subrange and the dissipation range are insensitive to the details of the boundary conditions at the integral scale, and are "universal" in this sense. We use the term "universality" to mean the property of insensitivity to boundary conditions at the integral scale. Kolmogorov (1941) found the striking result that the rate of dissipation is insensitive to the value of the viscosity, but is determined by the rate that the largest-scale flows feed the cascade. This behavior of the nonlinear flow "hides" the microscopic value of the viscosity. We use Kolmogorov theory to describe the flow in the range where universality holds.

Zoom In Zoom Out Reset image size

Direct Numerical Simulations (DNS) resolve the small scales at which dissipation happens, and can extend up to the inertial range, but not to stellar scales. Implicit Large Eddy Simulations (ILES) can extend from stellar (integral) scales down to the inertial range, but not to the dissipation range. Figure 1 illustrates both.

Landau objected to the notion of complete universality on the grounds that the largest scales were subject to boundary conditions that would be specific to the case in question (Landau & Lifshitz 1959; Frisch 1995). We will incorporate this idea by splitting the turbulent flow into two parts: the integral-scale motion and the turbulent cascade. As an aid to understanding the integrated properties of the integral-scale motion, we are guided by the simplest model of a convective roll, due to Lorenz (1963). This model contains the famous Lorenz strange attractor, and exhibits chaotic behavior. It also agrees surprisingly well with 3D simulations of turbulent convection associated with oxygen burning prior to core collapse (Meakin & Arnett 2007b; Arnett & Meakin 2011b). This approximation does lack multi-mode behavior, as compared with the simulations, which are dominated by five low order modes (see Figure 1 in Arnett & Meakin 2011b); this may affect the accuracy of the representation of intermittency at large scales and of coherent structures.

Our challenge is to simplify this very complex problem, with time dependence and an astronomically large number of degrees of freedom, down to a feasible level for use in a stellar evolutionary code, without losing important features. Our approximation, 321D, is an attempt to increase physical realism at feasible cost in computational complexity. It is desirable to avoid astronomical calibration as far as possible, and base changes upon behavior quantified in laboratory and numerical experiments. In particular, we do not validate our approximation by how well it reproduces standard MLT results. By basing approximations on 3D ILES simulations that (1) exhibit turbulence, (2) have non-uniform composition, and (3) resolve dynamic boundary behavior, it is possible to remove some of the vagueness inherent in many theoretical treatments of convection.

We will compare the global properties of turbulent convection from numerical and analytical viewpoints in Section 2, examine the structure and nature of boundaries of convection zones in Section 3, and summarize our conclusions in Section 4. In the Appendix we provide a derivation from 3D fluid flow equations for some useful expressions.

Using a classical radiative viscosity (Mihalas & Mihalas 1984), the Reynolds number is $\mathrm{Re}\sim {10}^{8}$ at the base of the solar convection zone.⁹ Numerical simulations and laboratory experiments become turbulent for roughly $\mathrm{Re}\geqslant {10}^{3}$, so fluid flows in stars are strongly turbulent if, as we assume for the moment, rotational and magnetic field effects may be neglected.

For homogeneous, isotropic, and steady-state turbulence, the Kolmogorov relation (Frisch 1995) between the dissipation rate of TKE per unit mass ${\epsilon }_{t}$, velocity v, and length scale ℓ is

$\begin{eqnarray}&&{\epsilon }_{t}=\displaystyle \frac{4}{5}{v}^{3}/{\ell }.\end{eqnarray} \tag{ 1 }$

Arnett et al. (2009) found that ${\epsilon }_{t}=0.85\ {v}_{\mathrm{rms}}^{3}/{{\ell }}_{\mathrm{cz}}$, where ℓ_cz is the depth of the convective zone, and v_rms is the average convective velocity; see their Equation (6) and nearby discussion, and references to other studies which report such coefficients. For homogeneous, isotropic turbulence, Kolmogorov (1941) predicted a coefficient $4/5$ for a region well away from boundaries. This factor of 0.8 might change for a strongly stratified region, which would have flow better described by plumes than convective rolls.

Equation (1) is a global constraint, averaged over fluctuations, and applies to each length scale λ in the turbulent cascade, so

$\begin{eqnarray}&&{\epsilon }_{t}\sim {({\rm{\Delta }}{v}_{\lambda })}^{3}/\lambda ,\end{eqnarray} \tag{ 2 }$

for all scales λ, or,

$\begin{eqnarray}&&{\rm{\Delta }}{v}_{\lambda }\sim {({\epsilon }_{t}\lambda )}^{\displaystyle \frac{1}{3}}.\end{eqnarray} \tag{ 3 }$

so that the velocity variation across a scale λ is ${\rm{\Delta }}{v}_{\lambda }$, which increases as ${\lambda }^{\displaystyle \frac{1}{3}}$. The largest scales have the largest velocities, and are dominated by advective transport (macroscopic mixing).

The velocity gradient across the scale λ is

$\begin{eqnarray}&&{\rm{\Delta }}{v}_{\lambda }/\lambda \sim {\epsilon }_{t}^{\displaystyle \frac{1}{3}}/{\lambda }^{\displaystyle \frac{2}{3}},\end{eqnarray} \tag{ 4 }$

and increases with decreasing λ. The smallest scales have the largest velocity gradients, and are eventually dominated by microscopic mixing (ionic diffusion, radiative diffusion, and viscosity). A description of the cascade needs both large and small scales; Equation (3) implies that the largest (integral) scales have most of the kinetic energy and momentum, while Equation (4) implies that the smallest scales have the fastest relaxation times, which is consistent with simulations (e.g., Arnett et al. 2009).

Landau & Lifshitz (1959), Section 32, estimated the number of degrees of freedom in a region of turbulent flow to be $N\sim {(\mathrm{Re})}^{9/4}$. Laminar flows with free boundaries become unstable at roughly $\mathrm{Re}\sim {10}^{3}$. A DNS would require well over 10⁸ zones to resolve the cascade for this marginally unstable case. Using $\mathrm{Re}\sim {10}^{8}$ (see Section 2.1), implies a need for more than 10¹⁸ zones for the Sun, far beyond current computer capacity.

There may be a smarter way. Kolmogorov's great insight is that turbulence hides the details of the viscous dissipation by the nonlinear interactions of the cascade, so that the dissipation rate is determined by macroscopic parameters. Simulations show a multimode behavior (Meakin & Arnett 2007b), but only $N\sim 5$ dominant modes¹⁰ for $\sim {10}^{8}$ zones. This dramatic reduction in complexity suggests the use of ILES, (see Figure 1 and Boris 2007) which approximate small scale behavior by a Kolmogorov cascade. Our approach is to assume that this simplification holds for very large Reynolds numbers, and to examine the consequences. Simulations which are presently feasible have effective Reynolds numbers limited by numerical resolution, but are sufficiently high to give truly turbulent solutions. State of the art simulations, with both improved algorithms and more powerful computers, support this approach (Porter & Woodward 2000; Herwig et al. 2014; S. Campbell et al. 2015, in preparation).

As an aid to the reader, Table 1 gives the correspondence of selected variables in three different theoretical approaches to turbulent convection: MLT, the Lorenz model, and the RANS formulation. MLT is 1D (radial), the Lorenz model is 2D (radial and transverse), while the RANS analysis is 3D projected to 1D. MLT is static, the Lorenz model and the RANS equations are time dependent. MLT is local (no spatial derivatives of velocity) while the Lorenz model is mildly nonlocal (it uses global derivatives over the roll), and the RANS equations are non-local. Comparison of MLT and Lorenz gives a sense of transverse versus radial properties.

Table 1.  Correspondence of Some Variables in MLT, Lorenz, and RANS

Quantitya	MLTb	Lorenzc	RANSd	Comment
Dissipation length	$({\alpha }^{2}/8){H}_{P}$	⋯	${{\ell }}_{{\rm{d}}}\approx 0.8{{\ell }}_{\mathrm{CZ}}$	${{\ell }}_{\mathrm{CZ}}$ is convection zone depth
	⋯	⋯	⋯	Meakin & Arnett (2007b)
Horizontala gradient	${\rm{\Delta }}{\rm{\nabla }}={\rm{\nabla }}-{{\rm{\nabla }}}_{e}$	$\left(\displaystyle \frac{2{H}_{P}}{{\ell }}\right){T}_{3}/{T}_{0}$	⋯	Smith & Arnett (2014)
Radial gradient	${{\rm{\nabla }}}_{e}-{{\rm{\nabla }}}_{a}$	$\left(\displaystyle \frac{2{H}_{P}}{{\ell }}\right){T}_{2}/{T}_{0}$	⋯	⋯
Imposed gradient	${{\rm{\nabla }}}_{r}-{{\rm{\nabla }}}_{a}$	$\left(\displaystyle \frac{2{H}_{P}}{{\ell }}\right){T}_{1}/{T}_{0}$	⋯	⋯
Convective velocity	Equation (5)	u	u'	Algebraic (MLT) versus ode
	⋯	⋯	⋯	Local (MLT) versus nonlocal
Turbulent heating	None	Ignored or	$\langle {(u^{\prime} )}^{3}\rangle /{{\ell }}_{{\rm{d}}}$	Arnett et al. (2009)
	⋯	${u}^{2}| u| /{{\ell }}_{{\rm{d}}}$	⋯	⋯
Kinetic energy flux	Assumed	Assumed	$\langle {\rho }^{\prime }u^{\prime} {\boldsymbol{u}}\cdot {\boldsymbol{u}}/2\rangle$	Meakin & Arnett (2007b, 2010)
	cancellation	cancellation	No cancellation,	⋯
	by symmetry	by symmetry	asymmetry	⋯
Buoyancy flux	$u{\beta }_{T}g{\rm{\Delta }}{\rm{\nabla }}$	$\displaystyle \frac{1}{2}{\beta }_{T}{{guT}}_{3}/{T}_{0}$	$-g\langle {\rho }^{\prime }u^{\prime} \rangle /{\rho }_{0}$	MLT ignores composition gradients
Enthalpy flux	$\rho {{uC}}_{P}T{\rm{\Delta }}{\rm{\nabla }}$	$\displaystyle \frac{1}{2}\rho {C}_{P}{{uT}}_{3}$	$\rho {C}_{P}\langle u^{\prime} T^{\prime} \rangle$	Viallet et al. (2013)
Acoustic energy flux	None	None	$\langle P^{\prime} u^{\prime} \rangle$	Small for low-mach flow
Compositione flux	Undefined	None	$\rho \langle Y^{\prime} u^{\prime} \rangle$	Arnett (1996)
Ye flux	Undefined	None	$\rho \langle {Y}_{e}^{\prime }u^{\prime} \rangle$	⋯

Notes.

^aThe MLT variables are all defined in the radial direction. RANS projects a 3D average onto the radial direction. The Lorenz model has both radial and horizontal gradients (Arnett & Meakin 2011b; Smith & Arnett 2014). ^bSmith & Arnett (2014). ^cArnett et al. (2009), Arnett & Meakin (2011b); ℓ is the roll diameter. ^dMeakin & Arnett (2007b), Viallet et al. (2013). ^eArnett (1996), $Y={Y}_{e}+{{\rm{\Sigma }}}_{i}{Y}_{i}={Y}_{e}+1/\bar{A}$ and ${Y}_{e}={{\rm{\Sigma }}}_{i}{Z}_{i}{Y}_{i}$.

Download table as:  ASCIITypeset image

In MLT the buoyant acceleration is approximately integrated over a mixing length ${{\ell }}_{\mathrm{MLT}}$ to obtain an average velocity u (e.g., Vitense 1953; Böhm-Vitense 1958; Kippenhahn & Weigert 1990),

$\begin{eqnarray}&&{u}^{2}=g{\beta }_{T}{\rm{\Delta }}{\rm{\nabla }}\left(\displaystyle \frac{{{\ell }}_{\mathrm{MLT}}^{2}}{8{H}_{P}}\right).\end{eqnarray} \tag{ 5 }$

The superadiabatic excess ${\rm{\Delta }}{\rm{\nabla }}$ is defined in Table 1 and Section 2.4. Here g is the gravitational acceleration, ${\beta }_{T}=-{(\partial \mathrm{ln}\rho /\partial \mathrm{ln}T)}_{P}$ is a thermodynamic variable (for uniform composition; see Section 2.4 for the nonuniform case), H_P is the local pressure scale height, and ${{\ell }}_{\mathrm{MLT}}$ is an adjustable length scale (the mixing length).

Equation (5) requires that ${\rm{\Delta }}{\rm{\nabla }}\geqslant 0$ for the velocity u to be a real number. The velocity depends only on the local value of the superadiabatic gradient ${\rm{\Delta }}{\rm{\nabla }}$. There are obvious problems with regions in which such integration extends past a boundary.

There have been a number of attempts to generalize MLT (e.g., Unno (1961), Gough (1967, 1977), Arnett (1969), Stellingwerf (1976), Kuhfuss (1986) ,Xiong (1986), Deng et al. (1996, 2006), Xiong et al. (1997), Hansen & Kawaler (1994), etc.). Working backward, Equation (5) may be expressed as a co-moving acceleration equation for a vector field ${\boldsymbol{u}}$:

$\begin{eqnarray}&&d{\boldsymbol{u}}/{dt}={\mathcal{B}}-{\mathcal{D}},\end{eqnarray} \tag{ 6 }$

where ${\mathcal{B}}$ is a generalized driving term and ${\mathcal{D}}$ a corresponding drag term (Prandtl & Tietjens 1934, Ch. V). A hydrostatic background will be assumed; see Appendix Section A. Similar equations result from (1) study of the nonlinear development of the Rayleigh–Taylor instability (RTI), and from (2) applications of RANS analysis to 3D simulations of turbulent convection.

If the driving is due to buoyancy alone, (see Section 2.4 for nonuniform composition), $-g(\delta \rho /\rho )\approx {\boldsymbol{g}}{\beta }_{T}{\rm{\Delta }}{\rm{\nabla }}$, then ${\mathcal{B}}\approx {\boldsymbol{g}}{\beta }_{T}{\rm{\Delta }}{\rm{\nabla }}$. If the drag is represented by ${\mathcal{D}}\approx {\boldsymbol{u}}/\tau$, where $\tau ={{\ell }}_{{\rm{d}}}/| u|$, then we have

$\begin{eqnarray}&&d{\boldsymbol{u}}/{dt}=\partial {\boldsymbol{u}}/\partial t+({\boldsymbol{u}}\cdot {\rm{\nabla }}){\boldsymbol{u}}={\boldsymbol{g}}{\beta }_{T}{\rm{\Delta }}{\rm{\nabla }}-{\boldsymbol{u}}/\tau .\end{eqnarray} \tag{ 7 }$

This is basically a statement of Newtonian mechanics, with driving by buoyancy and damping by drag. Gough (1977) gives a historical context going back to Prandtl (1925) and to Biermann (1932). The early attempts, and many of the recent ones, have used a kinetic theory model, in which the mixing length was a sort of mean free path. In contrast, we interpret Equation (6) as a model of the momentum equation for fluid dynamics, involving structures such as waves, convective rolls, or plumes. Because it is non-local, Equation (7) allows formally stable regions to be convective, unlike MLT, because of finite velocities. This may be relevant for composition mixing in weakly stable regions, and the mass contained in convective regions.

Taking the dot product of Equation (7) with ${\boldsymbol{u}}$ gives a kinetic energy equation,

$\begin{eqnarray}&&d({u}^{2}/2)/{dt}={\boldsymbol{u}}\cdot {\boldsymbol{g}}{\beta }_{T}{\rm{\Delta }}{\rm{\nabla }}-{u}^{2}/\tau ,\end{eqnarray} \tag{ 8 }$

for which the steady-state solution¹¹ is Equation (5), with ${{\ell }}_{{\rm{d}}}={{\ell }}_{\mathrm{MLT}}^{2}/8{H}_{P}$, and ${\rm{\Delta }}{\rm{\nabla }}\gt 0$. In Equation (8), negative values of ${\rm{\Delta }}{\rm{\nabla }}$ are allowed; this permits buoyant deceleration (Brummell et al. 2002). The singularities in MLT at the convective zone boundaries (Section 9 in Gough 1977), and in boundary layers (Section 40 in Landau & Lifshitz 1959) are removed.¹²

The flow is relative to the grid of the background stellar evolution model, so the co-moving time derivative of TKE leads to

$\begin{eqnarray}&&d({u}^{2}/2)/{dt}={\partial }_{t}({\boldsymbol{u}}\cdot {\boldsymbol{u}})/2+{\boldsymbol{\nabla }}\cdot {{\boldsymbol{F}}}_{K},\end{eqnarray} \tag{ 9 }$

where ${{\boldsymbol{F}}}_{K}=\rho {\boldsymbol{u}}({\boldsymbol{u}}\cdot {\boldsymbol{u}})/2$ is a flux of kinetic energy. The generation of the divergence of a kinetic energy flux in this way is robust for dynamic models; it occurs in the more precise RANS approach (Equation (18) as well as Equation (8)).

We may write Equation (6) as

$\begin{eqnarray}&&{\partial }_{t}({\boldsymbol{u}}\cdot {\boldsymbol{u}})/2+{\boldsymbol{\nabla }}\cdot {{\boldsymbol{F}}}_{K}={\boldsymbol{u}}\cdot ({\mathcal{B}}-{\boldsymbol{D}}).\end{eqnarray} \tag{ 10 }$

In a steady state, the divergence of TKE flux is zero only if there is a local balance between the driving and the drag terms. Otherwise TKE flux may be non-negligible. The TKE flux smooths the distribution of TKE between regions in which it is generated in excess, and the whole turbulent region. The drag term is usually relatively smooth in comparison to the driving term, which can be strongly peaked. TKE transport is especially important if convection is driven by cooling near the photosphere, so that the (negative) buoyancy is localized and the stratification is strong. Meakin & Arnett (2010) have shown that stratification enhances the asymmetry in convective kinetic energy flux for driving from the top, and reduces it for driving at the bottom; see also Stein & Nordlund (1989), Cattaneo et al. (1991). This asymmetry is small for shallow convective zones, growing with stratification.

This behavior does not occur in MLT, which enforces an exact symmetry between up-flows and down-flows so that ${\boldsymbol{\nabla }}\cdot {{\boldsymbol{F}}}_{K}=0$. Although simulations of 3D atmospheres exhibit strong downward (negative) net fluxes of kinetic energy, such information was not included in MLT fits for such atmospheres (Trampedach 2007; Trampedach & Stein 2011; Magic et al. 2014). Simulations of 3D red-giant atmospheres by Kudwig & Kucinskas (2012) indicate that the fits to MLT require at least a two parameter family, as have simulations of deeper convection. In the red giant model in Viallet et al. (2013), the downward directed kinetic energy flux reaches 35% of the maximum enthalpy flux. Stein & Nordlund (1998) find that their solar model has a downward directed kinetic energy flux which is 10% of the enthalpy flux. This downward kinetic energy flux must be compensated for by a larger (outward) enthalpy flux. This kinetic energy flux is accompanied by a momentum flux, which affects the convective boundary, as shown in Section 3.8. These are nontrivial differences relative to MLT, and may have implications which are detectable with asteroseismology as deviations from the predictions of MLT models.

At present, stellar evolution theory has no turbulent heating term. This is inconsistent¹³ with Kolmogorov theory, which states that TKE is fed back into the thermal bath at the rate given by Equation (1). From the viewpoint of a dynamic model (e.g., Equation (6)), this is a "frictional" cost of moving energy by convection. Arnett et al. (2009) show that energetic self-consistency requires that the usual stellar evolution equations must be modified to include such a heating term, or equivalently, to explicitly include terms for heating by buoyancy work and divergence of kinetic energy and acoustic fluxes (see Arnett et al. 2009, Equations (20)–(22); Mocák et al. 2014, Sections 21.5, 21.6). The Kolmogorov term appears as heating in the internal energy equation and cooling (damping) in the TKE (acceleration) equation. Total energy is conserved; TKE is transformed into heat.

It may be more convenient to apply the heating term directly, rather than use the buoyancy work and divergence of TKE and acoustic fluxes, as the velocity is available from solution of Equation (7). Turbulent heating (and divergence of kinetic energy flux) may have implications for the standard solar model and solar abundances.¹⁴ Such heating may also be important for the motion of convective burning shells into electron-degenerate fuel.

In the local, steady-state, limiting case, the left-hand side of Equation (8) vanishes, and an equation similar to Equation (5) results, but with a turbulent damping length instead of a mixing length. In simulations this is the lesser of the depth of the convective zone or 4 pressure scale heights¹⁵ (Arnett & Meakin 2011b). With this change, the cubic equation of Böhm-Vitense may be derived (Smith & Arnett 2014), and we recover a form of MLT.

Had it been available, Böhm-Vitense might have identified the mixing length with the Kolmogorov damping length (Equation (1)). However, Kolmogorov found the damping length ℓ_d to be the depth of the turbulent region, so that it is not a free parameter, unlike MLT. There is a further issue: ${\epsilon }_{t}$ is the average dissipation rate, not the instantaneous local value (${u}^{3}/{{\ell }}_{{\rm{d}}}$) which fluctuates over time and space (see Figure 4 in Meakin & Arnett 2007b); that is, $u\ne v$ except on average. This is reminiscent of the RANS approach (Sections 2.6 and 2.7).

Suppose we assume that the integral scale motion is that of a 2D convective roll, where $d{\boldsymbol{u}}/{dt}$ is given by Equation (7). Using this and a corresponding thermal energy equation, we obtain a form of the classic Lorenz equations, but with a nonlinear damping term provided by the Kolmogorov cascade (Arnett & Meakin 2011b). Because of the time lag, as implied by the time needed to traverse the cascade from integral to Kolmogorov scales, the modified equations are even more unstable than the original ones, and have chaotic behavior.¹⁶

In Equation (7) it was assumed that the density fluctuation which drives the buoyancy could be represented by $-g(\delta \rho /\rho )\approx {\boldsymbol{g}}{\beta }_{T}{\rm{\Delta }}{\rm{\nabla }}$, involving only a fluctuation in temperature. This is only true for uniform composition and mild stratification. The formulation makes use of the expansion of pressure fluctuation,

$\begin{eqnarray}P^{\prime} & = & {(\partial P/\partial T)}_{\rho ,Y}T^{\prime} +{(\partial P/\partial \rho )}_{T,Y}{\rho }^{\prime }\\ & & +{(\partial P/\partial Y)}_{T,\rho }Y^{\prime} ,\end{eqnarray} \tag{ 11 }$

which may be written as

$\begin{eqnarray}{\rho }^{\prime }/\rho & = & -{\beta }_{T}(T^{\prime} /T)-{\beta }_{Y}(Y^{\prime} /Y)\\ & & +(P/\rho {s}^{2})P^{\prime} /P,\end{eqnarray} \tag{ 12 }$

where

$\begin{eqnarray}&&{\beta }_{T}=-{(\partial \mathrm{ln}\rho /\partial \mathrm{ln}T)}_{Y,P},\end{eqnarray} \tag{ 13 }$

$\begin{eqnarray}&&{\beta }_{Y}=-{(\partial \mathrm{ln}\rho /\partial \mathrm{ln}Y)}_{T,P},\end{eqnarray} \tag{ 14 }$

$\begin{eqnarray}&&{s}^{2}={(\partial P/\partial \rho )}_{T,Y}.\end{eqnarray} \tag{ 15 }$

Here s is the sound speed. The composition variable Y denotes the number of free particles per baryon (Arnett 1996), and is essentially the inverse of the mean molecular weight μ (Kippenhahn & Weigert 1990; Hansen & Kawaler 1994). An illustrative and simple example is the ideal gas, $P={\mathcal{R}}\rho {YT}$. For subsonic flows, $| P^{\prime} /P| \sim {(u/s)}^{2}$, where $u/s$ is the Mach number of the flow, and is small.¹⁷ In MLT, the pressure fluctuation is assumed zero (no acceleration by pressure dilatation), so

$\begin{eqnarray}&&{\beta }_{T}(T^{\prime} /T)+{\beta }_{Y}(Y^{\prime} /Y)\approx -({\rho }^{\prime }/\rho ),\end{eqnarray} \tag{ 16 }$

and it is further assumed that $Y^{\prime} =0$ to obtain Equation (5). Even in the limit of negligible pressure fluctuations, variations in Y enter in a way similar to variations in T, so even small composition variations can be significant when superadiabatic temperature variations are also small. Many of the difficulties found using MLT are related to situations in which $Y^{\prime} \ne 0$: overshooting, semi-convection, and entrainment.

There seems to be a deep connection between Equation (7), RTIs, and turbulent mixing. An almost identical equation (Equation (4.1) in Abarzhi 2010) is used to describe the nonlinear development of the RTI into the turbulent mixing regime. Unlike canonical Kolmogorv turbulence, the RT turbulent mixing is statistically unsteady, and involves the transport of potential and kinetic energies as well as enthalpy. Because of its importance in a variety of high energy–density conditions (Kane et al. 1997; Remington et al. 1999, 2006; Zeldovich & Razier 2002; Dimonte et al. 2004; Drake 2009; Kuranz et al. 2011), much experimental effort for its study as well as an extensive literature have developed.

The RTI happens when a heavier fluid overlays a lighter one, proceeding from linear instability of perturbations (Chandrasekhar 1961), to mildly nonlinear motion of bubbles and spikes, and then to nonlinear turbulent mixing (Abarzhi 2010). The initial acceleration is 1D, but as instability develops, the motion breaks symmetry and approaches isotropy (as seen in a co-moving frame), much like the cascade in steady turbulence (Frisch 1995). The essential difference between stellar convection and RTI is that the RTI is not contained, while convection operates within a definite and slowly varying volume. This means that the vertical and the transverse scales are causally connected in convection, but may be independent in the RTI (Abarzhi 2010).

Inconsistency between experimental and numerical investigation of the RTI in the nonlinear regime led to the ${\alpha }_{{\rm{b}}}$ problem (Dimonte et al. 2004). The RTI in the limit of strong mode-coupling can be initiated to have self-similar evolution, so that the amplitude (diameter of the bubble ${D}_{{\rm{b}}}\propto {h}_{{\rm{b}}}$) evolves as ${h}_{{\rm{b}}}\sim {\alpha }_{{\rm{b}}}{\rm{A}}{{gt}}^{2}$, where A is the Atwood number (density ratio, Chandrasekhar 1961), g is gravity, and t the elapsed time. The simulation value ${\alpha }_{{\rm{b}}}\sim 0.025\pm 0.003$ is smaller than the experimental value ${\alpha }_{{\rm{b}}}\sim 0.057\pm 0.008$. This discrepancy seems to have been resolved by the idea that unquantified errors in the experimental initial conditions were the cause. To the extent that such uncertainties cannot be precisely known, this suggests a statistical approach, and illustrates the need for combined theoretical, experimental, and numerical studies.

Meakin & Arnett (2007b) found that regions of their simulated convection zone exhibited recurring "bursts" of convection (see their Figure 4). These bursts, although multi-modal ($n\sim 5$), seem to share the chaotic behavior of the Lorenz (1963) model of a single-mode convective roll (Arnett & Meakin 2011b). This encourages the use of Equation (7), which is related to the momentum-driven model of RTI (Abarzhi 2010), for timescales less than or of order of the transit time. For longer, evolutionary timescales (stellar convection) we need to average over fluctuations, which means averaging over several transit times for the convective roll (see Equation (18) below). These bursts result from underlying physics similar to that in the RTI; their short timescale behavior may be relevant for stellar pulsations and eruptions (the τ-mechanism, Arnett & Meakin 2011b, or equivalently, stochastic excitation of oscillations, Goldreich & Kumar 1994; Goldreich et al. 1994; Aerts et al. 2010).

The weak coupling between driving at the large scale, and dissipation at the small scale, allows time dependent fluctuations of significant amplitude in luminosity and turbulent velocity. The term $\partial {\boldsymbol{u}}/\partial t$ (Equations (6) and (17)) is needed for chaotic fluctuations and wave generation. These fluctuations have a cellular structure in space and time; if there are many cells, with random phases, the fluctuations in the average total luminosity are reduced by cancellation (Arnett & Meakin 2011b).

Fluctuations are fundamental features of turbulence and mixing. Because of sensitivity to initial conditions which can never be known with complete accuracy, descriptions of turbulence should be statistical in nature, even though the equations are deterministic (Frisch 1995). Turbulent simulations can be said to be numerically converged only in a statistical sense. Eventually trajectories will diverge. Lyapanov exponents characterize this divergence, a feature characteristic of turbulence (Manneville 2010) that makes turbulent mixing so effective. Unlike the diffusion picture, in which a stellar mixing front moves radially, limited by the random walk of mean-free-path strides, turbulent mixing involves a network of trajectories throughout the space of the turbulent region, laced with inhomogeneities, which finally disappear at the Kolmogorov scale.

In stratified regions, mass conservation constrains the flow, but it tends to change the cross-sectional area of the plumes as opposed to limiting their range. Although the flow is locally wild with fluctuations, these tend to cancel upon horizontal and time averaging, leaving a much more placid behavior due to the cancellation of random phases. Figure 2 illustrates this for a particular but representative case; the velocity in the theta direction, ${v}_{\theta }$, is shown as a function of radius, from the oxygen burning data set in Viallet et al. (2013). The top panel shows the instantaneous value of ${v}_{\theta }$ (in units of cm s⁻¹) for a sequence of time steps $\delta t\sim 0.5\;{\rm{s}}$. The bottom panel shows the running average (a horizontal average, i.e., over a spherical surface of radius r) of the same variable over 300 such time steps (150 s), stepping forward over 20 time steps (10 s) at a stride, on the same velocity scale. The amplitude in the bottom panel is much reduced by cancellation; what does remain is the larger length scale, as suggested by the cascade idea discussed in Section 2.1. The cancellation does not work for quadratic terms; they remain non-zero, e.g., contributing to the rms velocity in this case (see Section 2.7). The product of fluctuations in velocity and temperature give rise to the enthalpy flux; those in velocity and composition give rise to the composition flux.

Zoom In Zoom Out Reset image size

A stellar evolution code must step over the shorter turnover time scales (weather) to solve for the evolutionary times (climate). How can this be done? It requires an average over active and inactive cells. The steady-state limit of the Lorenz equation seems to give a reasonable approximation to its average behavior, filtering out the chaotic fluctuations (Arnett & Meakin 2011b). Instead of $d{\boldsymbol{u}}/{dt}=0$, we use

$\begin{eqnarray}d{\boldsymbol{u}}/{dt} & = & \partial {\boldsymbol{u}}/\partial {\text{}}t+({\boldsymbol{u}}\cdot {\boldsymbol{\nabla }}){\boldsymbol{u}},\\ & & \to ({\boldsymbol{u}}\cdot {\boldsymbol{\nabla }}){\boldsymbol{u}}.\end{eqnarray} \tag{ 17 }$

We apply the same approximation (Equation (17)) to Equation (7) for slow stages of stellar evolution. This allows non-local behavior, will prove important for our discussion of convective boundaries later in Section 3, and can represent ram pressure (Reynolds stress) and the flux of TKE; see also Section 3.2 in Porter & Woodward (2000), for a discussion of ram pressure in 3D simulations relative to MLT.

Now we have established connections between an acceleration equation (Equation (6)) and (1) MLT, (2) historical attempts to extend MLT, (3) modern research on RTI (Abarzhi 2010), (4) the important advances of Kolmogorov (1941, 1962) and Lorenz 1963, and (5) a rational way to step over fluctuations for stellar evolution.

A more rigorous alternative is to use the RANS approach, which directly averages the fluctuations over space and time. This has been explored by Canuto (Canuto 2012a, 2012b, 2012c, 2012d, 2012e), see also Xiong et al. (1997), Deng et al. (2006); a detailed comparison with their work, while desirable, is beyond the scope of this paper. Canuto uses simulations and experiments from geophysics to effect a closure of the RANS equations, while in contrast, our closure of the RANS is based on our 3D simulations.

The TKE equation is obtained by a Reynolds decomposition of the velocity, density, and pressure (detailed discussion may be found in Meakin & Arnett 2007b; Arnett et al. 2009; Viallet et al. 2013; Mocák et al. 2014). In principle the TKE is exact; errors arise from closure, i.e., our analytical approximations to the terms in the RANS equations are at fault. Well-resolved 3D ILES simulations show excellent agreement with the TKE (Viallet et al. 2013), and allow the dominant terms to be identified. Being more general than the simpler approximations discussed above, the TKE allows us to identify and quantify neglected terms. Most importantly, it allows an enormous simplification and compaction of the 3D numerical data, while that data in turn allows a closure of the RANS procedure.

The TKE may be written as (Meakin & Arnett 2007b):

$\begin{eqnarray}&&{\partial }_{t}\langle \bar{\rho {E}_{K}}\rangle +{\boldsymbol{\nabla }}\cdot \langle \bar{\rho {E}_{K}{{\boldsymbol{u}}}_{0}}\rangle \\ &&\quad =-{\boldsymbol{\nabla }}\cdot \langle \bar{{{\boldsymbol{F}}}_{P}+{{\boldsymbol{F}}}_{K}}\rangle +\langle \bar{P^{\prime} {\boldsymbol{\nabla }}\cdot u^{\prime} }\rangle \\ &&\qquad +\langle \overline{{\rho }^{\prime }{\boldsymbol{g}}\cdot {\boldsymbol{u}}^{\prime} }\rangle -\rho {\epsilon }_{{\rm{d}}},\end{eqnarray} \tag{ 18 }$

We use $\langle q\rangle$ and $\bar{q}$ to denote angular and time averages of a quantity q. Primes refer to fluctuating quantities; for example ${\boldsymbol{u}}={{\boldsymbol{u}}}_{0}+{\boldsymbol{u}}^{\prime}$, and $\langle {\boldsymbol{u}}\rangle ={{\boldsymbol{u}}}_{0}$, and similarly for the time average. The TKE per unit mass is ${E}_{K}=\frac{1}{2}({\boldsymbol{u}}^{\prime} \;{\boldsymbol{\cdot }}\;{\boldsymbol{u}}^{\prime} )$, a measure of the rms turbulent velocity. The acoustic and turbulent kinetic fluxes are ${{\boldsymbol{F}}}_{P}=P^{\prime} {\boldsymbol{u}}^{\prime}$ and ${{\boldsymbol{F}}}_{K}=\rho {E}_{K}{\boldsymbol{u}}^{\prime}$. The dissipation may be written as

$\begin{eqnarray}&&{\epsilon }_{{\rm{d}}}=\langle \bar{{\boldsymbol{u}}^{\prime} \cdot {\boldsymbol{u}}^{\prime} | u^{\prime} | }\rangle /{\ell },\end{eqnarray} \tag{ 19 }$

a form which we identify with Equation (1), the expression of Kolmogorov (1941, 1962); notice that it involves averages of powers of the velocity fluctuation, not the instantaneous values.

Using the RANS approach is equivalent to using the bottom panel in Figure 2 rather than the top; it removes the fluctuating activity which cancels (has no net effect), while keeping what does not cancel.

To better understand the implications of the TKE, consider (1) a steady state (${\partial }_{t}\langle \bar{\rho {E}_{K}}\rangle =0$) with (2) no background motion (${{\boldsymbol{u}}}_{0}=0$). Then the TKE reduces to the divergence of the fluxes ${\boldsymbol{\nabla }}\cdot \langle \bar{{{\boldsymbol{F}}}_{P}+{{\boldsymbol{F}}}_{K}}\rangle ,$ balancing the net result of two source terms $\langle \bar{P^{\prime} {\boldsymbol{\nabla }}\cdot {\boldsymbol{u}}^{\prime} }\rangle$ and $\langle \overline{{\rho }^{\prime }{\boldsymbol{g}}\cdot {\boldsymbol{u}}^{\prime} }\rangle$, and a damping term $-\rho {\epsilon }_{{\rm{d}}}$:

$\begin{eqnarray}{\boldsymbol{\nabla }}\cdot \langle \bar{{{\boldsymbol{F}}}_{P}+{{\boldsymbol{F}}}_{K}}\rangle & = & \langle \bar{P^{\prime} {\boldsymbol{\nabla }}\cdot {\boldsymbol{u}}^{\prime} }\rangle \\ & & +\langle \overline{{\rho }^{\prime }{\boldsymbol{g}}\cdot {\boldsymbol{u}}^{\prime} }\rangle -\rho {\epsilon }_{{\rm{d}}}.\end{eqnarray} \tag{ 20 }$

This may be integrated over the convection zone (taking the surface fluxes to be zero or small at the boundaries), and if we ignore the pressure dilatation $\langle \bar{P^{\prime} {\boldsymbol{\nabla }}\cdot {\boldsymbol{u}}^{\prime} }$ for the moment, gives an expression for the damping length ℓ_d,

$\begin{eqnarray}&&{{\ell }}_{{\rm{d}}}={\displaystyle \int }_{\mathrm{CZ}}\langle \bar{{({\boldsymbol{u}}^{\prime} \cdot {\boldsymbol{u}}^{\prime} )}^{\displaystyle \frac{3}{2}}}\rangle {dm}/\left({\displaystyle \int }_{\mathrm{CZ}}\langle \overline{\displaystyle \frac{\rho ^{\prime} }{{\rho }_{0}}{\boldsymbol{g}}\cdot {\boldsymbol{u}}^{\prime} }\rangle \right){dm},\end{eqnarray} \tag{ 21 }$

which is a global condition that must be satisfied to be consistent with Kolmogorov damping, which also requires that ℓ_d is approximately the depth of the turbulent region. This characteristic length scale is a fundamental property of turbulence, and is generated robustly in the numerical simulations.

Equation (21) might be regarded as a generalization of the Roxburgh (1989, 1992) integral constraint to include damping by turbulence. Notice that ℓ_d, which appears in both Equations (7) and (21), must be solved for consistently; it tends to be a slowly varying function, of order of the convective zone depth. Equation (21) involves some of the important "bulk" properties discussed by Canuto (1992), and is a statement of a global balance between driving and damping.

What approximations would be necessary to make the TKE equation equivalent to MLT? In MLT, (1) the net flux of TKE ${{\boldsymbol{F}}}_{K}$ is defined to be zero by symmetry, (2) pressure fluctuations are ignored so the acoustic flux ${{\boldsymbol{F}}}_{P}$ and pressure dilatation $\langle \bar{P^{\prime} {\boldsymbol{\nabla }}\cdot {\boldsymbol{u}}^{\prime} }\rangle$ are zero, and (3) the damping length ℓ_d is taken to be an arbitrary adjustable parameter. Enforcing these gives

$\begin{eqnarray}&&\langle \overline{{\rho }^{\prime }{\boldsymbol{g}}\cdot {\boldsymbol{u}}^{\prime} }\rangle =\rho {\epsilon }_{{\rm{d}}}.\end{eqnarray} \tag{ 22 }$

This is the local version of the global balance in Equation (21); it is equivalent to the Böhm-Vitense cubic equation of MLT for the appropriate choice of mixing length (Smith & Arnett 2014).

This approximation leads to a series of errors: (1) Symmetry between up-flows and down-flows is broken by stratification, so that TKE fluxes are not generally zero (Stein & Nordlund 1989; Cattaneo et al. 1991; Canuto 1992). This is a qualitative error. (2) Pressure fluctuations may not be ignored for strongly stratified convection zones. This is a quantitative error. Viallet et al. (2013) find that acceleration by the pressure dilatation term is comparable to that from buoyancy. (3) The damping length may not be freely adjusted if the relation of Kolmogorov (1941, 1962) is to be satisfied. Such adjustments are usually necessary to compensate for a lack of non-locality in atmospheres due to the lack of ram pressure, and deeper into interiors due to a lack of kinetic energy flux (the two parameters discussed in regard to 3D atmospheres in Section 2.3).

Our efforts have been three-fold: (1) construction of accurate numerical solutions of the Navier–Stokes equations which exhibit turbulence, (2) theoretical analysis of these solutions in the RANS framework to determine the most important features, and (3) invention of simpler analytic representations which capture the essential features of the numerical solutions. Pasetto et al. (2014) have presented a novel analytical theory of convection in stars which does not contain a mixing-length parameter; this is an alternative to (3) above, and it is of interest to compare how well it agrees with both our numerical solutions (1 and 2), and our analytic approximations (3).

As we have shown in Section 2.7, the natural length scale for convection is the dissipation length for the turbulent cascade. Part of the foundation of the model of Pasetto et al. (2014) is the use of potential flow and the Bernoulli equation (Landau & Lifshitz 1959, Equation (10.7) in Section 10), which result from the Euler equation, not the Navier–Stokes equation. Their theory seems to be equivalent to assuming the process occurs on a scale much less than the size of the convective region, so that there is no way to define a length scale for turbulent dissipation. In contrast, following Kolmogorov (Section 2.1), the length scale in our theory is the size of the turbulent region, which is not arbitrary but determined by the turbulent flow. Our length scale is not an assumption (as in MLT) but a consistent and robust result of our simulations. It is the length scale over which driving and damping of turbulence balance (Section 2.3). In order to describe the turbulent cascade, a complete theory must deal with the whole turbulent region.

Is the theory of Pasetto et al. (2014) physically correct? Stellar interior convection is extremely turbulent, so the question becomes: what are the errors introduced by ignoring turbulence? Landau & Lifshitz (1959) give a careful discussion of the applicability of potential flow (their Section 9), and they note that the validity of Bernoulli's equation is limited because of the formation of boundary layers in which viscous effects must be included (see also Prandtl & Tietjens 1934). Stars have large Reynolds numbers, so that turbulent boundary layers form (Landau & Lifshitz 1959, Chap. III), as they do in our simulations (Figure 3). The Pasetto theory, like MLT, ignores boundary layers and turbulence, as well as composition gradients.

Zoom In Zoom Out Reset image size

A basic assumption of the Pasetto et al. (2014) theory is that velocities of lateral expansion are much larger than those of the vertical rise of convective elements (their Section 4.2). However, the simulations show average velocities in the turbulent region which are not strongly biased toward the laterial directions; this was already clear in Meakin & Arnett (2007b; their Figure 6), and has held true for subsequent simulations with refined resolution (Viallet et al. 2013; S. Campbell et al. 2015, in preparation). The rms velocity in the radial direction is actually larger than the lateral rms velocities, rather than smaller (Arnett et al. 2009).

A key test presented in Pasetto et al. (2014) of their theory is a comparison with MLT¹⁸ at $r=0.98{R}_{\odot }$, well inside the super-adiabatic region (SAR) at $r\sim 0.9985{R}_{\odot }$ in the Sun. It is the inefficient convection in the SAR which determines the solar radius in calibrations of stellar evolutionary codes, so that a test in the SAR would be instructive. Pasetto et al. (2014) state "Convective elements in this region have low thermal capacity, so that the super-adiabatic approximation can no longer be applied, and the temperature gradient of the elements and surrounding medium must be determined separately." The theory in its present form may not yet be applicable to the SAR.

The value of the Pasetto theory may prove to lie in its significant conceptual differences from MLT, and in its use as a null case to provide insight into the effects of turbulence.

For the oxygen-burning shell, the temperature T has an abrupt jump inside the mixing region (radius $r\sim 4.3\times {10}^{8}\ \mathrm{cm}$ in Figure 3). Pressure is continuous through the boundary containing this transition, so that the density curve has a corresponding dip; see Figure 2 in Meakin & Arnett (2007b) or Figure 5 in Viallet et al. (2013). This implies a steep increase in entropy; as evolution continues this entropy jump grows, and the transition region narrows. Such steep gradients in T are a consequence of cooling by neutrinos. They are not seen in earlier, photon-cooled stages of evolution and can only be supported for times short compared with timescales for thermal diffusion and electron heat conduction. This is easily the case for oxygen burning because of high opacity and short evolutionary times ($\sim {10}^{5}\;{\rm{s}}$).

The Péclet number is defined as the ratio of the advective transport rate to the diffusive transport rate of the physical quantity being transported, which here we take to be thermal energy, so

$\begin{eqnarray*}&&\mathrm{Pe}=\displaystyle \frac{\mathrm{thermal}\ \mathrm{advection}\ \mathrm{rate}}{\mathrm{thermal}\ \mathrm{diffusion}\ \mathrm{rate}}.\end{eqnarray*}$

In oxygen burning, radiative diffusion is slow while advection occurs rapidly, giving large Péclet numbers (formally infinite since radiative diffusion was small enough to be neglected in some simulations; the infinity results from the denominator in the definition being a negligible term, not from any exceptional behavior of the physics).

This contrasts with the situation in stellar atmospheres, in which the radiative diffusion becomes faster than advective transport, so that $\mathrm{Pe}\lt 1$. This difference in Péclet numbers suggests the possibility of a fundamental flaw in the notion that observations of stellar atmospheres may be sufficient to define the nature of deep stellar convection. See the discussion in Zahn (1991) and Viallet et al. (2015).

Peter Eggleton took an early step in dealing with steep gradients in composition, with the introduction of a diffusion operator which he stressed was ad hoc (Eggleton 1973). This numerically advantageous procedure has been widely adopted for stellar evolution, even though it has the potentially worrisome mathematical property that it increases the order of the spatial derivatives in the equations to be solved. The Eggleton (1973) equation is

$\begin{eqnarray}&&\displaystyle \frac{d}{{dm}}\left(\sigma \displaystyle \frac{{dX}}{{dm}}\right)=\displaystyle \frac{{DX}}{{Dt}}-{\mathcal{R}}\end{eqnarray} \tag{ 23 }$

where X is the mass fraction, m is the lagrangian mass coordinate, $\sigma ={v}_{\mathrm{ML}}{{\ell }}_{\mathrm{ML}}{(4\pi {r}^{2}\rho )}^{2}$ is the effective diffusion coefficient, and ${\mathcal{R}}$ is the nuclear reaction network matrix (Arnett 1996). This is equivalent to modeling convection as "turbulent diffusion." The left-hand side is the heuristic diffusion operator; the right hand side is the reaction network operator. The actual composition flux is related to the co-moving derivative on the right-hand side; see Arnett (1996), Section 4.6. Eggleton integrates over the convection zone to eliminate that spatial derivative; usually it is simply ignored in stellar codes.

The Eggleton approach is equivalent to approximating the composition flux

$\begin{eqnarray}&&{F}_{Y}=\rho {\mathcal{A}}\overline{\langle u^{\prime} Y^{\prime} \rangle }\end{eqnarray} \tag{ 24 }$

by a "down-gradient" expression (critically discussed by Canuto 1992),

$\begin{eqnarray}&&{F}_{Y}\to \rho {\mathcal{A}}\ u(-{\ell }\partial Y/\partial r).\end{eqnarray} \tag{ 25 }$

Direct comparison with simulations shows that this can be qualitatively wrong (by two orders of magnitude). For a contact discontinuity (Landau & Lifshitz 1959, Section 81), ${F}_{Y}\to \rho {\mathcal{A}}u{\rm{\Delta }}Y$, as in Equation (24), not $\rho {\mathcal{A}}u(-{\ell }/{\rm{\Delta }}r){\rm{\Delta }}Y\to \infty$, as in Equation (25). Proper scaling requires that ${\ell }\to {\rm{\Delta }}r$ at a boundary if Equation (25) is used.

As Eggleton intended, the algorithm smooths steep gradients, but sometimes faster than real physical processes do, as Eggleton warned. To the extent that gradients in abundance need to be correctly represented (e.g., for ionic diffusion, or density structure), the down-gradient approximation (in Equations (23) and (25)), is questionable. In particular, fluxes directly computed in simulations (Meakin & Arnett 2007b; Viallet et al. 2013) show that the down-gradient approximation fails in boundary layers (Mocák et al. 2014).

While real atomic (ionic) diffusion is thought to be slow in stars, the diffusion operator is second order in space derivatives, so that it becomes important in steep composition gradients, i.e., boundaries. Georges Michaud has led in the application of true diffusion processes and radiative levitation to stellar evolution (Michaud 1970, 1991; Michaud et al. 2007). Recently these processes have been applied to horizontal branch and subdwarf B (sdB) stars (Michaud et al. 2005, 2007, 2011; Hu et al. 2008, 2009, 2010, 2011; Bloemen et al. 2014). Gravitational settling (Hu et al. 2009) and radiative levitation (Hu et al. 2011) are important to (1) recover the iron-group opacity bump that excites the pulsations (Charpinet et al. 1997) in those stars, (2) obtain the correct position of the instability strip in the $\mathrm{log}g-{T}_{\mathrm{eff}}$ diagram, and (3) help in understanding their observed atmospheric abundances (Michaud et al. 2011).

Because the Eggleton (1973) diffusion uses a difference operator similar to that for ionic diffusion (second order in space), and may reduce the gradients which drive that diffusion, care should be taken that the algorithmic diffusion does not cause errors in the real diffusion (e.g., see Schindler et al. 2015).

In stellar physics, the idea of semi-convection has spawned various algorithms (e.g., Schwarzschild & Härm 1958; Stothers 1963; Castellani et al. 1971a, 1971b; Demarque & Mengel 1972; Sweigart & Gross 1976; Dorman & Rood 1993), some of which seem to be physically and numerically inconsistent with others. The term "semi-convection" refers to a mixing process which occurs in a region that is stable according to the Ledoux criterion but unstable according to the Schwarzschild criterion. It generally is thought to involve mixing of composition, but not significant enthalpy. The composition profile may be adjusted to marginal stability according to the Ledoux criterion.

Semi-convection is also often discussed as a double diffusive instability, involving an interaction between radiative diffusion and ionic diffusion (Spruit 2013; Lattanzio et al. 2014). Although both radiative and ionic diffusion may be included in a 1D stellar code, this does not capture their interaction and 3D dynamics. Semi-convection may be related to oceanic phenomena (thermohaline mixing) in which heat flow and salt concentration play the doubly diffusive roles, and which have a long history of study (e.g., see Chap. 8 in Turner 1973; Gill 1982). Rosenblum et al. (2011) and Wood et al. (2013) give an extensive discussion with numerical simulations based on the oceanic model, and conclude that, while the problem can be solved in the planetary range of parameter space, the stellar case requires a large extrapolation. This difficulty may be further exacerbated by the indication that many such regions in stars are bathed in a flux of g-mode waves (Meakin & Arnett 2007b), which are a nonlocal effect that may complicate the analysis in a nontrivial way (Mocák et al. 2014).

Even with these uncertainties, there are energetic constraints (see Equation (30)) which must be obeyed. The amount of mixing possible is limited by the energy available to mix, which is generally taken to be related to the excess ${{\rm{\nabla }}}_{r}-{{\rm{\nabla }}}_{a}$, so that luminosity is used to supply the energy required to mix.

MLT, as a local theory, must be supplemented by additional assumptions about behavior at the boundaries of the convection zone (Spiegel 1971, 1972). These are usually discussed in terms of linear stability theory, i.e., in terms of the Ledoux and the Schwarzschild criteria (Kippenhahn & Weigert 1990) being positive. The Schwarzschild criterion for convective instability is defined by

$\begin{eqnarray}&&{\mathcal{S}}={{\rm{\nabla }}}_{r}-{{\rm{\nabla }}}_{a}.\end{eqnarray} \tag{ 26 }$

Here ∇_r is what the dimensionless temperature gradient would be if all the luminosity were carried by radiative diffusion and ∇_a is the adiabatic gradient (see the Appendix). The Ledoux criterion for convective instability has a composition dependence, and is defined by

$\begin{eqnarray}&&{\mathcal{L}}={{\rm{\nabla }}}_{r}-{{\rm{\nabla }}}_{a}-\displaystyle \frac{{\beta }_{Y}}{{\beta }_{T}}{{\rm{\nabla }}}_{Y}.\end{eqnarray} \tag{ 27 }$

The last term is written as $\frac{\phi }{\delta }{{\rm{\nabla }}}_{\mu }$ by Kippenhahn & Weigert (1990), Section 6.1, their Equation (6.12). The β factors are defined as in Section 2.4 above. Notice that positive $\frac{{\beta }_{Y}}{{\beta }_{T}}{{\rm{\nabla }}}_{Y}$ and positive $\frac{\phi }{\delta }{{\rm{\nabla }}}_{\mu }$ both inhibit mixing.

Neither of these choices seems satisfactory. They have no dependence upon the vigor of the flow on the unstable side of the boundary, which clearly must make a difference.

Linear perturbation theory examines the instability of a stable region, treating both sides of the boundary equally. In reality they differ: one side is convective. The stiffness of the non-convective side is measured by the Brunt-Väisälä (buoyancy) frequency N, (see Equation (6.18) in Kippenhahn & Weigert 1990, and Equation (3.73) in Aerts et al. 2010), where

$\begin{eqnarray}&&{N}^{2}=-\displaystyle \frac{{\beta }_{T}g}{{H}_{P}}({{\rm{\nabla }}}_{e}-{{\rm{\nabla }}}_{a}-({\beta }_{Y}/{\beta }_{T}){{\rm{\nabla }}}_{Y})=-\displaystyle \frac{{\beta }_{T}g}{{H}_{P}}{\mathcal{L}}.\end{eqnarray} \tag{ 28 }$

N is the frequency of elastic rebound from a perturbation; it is imaginary in convective regions. Here ∇_e is the dimensionless temperature gradient relevant¹⁹ to the perturbed element. On the non-convective side of the boundary, it may be the same as ∇_r above, giving the second equality, which refers to the tendency to restore stability in the radiative region.

A delicate point is the value of ${{\rm{\nabla }}}_{Y}$ near the boundary (Gabriel et al. 2014). By what mechanism does mixing occur? What is the structure of the partially mixed region of transition between well-mixed and unmixed? Present practice in stellar evolution is to use the Schwarzschild criterion, which has no ${{\rm{\nabla }}}_{Y}$, so that these issues may be ignored, or to use the Ledoux criterion with one of the prescriptions for semi-convective mixing (see Section 3.4).

Such interfacial issues have long been studied in the fluid dynamics and geophysics communities; see Turner (1973) for an extensive discussion. The Richardson number is defined as some measure of

$\begin{eqnarray*}&&\mathrm{Ri}=\displaystyle \frac{\mathrm{potential}\ \mathrm{energy}\ \mathrm{needed}\ \mathrm{to}\ \mathrm{mix}}{\mathrm{kinetic}\ \mathrm{energy}\ \mathrm{available}\ \mathrm{to}\ \mathrm{mix}}.\end{eqnarray*}$

The linear condition for ability of a layer to resist shear is the "gradient" Richardson number Ri.

$\begin{eqnarray}&&\mathrm{Ri}={N}^{2}/{(\partial u/\partial r)}^{2}\gt \displaystyle \frac{1}{4}\end{eqnarray} \tag{ 29 }$

is stable; larger stiffness (N²) and less swirling (${(\partial u/\partial r)}^{2}$) tend toward stability. In their discussion of entrainment, Meakin & Arnett (2007b) used a "bulk" (i.e., non-local and nonlinear) Richardson number which involved an integral over the region around the boundary.

In the absence of global rotation, a layer having constant total entropy²⁰ is energetically neutral with regard to mixing. If after a mixing episode, the luminosity returns to its value for radiative balance (∇_r is unchanged), then the additional energy²¹ required to remove the stable compositional stratification is

$\begin{eqnarray}&&{E}_{\mathrm{mix}}={{gH}}_{P}{\beta }_{T}({\mathcal{L}}-S)={{gH}}_{P}(-{\beta }_{Y}{{\rm{\nabla }}}_{Y}).\end{eqnarray} \tag{ 30 }$

Both ${\beta }_{T}$ and ${\beta }_{Y}$ are intrinsically negative in stars. If this energy E_mix changes sign, mixing may occur which is driven by the gradient in composition (Mocák et al. 2011a). Using a specific kinetic energy of $\frac{1}{2}{u}^{2}$, a Richardson number may be constructed,

$\begin{eqnarray}&&\mathrm{Ri}=2\;{{gH}}_{P}(-{\beta }_{Y}{{\rm{\nabla }}}_{Y})/{u}^{2}.\end{eqnarray} \tag{ 31 }$

Here the traditional $\mathrm{Ri}\gt 1/4$ is a plausible condition for stability, at least roughly.

In their pioneering work on solar convection, Stein & Nordlund (1989) carefully explored the topology of convective flow below the photosphere: converging, cool downdrafts being dominant, with radiative cooling providing the entropy deficit which drives the circulation. Freytag et al. (1996) examined shallow (weakly stratified) convection, driven by atmospheric cooling, and emphasized the importance of the atmosphere in determining the nature of the convection zone. As deep interior convection (Arnett 1994; Bazàn & Arnett 1994) has no atmosphere, atmospheric physics can have no strong role there (the circulation is driven by nuclear burning). Furthermore, the bottom boundary, which could be ignored in the simulations of Stein & Nordlund (1989), may be important for the detailed effects of solar convection on the interior.

Schwarzschild (1958), Section 11, showed that, for stellar interior models, the atmosphere could be represented by an entropy jump between the photosphere and the adiabatic (deep) convective region. This entropy jump is a primary parameter for determining the depth of the convection zone. The atmospheric model is crucial for predicting spectral features for a given entropy jump, but has a weak influence on that entropy jump itself (Tanner et al. 2012, 2014).

Many features of the atmospheric and deep interior simulations are similar, leading to the idea that atmospheric physics, however crucial for spectral formation (Stein & Nordlund 1998; Magic et al. 2013, 2014), may be treated as a boundary condition issue rather than a key feature of deep turbulent convection. Meakin & Arnett (2010) showed that the general characteristics of the flow in solar convection (narrow, fast down-flows with broad, slow up-flows and acceleration by pressure dilatation; Stein & Nordlund 1989; Viallet et al. 2013), require only localized top cooling and stratification. Global simulations of the solar convection zone are necessarily less well resolved for comparable computational resources; the simulations of Miesch et al. (2007) are beginning to show turbulence, but may require finer zoning to deal with some details of the turbulent flow (e.g., Hanasoge et al. 2012; Brandenburg 2015).

The simplest of stellar convection zones are cooled by the local processes (cooling by neutrino emission and heating by nuclear burning), rather than the non-local processes (radiative transfer), giving a cleaner example of the dynamics of boundaries for deep convection. A slightly more complex case is a convection zone with heat conduction by radiative diffusion; Viallet et al. (2013) consider both. These two cases cover almost all of the conditions relevant to stellar evolution, except the outer layers simulated in 3D atmospheres.

For the oxygen-burning shell, some integral properties of the main convective region and the braking layers are summarized in Table 2. About 14% of the mass and 15% of the thickness of the total convection zone are in the boundary layers (upper BL and lower BL), as is 8.5% of the TKE. These boundary regions provide deceleration (braking) of the vertically directed flow, allowing it to remain bounded by the convective volume. If the buoyancy flux is $q=-g\langle {u}_{z}^{\prime }\rho ^{\prime} \rangle /{\rho }_{0}$, then the rate at which TKE increases due to buoyancy in a region a, is

$\begin{eqnarray}&&{L}_{\mathrm{buoy}}(a)={\displaystyle \int }_{a}q\ {dm},\end{eqnarray} \tag{ 32 }$

which is positive in the middle region, but negative in the boundary regions. These regions of negative buoyancy are a robust qualitative feature of the simulations, dating back to early 2D work (Hurlburt et al. 1984; Arnett 1994). In the oxygen-burning shell they reduce the driving of TKE by only 1.8%.

Table 2.  Integral Properties of Convection Zone Regions

Variable	Symbol	Total (CZ+BL+BL)	Lower BL	Upper BL
Mass	${\rm{\Delta }}m/{M}_{\odot }$	0.9205	0.0161	0.1150
Depth	${\rm{\Delta }}r/{10}^{8}\;\mathrm{cm}$	4.460	0.078	0.587
Kinetic energy	KE/1046 erg	8.608	0.255	0.561
Buoyancy luminosity	${L}_{\mathrm{buoy}}/{10}^{45}$ erg s−1	4.576	−0.0342	−0.0492
Pressure	${\rm{\Delta }}\mathrm{ln}P$	2.032	0.046	0.228
Number of zonesa	${\rm{\Delta }}i$	236	8	23

Note.

^aThe total number of zones in the radial direction was 400 in Meakin & Arnett (2007b) (this table, medium resolution), 800 in Viallet et al. (2013) (high resolution), and 1536 in Figure 4. The basic features appear even at lower resolutions.

Download table as:  ASCIITypeset image

Table 2 shows the depth of each region in pressure scale heights (${\rm{\Delta }}\mathrm{ln}P$). The depth of the boundary zones is not a universal constant in $\mathrm{ln}P$, but varies by a factor of 5 between top and bottom. The last line gives the number of zones in each region for "medium" resolution (Meakin & Arnett 2007b); the lower boundary region is most demanding, having a steep transition from convective to stable stratification.

Little of the kinetic energy is lost in the boundary regions, so L_buoy provides a good first estimate of the rate of generation of TKE. These regions contain 17% of the mass in the "convection zone"; most of this comes from the upper layer, which has less extreme stratification.

Figure 3 shows the buoyancy flux versus radius, averaged over 100 s, for the oxygen-burning shell simulation (OB); more detail may be found in (Meakin & Arnett 2007b; Arnett et al. 2009; Viallet et al. 2013). The buoyancy flux, $-{\boldsymbol{u}}\cdot {\boldsymbol{g}}{\rho }^{\prime }/{\rho }_{0},$ is the rate of work done by gravity (Zahn 1991). It is the rate of flow of buoyancy, $-{\boldsymbol{g}}{\rho }^{\prime }/{\rho }_{0},$ and has units of energy per unit mass per unit time (e.g., erg g⁻¹ s⁻¹). Over most of the convective region it is proportional to the enthalpy flux (Arnett et al. 2009).

Figure 3 shows that the convective zone simulation is naturally split into three regions, separated by two boundaries. The regions above and below are stable. The middle region is relatively uninfluenced by the boundaries; it is characterized by positive fluxes of buoyancy and of enthalpy, that is, a positive "superadiabatic gradient" ${\rm{\Delta }}{\rm{\nabla }}$. It is convectively unstable according to both the Schwarzschild and the Ledoux criteria. With an appropriate²² choice of mixing length, this middle region can be reasonably well approximated by MLT.

MLT works poorly for the bottom and top boundary layers, which have negative values of ${\rm{\Delta }}{\rm{\nabla }}$. While the central region is defined by positive buoyancy, and positive enthalpy flux, outside the convective zone these quantities are zero, and in the boundaries they are negative. In MLT this is impossible because it would imply that the velocity in Equation (5) is imaginary, but in Equation (7) merely implies buoyancy braking, hence the labels "braking" in Figure 3.

Zahn (1991) has summarized²³ the issue of negative buoyancy and convective flux in connection with penetrative convection. Schmitt et al. (1984) have discussed the overshoot at the bottom of the solar convection zone in the context of convective plumes and magnetic dynamos, and Spiegel & Zahn (1992) have discussed this in the context of solar rotation and the tachocline. In stellar evolution theory (i.e., MLT) the existence of these braking regions is obscured by use of the Schwarzschild (or Ledoux) linear stability criterion. These braking layers are related to issues of overshoot and penetrative convection (Veronis 1963; Massaguer et al. 1984; Hurlburt et al. 1986). The braking layers are not a part of MLT but, as we shall see (Section 3.8), arise naturally from Equation (7).

Figure 4 shows the inner braking zone (the region of negative buoyancy work) at $r\sim (0.433$ to $0.445\times {10}^{9}$ cm). The "hi-res" case of Viallet et al. (2013) ($768\times {512}^{2}$ zones) and a still higher-resolution case of S. Campbell et al. (2015, in preparation) (1536 × 1024² zones) are shown. In comparison with Figure 3, the negative "spike" is now well-resolved. A detailed analysis of these simulations will appear elsewhere. The degree of numerical convergence is promising, and we conclude that such braking zones are a robust feature of well-resolved simulations of neutrino-cooled stellar convection.

Zoom In Zoom Out Reset image size

The radial velocity becomes small in the braking region, while the transverse velocity extends deeper before it also becomes small. The convective motion turns, and a small (mostly g-mode) wave velocity remains. The composition gradient is steeper than would be predicted by algorithmic diffusion (Equation (23)), and begins at the bottom of the braking region. The boundary composition profiles are smooth and self-similar when time-averaged. This suggests that the turbulent spectrum has a consistent net effect on the composition profiles and on the mixing, and therefore this interface should be amenable to approximation over time-steps in 1D evolutionary calculations.

For oxygen burning, the composition gradient in the boundary layer is not well-represented by conventional turbulent diffusion theory which requires a span of many "turbulence mean-free-paths" per density scale height (Amsden & Harlow 1968) for validity.²⁴ In MLT, the span is a fraction of a scale height (see ${\rm{\Delta }}\mathrm{ln}P$ in Table 2) for oxygen burning. The small length scales are accompanied by small time scales for change, so that a steady state model may be appropriate.

Fluid motion in a star may be separated into two fundamentally different flows (Landau & Lifshitz 1959): solenoidal flow (divergence free: ${\boldsymbol{\nabla }}\cdot \rho {\boldsymbol{u}}=0$) and potential flow (curl free: ${\boldsymbol{\nabla }}\times \rho {\boldsymbol{u}}=0$), which together represent the Helmholtz decomposition of an arbitrary vector field. Potential flow is associated with wave motion and solenoidal flow (vorticity) is a feature of turbulence. A striking separation in the nature of the flow is visible at boundaries between these types of flow; see the discussion of boundary layers in Prandtl & Tietjens (1934), Landau & Lifshitz (1959), and Figure 19 in Viallet et al. (2013). This separation in types of flow is closely related to wave generation and propagation (Press 1981; Press & Rybicki 1981; Goldreich & Kumar 1994; Goldreich et al. 1994).

The structure and nature of these boundary layers is important for estimation of the rate at which turbulent flow moves into or from non-turbulent regions—the growth and recession of convective zones. Meakin & Arnett (2007b) had about 8 zones across the lower boundary layer for "medium" resolution; see also Herwig et al. (2014). Viallet et al. (2013) had double the resolution across the convective zone (twice as many radial points), but the boundary layer became physically narrower. Recent simulations at still higher resolution (see Figure 4 and S. Campbell et al. 2015, in preparation) show that the lower boundary layer has about 20 zones and the same physical depth. The computed entrainment rate may be affected by numerical viscosity, so that lower resolution simulations will give overestimates.

The "medium" resolution of Meakin & Arnett (2007b) was sufficient to give numerical viscosity (Reynolds number) similar to that of laboratory experiments on entrainment, but not of stars. Coarse resolution in those simulations may have been a partial cause of the difficulties found by Staritsin (2013) in an attempt to apply the entrainment rates of Meakin & Arnett (2007b) for oxygen burning directly to main sequence stars. The real entrainment rates for stars should be smaller. Another issue is that oxygen burning and hydrogen burning have very different Péclet numbers (Viallet et al. 2015), which can affect the entrainment rate (see below).

Here we construct a simple but dynamically consistent picture of a convective boundary. This is illustrated in Figure 5, which shows the driving, turning, shear and stable regions. At its most elemental level, the velocity vector must turn at boundaries; that is, flow must turn back to stay inside the convective region. We do not assume that "blobs" disappear (like MLT). Most of the momentum is contained in the largest scales, so we focus on the average dynamics at these scales, and the simplest flow patterns.

Zoom In Zoom Out Reset image size

The magnitude of the acceleration required to turn the flow is just the centrifugal value ${u}^{2}/b$ where b is the radius of the turning region and u the relevant velocity. Using Equation (7) in the steady state limit, and taking $b\sim {\rm{\Delta }}r\ll {\ell }$, the radial component of the acceleration equation becomes

$\begin{eqnarray}&&{u}_{r}\partial {u}_{r}/\partial r\sim {\rm{\Delta }}(\displaystyle \frac{1}{2}{u}^{2})/{\rm{\Delta }}r\sim {\mathcal{B}},\end{eqnarray} \tag{ 33 }$

where ${\mathcal{B}}$ is the acceleration due to buoyancy and pressure fluctuations (Equation (6), and Section A.2). So far we have considered the top of a convective zone; the bottom of a convection zone behaves similarly if care is taken with signs.

Simulations (Meakin & Arnett 2007b; Viallet et al. 2013) show a consistent pattern in velocity and composition structure in the boundary layers. Moving toward the boundary from the interior of the convection zone, we find (1) the radial velocity u_r decreases, (2) the pressure fluctuations P' increase, and (3) the transverse velocity u_h increases to a maximum and then decreases, joining on to a finite and small rms velocity due to wave motion. The transition to small rms velocity occurs at about the same point that the composition changes from being well-mixed to supporting a radial composition gradient. This pattern holds for both top and bottom boundaries.

The dynamical equations we use are derived in Appendix A. We use Section A.2, the same quasi-steady state and thin shell ($b\ll {\ell }$) approximations, and choose an inertial frame in which a hydrostatic background is assumed. Near the boundary, the radial component of the acceleration is essentially just

$\begin{eqnarray}&&{\mathcal{B}}=-\left(\displaystyle \frac{{\rho }^{\prime }}{{\rho }_{0}+{\rho }^{\prime }}\right)g-\displaystyle \frac{1}{{\rho }_{0}+{\rho }^{\prime }}\partial P^{\prime} /\partial r.\end{eqnarray} \tag{ 34 }$

The buoyancy force (the first term on the rhs) is parallel to the gravity vector ${\boldsymbol{g}}$, which is radial, and provides no transverse acceleration. Baryon conservation implies that this reduction in the radial velocity alone will give an increase in density (matter accumulates), which gives an increase in the pressure fluctuation P' as the boundary is approached. The two transverse components of velocity satisfy

$\begin{eqnarray}&&{u}_{h}\partial {u}_{h}/\partial h=-\displaystyle \frac{1}{{\rho }_{0}+{\rho }^{\prime }}\partial P^{\prime} /\partial h.\end{eqnarray} \tag{ 35 }$

The transverse motion requires a transverse acceleration which is provided by a pressure excess (see also Stein & Nordlund 1989) at the point of contact of the plume with the boundary (note the similarity to the RTI, Section 2.5; and Schmitt et al. 1984).

This same pressure excess also implies a radial acceleration of the boundary, making the boundary undulate (Meakin & Arnett 2006, 2007b). In addition to the horizontal force from the pressure excess, the buoyancy force is negative, so the net effect on the flow is to complete the turn. The turning region has a width $b={r}_{2}-{r}_{1}$; this material is well-mixed because it moves back into the convective region after it turns. Thus the region ${r}_{2}-{r}_{1}$ might be termed the "over-shoot" region, and we are discussing the dynamics of "overshoot."

Figure 4 shows our highest resolution simulation of the most demanding boundary; does this simple model of boundary dynamics work for it? The orientation is reversed for the bottom boundary, so ${r}_{\mathrm{mix}}\lt {r}_{2}\lt {r}_{1}$ in this case. The steep drop in buoyancy work at $r\sim 0.433\times {10}^{9}\ \mathrm{cm}$ corresponds to r_mix and the "shear" region in Figure 5, which can maintain a composition gradient because the velocity is due to wave motion. At the radius r₂, at which the radial component of the velocity is ${u}_{r}\sim 0$, the flow is transverse to the radial coordinate (${u}_{h}\ne 0$), so there is a shear layer at this surface which will be unstable to the Kelvin–Helmholtz (KH) instability (Chandrasekhar 1961). The partial mixing layer extends to radius r_mix (at which ${u}_{h}\sim 0$) and contains this KH layer. The linear condition for ability of a layer to resist shear (stability against mixing) is the "gradient" Richardson number, $\mathrm{Ri}\gt 1/4$. The Brunt-Väisälä frequency $N\sim 3\ {{\rm{s}}}^{-1}$ is evaluated in the stable region, near the boundary, and may be sensitive to resolution. The shear velocity is ${u}_{h}\leqslant 0.8\times {10}^{7}\;\mathrm{cm}\;{{\rm{s}}}^{-1}$, and from this crude estimate ${r}_{\mathrm{mix}}-{r}_{2}\sim {u}_{h}/2N\sim {10}^{6}\ \mathrm{cm}$. This small length is consistent with the steep "cliff" in Figure 4.

Both terms in ${\mathcal{B}}$ (Equation (34)) act to turn the flow, and are comparable in magnitude. A crude but interesting estimate follows if we take ${\mathcal{B}}\sim g{\beta }_{T}{H}_{P}{\rm{\Delta }}{\rm{\nabla }}$, where the ${\rm{\Delta }}{\rm{\nabla }}$ is an average value over ${r}_{2}-{r}_{1}$. The turning radius in units of local pressure scale height is then

$\begin{eqnarray}&&b{/H}_{P}\sim {\rm{\Delta }}r{/H}_{P}\sim {\rm{\Delta }}(\displaystyle \frac{1}{2}{u}^{2})/g{\beta }_{T}{H}_{P}{\rm{\Delta }}{\rm{\nabla }},\end{eqnarray} \tag{ 36 }$

which is related to the inverse of a Richardson number; compare to Equations (29) and (31). Both ${\rm{\Delta }}(\frac{1}{2}{u}^{2})$ and ${\rm{\Delta }}{\rm{\nabla }}$ are negative here, giving a positive ratio. The use of Equation (7) automatically leads to an approximate Richardson number criterion for the edge of the convective region, without the need of an additional imposed boundary condition beyond the requirement that u² becomes small (see Section 3.5).

The minimum in buoyancy work at $r\sim 0.437\times {10}^{9}\ \mathrm{cm}$ corresponds to r₂, the edge of the braking region and the "turn" in Figure 5. At $r\sim 0.443\times {10}^{9}\ \mathrm{cm}$ the buoyancy work becomes positive, so that this corresponds to r₁ and the beginning of the "driving" region, at which ${\rm{\Delta }}{\rm{\nabla }}$ changes sign. Contrary to MLT, the radius r₁, at which the Schwarzschild criterion is zero, is not at the boundary of zero convective motion.

How does this braking region develop a negative buoyancy? Suppose the region r₂ to r₁ is well mixed, to uniform composition and entropy. There is no braking, so convective flow is unabated to the composition gradient beginning at r₂. Vigorous entrainment erodes the boundary, causing a thin layer of partially mixed matter, which contains the heavier nuclei from below the oxygen burning shell. This makes the buoyancy more negative, establishing a braking layer and reducing the rate of entrainment. The braking layer grows until the entrainment rate balances the rate of mixing into the edge of the convection zone. If the braking layer is too large, such mixing will reduce it; there is negative feedback. The braking layer is thinner than the convective zone, so the time scale is shorter than the turnover time (Section 2.1), and a quasi-steady state can be set up. This simplistic analysis (which ignores fluctuations) indicates some of the dynamics involved with the braking layers and composition boundaries. Further analysis with the new higher resolution simulations (S. Campbell et al. 2015, in preparation, C. Meakin et al. 2015, in preparation) is in progress.

This limiting case ("elastic collision") is a reasonable approximation for the time averaged behavior of the oxygen burning shell (Meakin & Arnett 2007b), in which radiative diffusion (and electron heat conduction) are slow; here ${\tau }_{\mathrm{turn}}\sim 0.6\;{\rm{s}}$, while the radiative diffusion time is ${\tau }_{\mathrm{diff}}\sim 3\times {10}^{7}\;{\rm{s}}$. A measure of the heat lost during the turn is a small number ($\sim 2\times {10}^{-8}$) for oxygen burning, and is roughly the inverse of the Péclet number. Even within the narrow braking layer, there is little heat flow by radiative diffusion during oxygen burning.

This discussion underestimates mixing because it ignores turbulent fluctuations (Section 2.6); larger fluctuations do more mixing than average, and mixing is irreversible. Turbulent kinetic energies fluctuate by factors ∼2, so the mixing estimates should be increased accordingly. Flow velocities do not go to zero at the convective boundaries, but become small and oscillatory (Press 1981; Press & Rybicki 1981; Goldreich & Kumar 1994; Goldreich et al. 1994). As convective plumes hit the boundary, and rebound, the boundary moves in response; how elastic this is depends upon heat flow (the Péclet number).

This "adiabatic" limit breaks down as the turnover time ${\tau }_{\mathrm{turn}}\sim b/u$ approaches the radiative diffusion time for the turn ${\tau }_{\mathrm{diff}}\sim {b}^{2}/\lambda c$. For larger radiation mean-free-paths, the Péclet number decreases. No sharp temperature gradients can persist. This gives an "inelastic collision" of the flow with the boundary. This is the case for stars in photon-cooled stages of evolution; even with relatively large Péclet numbers for the whole convective region, the narrow boundary layers may still have significant energy flow by radiative diffusion. The previous discussion of the effect of excess pressure P' still holds, but because of thermal diffusion P' becomes increasingly dominated by density excess $\rho ^{\prime}$ rather than the temperature excess T'.

The red giant model of Viallet et al. (2011) provides an example of a boundary layer (the bottom) in which there is significant radiative diffusion; Viallet et al. (2013) analyze this in detail (their Section 4.6). As the boundary is approached from above, the down-flows are accelerated by pressure dilatation. These down-flows have an entropy deficit, so that they are heated by radiative diffusion from the surrounding material. In the braking region, compression causes a "hot spot" to develop. The flow is turned to a non-radial direction, and is now cooled by radiative diffusion (see Figure 7 in Viallet et al. 2013).

Such behavior differs from that obtained by present stellar evolution algorithms. The turning of the down-flow forces the mixed region to extend beyond that implied by the Schwarzschild criterion, and heating/cooling by radiative diffusion modifies the structure. While modest, such differences can be important for detailed models. In compensation for such changes, a standard solar model requires less opacity to have the same convection zone depth; this implies a lower metallicity. These changes in the solar model provide a means to reduce the disagreement with helioseismology (Christensen-Dalsgaard et al. 2011; Zhang et al. 2012). Deng & Xiong (2008) gave a justification for compositional smoothing, as did simulations (Meakin & Arnett 2007b; Viallet et al. 2013). The thermal characteristics needed (Christensen-Dalsgaard et al. 2011) follow from the analysis given above, which was not designed for the solar problem, and involved no solar or stellar calibration. A more physically correct convective boundary condition tends to improve agreement with abundances inferred from 3D stellar atmospheres (Asplund 2005) and the standard solar model.

If heat flow processes are included, the "inelastic collision" with the boundary allows the loss of heat so that the entropy decreases for the downward flow, enhancing the downward acceleration. This effect tends to drive motion in convective envelopes. Heating at the bottom also tends to drive convective flow. However, cooling at the bottom (as with URCA-shells, Arnett 1996) or heating at the top (downwardly entrained, burning fuel) both tend to halt the flow. Such halting processes can cause convective zones to split (Mocák et al. 2011b).

There may be observational evidence supporting this description of boundaries of convection which are deep in stellar interiors. Detection of g-mode pulsations in sdB stars allows an asteroseismic estimation of the size of the He-burning cores, which are significantly larger than predicted by the Schwarzschild criterion and standard stellar evolution theory (see Schindler et al. 2015 for discussion and references). Similar issues apparently are general for core helium burning stars observed by Kepler (Constantino et al. 2014; Mosser et al. 2014).

Finally, the origin (r = 0), in a 1D stellar evolutionary code using MLT, is a boundary as well. The use of Equation (5) with adequate zoning implies that the convective velocity becomes very small due to symmetry (derivatives go to zero at the origin). This is a false braking layer caused by MLT being a local theory. Use of Equation (6) allows flow through the origin provided that a counter flow gives conservation of linear momentum (e.g., a toroidal roll). At the origin in a turbulent convective core, this projects onto 1D as a finite rms velocity, with a zero radial gradient. MLT has problems with velocity at r = 0.

The enormous simplification, from 3D turbulent simulations requiring terabytes of storage down to a single additional ordinary differential equation (e.g., Equation (6)), means that much is missing. For some applications the missing items may be important. One might use the RANS equations directly in a stellar evolutionary code, with 3D simulations to guide closure (Mocák et al. 2014). We have presented a step toward that goal. Alternatively, one might add to the simple 321D as needed, using new models guided by RANS results. Probably both paths should be followed, given the complexity of the problem.

4.1.1. 321D algorithms

4.1.2. Further Simulations

Because of the fundamental importance of convection in stellar evolution theory, a replacement for MLT will have implications for many areas throughout astronomy and astrophysics. A few of the most striking are as follows.

4.2.1. Helioseismology

4.2.2. Asteroseismology

4.2.3. Convective Boundaries, Nucleosynthesis Yields and Pre-supernovae

4.2.4. Core Collapse

The vector form of the continuity equation (Landau & Lifshitz 1959) is

$\begin{eqnarray}&&\partial \rho /\partial t+{\boldsymbol{\nabla }}\cdot \rho {\boldsymbol{u}}=0,\end{eqnarray} \tag{ 37 }$

where ρ is the mass density and ${\boldsymbol{u}}$ is the fluid velocity. In the incompressible limit, for a steady flow, the net flux of mass into a region equals the mass flux out. In thin boundary layer, perpendicular to the radial direction r, the average velocities must satisfy

$\begin{eqnarray}&&\partial {u}_{r}/\partial r=-2\partial {u}_{h}/\partial h,\end{eqnarray} \tag{ 38 }$

where h is either of the symmetric transverse coordinates (i.e., locally cartesian), to avoid changing the density (as seen in the Eulerian frame).

Viallet et al. (2013) show (their Equation (28)), that for fluctuations against a steady background,

$\begin{eqnarray}&&{\boldsymbol{\nabla }}\cdot {\boldsymbol{u}}^{\prime} =\displaystyle \frac{{u}_{r}^{\prime }}{{H}_{\rho }},\end{eqnarray} \tag{ 39 }$

where ${H}_{\rho }$ is the density scale height, and ${u}_{r}^{\prime }$ is the radial component of the velocity fluctuation. This approaches zero (the incompressible limit) for shallow, subsonic convection (large density scale height and small radial velocity mach number, ${u}_{r}^{\prime }\ll s$, where s is the sound speed). This velocity "dilatation" is due to the vertical motion in the background stratification and becomes an important component in convective driving in deep convection zones (Viallet et al. 2013). Notice that rising plumes (${u}_{r}\gt 0$) expand and falling plumes contract (Stein & Nordlund 1989; Meakin & Arnett 2010).

The vector acceleration equation (Equation (6)) is

$\begin{eqnarray}&&\partial {\boldsymbol{u}}/\partial t+({\boldsymbol{u}}\;{\boldsymbol{\cdot }}\;{\boldsymbol{\nabla }}){\boldsymbol{u}}={\mathcal{B}}-{\boldsymbol{u}}/\tau \end{eqnarray} \tag{ 40 }$

where ${\boldsymbol{u}}$ is the velocity, $\tau =| u| /{{\ell }}_{{\rm{d}}}$ with ℓ_d is the Kolmogorov damping length, and the variable ${\mathcal{B}}$ is defined as in Section 2.3. If

$\begin{eqnarray}&&{\mathcal{B}}=-\displaystyle \frac{1}{\rho }{\boldsymbol{\nabla }}P-{\boldsymbol{g}},\end{eqnarray} \tag{ 41 }$

where P is pressure and ${\boldsymbol{g}}$ is gravitational acceleration, then Equation (40) is a Navier–Stokes description of the largest scales of turbulence, with a simplified damping term which is consistent with Kolmogorov (1962). Note that the usual formulation of hydrostatic equilibrium in stellar evolution theory is some variant of the condition ${\mathcal{B}}=0$. Projecting Equation (40) onto the radial coordinate, we have

$\begin{eqnarray}&&\partial {u}_{r}/\partial t+{u}_{r}\partial {u}_{r}/\partial r=-\displaystyle \frac{1}{\rho }\partial P/\partial r-g-{u}_{r}/\tau .\end{eqnarray} \tag{ 42 }$

The full equations in spherical coordinates are shown in Section 15, Landau & Lifshitz (1959) (see also Mihalas & Mihalas 1984 for a detailed discussion), with the bare viscosity terms rather than Komogorov's expression for integration of the turbulent cascade. In tensor form the momentum equation is

$\begin{eqnarray}&&\partial {u}_{i}/\partial t+{u}_{k}\partial {u}_{i}/\partial {x}_{k}=-\displaystyle \frac{1}{\rho }\partial P/\partial {x}_{i}-{g}_{i}+\displaystyle \frac{1}{\rho }\displaystyle \frac{\partial }{\partial {x}_{k}}\\ &&\qquad \times \left[\eta \left(\displaystyle \frac{\partial {u}_{i}}{\partial {x}_{k}}+\displaystyle \frac{\partial {u}_{k}}{\partial {x}_{i}}-\displaystyle \frac{2}{3}{\delta }_{{ik}}\displaystyle \frac{\partial {u}_{l}}{\partial {x}_{l}}\right)\right]+\displaystyle \frac{\partial }{\partial {x}_{i}}\left(\zeta \displaystyle \frac{\partial {u}_{l}}{\partial {x}_{l}}\right).\end{eqnarray} \tag{ 43 }$

Kolmogorov's four-fifths law (Frisch 1995) states an amazing simplification, that integration over the turbulent cascade reduces the last term in Equation (43) to $-{\boldsymbol{u}}/\tau$ (Equation (40)) on average, ignoring boundary effects (see Section 3).

To illustrate how turning happens at boundaries, it is sufficient to consider the simpler case of flows with θ and ϕ length scales small compared to r, so the transverse dimensions are quasi-cartesian (the inertial terms in $1/r$ are neglected; for convective cores, the more cumbersome full equations are needed because r cannot be large near the origin). Then the two transverse components are symmetric in this approximation and satisfy

$\begin{eqnarray}&&\partial {u}_{h}/\partial t+{u}_{h}\partial {u}_{h}/\partial h=-\displaystyle \frac{1}{\rho }\partial P/\partial h-{u}_{h}/\tau ,\end{eqnarray} \tag{ 44 }$

where dh is ${rd}\theta$ or $r\mathrm{sin}\theta d\phi$. We consider finite fluctuations about a static background, so that we substitute $\rho ={\rho }_{0}+\rho ^{\prime}$ and $P={P}_{0}+P^{\prime}$. We ignore variations in g (the Cowling approximation, Cox 1980). Using $-\partial {P}_{0}/\partial r={\rho }_{0}g$, the radial equation becomes

$\begin{eqnarray}\partial {u}_{r}/\partial t+{u}_{r}\partial {u}_{r}/\partial r & = & -\left(\displaystyle \frac{{\rho }^{\prime }}{{\rho }_{0}+{\rho }^{\prime }}\right)g\\ & & -\displaystyle \frac{1}{{\rho }_{0}+{\rho }^{\prime }}\partial P^{\prime} /\partial r-{u}_{r}/\tau .\end{eqnarray} \tag{ 45 }$

Convection is often described using only the buoyancy term; the pressure fluctuations are taken to be small, of order the mach number squared. However, near boundaries the pressure fluctuations provide the tangential acceleration which is necessary to turn the flow, and should not be neglected (see Nordlund 1985). The buoyancy term acts through the density fluctuation $\rho ^{\prime}$, and only in the direction parallel to the gravity vector. The transverse equation is

$\begin{eqnarray}&&\partial {u}_{h}/\partial t+{u}_{h}\partial {u}_{h}/\partial h=-\displaystyle \frac{1}{{\rho }_{0}+{\rho }^{\prime }}\partial P^{\prime} /\partial h-{u}_{h}/\tau .\end{eqnarray} \tag{ 46 }$

Note that the radial and transverse equations are coupled primarily by the pressure fluctuation term P', but also by $u/\tau$, because $\tau ={{\ell }}_{{\rm{d}}}/| u|$ where $| u{| }^{2}={u}^{2}={u}_{r}^{2}+2{u}_{h}^{2}$ (turbulence damps regardless of orientation of the large scale flow). The fluctuating pressure near convective boundaries insures the generation of waves.

Following Landau & Lifshitz (1959), Section 6, the equation of energy conservation is

$\begin{eqnarray}\displaystyle \frac{\partial }{\partial t}\left(\displaystyle \frac{1}{2}\rho {u}^{2}+\rho E+\rho \phi \right) & = & -{\boldsymbol{\nabla }}\cdot \left[\rho {\boldsymbol{u}}(\displaystyle \frac{1}{2}{u}^{2}+W+\phi )\right]\\ & & +T\displaystyle \frac{\partial \rho S}{\partial t},\end{eqnarray} \tag{ 47 }$

where ϕ is the gravitational potential and ${\boldsymbol{g}}=-{\rm{\nabla }}\phi$. If taken to both the steady state and adiabatic limits, this becomes the Bernoulli equation (Landau & Lifshitz 1959). The entropy change equation may be written as

$\begin{eqnarray}&&T\displaystyle \frac{\partial \rho S}{\partial t}=\rho {\epsilon }_{\mathrm{nuc}}+\rho {\epsilon }_{\mathrm{visc}}-{{\boldsymbol{F}}}_{\mathrm{rad}},\end{eqnarray} \tag{ 48 }$

where ${\epsilon }_{\mathrm{nuc}}$ is the net heating from nuclear and neutrino reactions, $\rho {\epsilon }_{\mathrm{visc}}$ is the Navier–Stokes viscous heating term as modified by Kolmogorov's four-fifth's law (see Equations (1), (40) and (43)), and F_rad is the energy flux due to radiative diffusion. The viscous term is missing from MLT and the Euler equation. Most of the TKE resides in the largest (integral) scale, while turbulent heating occurs at the small (Kolmogorov) scale. Then ${\epsilon }_{\mathrm{turb}}={\boldsymbol{u}}\cdot {\boldsymbol{u}}| u| /{{\ell }}_{{\rm{d}}}$ is the Kolmogorov heating from the turbulent cascade, and $T\partial \rho S/\partial t$, $\rho {\epsilon }_{\mathrm{nuc}}$ and F_rad are now the appropriate RANS averages (Viallet et al. 2013). One requirement for Bernoulli's equation to be valid, as assumed in Pasetto et al. (2014) (see Section 2.8), is that the rhs of Equation (48) must be zero (Landau & Lifshitz 1959, Ch. I). This is found not to be generally true, either in the 3D simulations (Viallet et al. 2013; Mocák et al. 2014), or experimentally in turbulent flows (Tennekes & Lumley 1972; Davidson 2004). Heating is an essential feature of 3D turbulence, which converts large scale, ordered velocities to disordered ones.