TURBULENT CONVECTION IN STELLAR INTERIORS. II. THE VELOCITY FIELD

We start with the equations we used in our 3D simulations (Meakin & Arnett 2007b). We use Reynolds decomposition of relevant flow and thermodynamic variables, which separates the fluctuating components from the slowly varying components (Tennekes & Lumley 1972). Numerical simulations dissipate features with wavelengths at or below the grid scale, and we will identify this with dissipation in the turbulent cascade (Kolmogorov 1941, 1962).

In the process of approximation of partial differential equations by finite methods, there is inevitably a loss of information at scales smaller than the grid size. A single volume element in space is approximated as a homogeneous entity; this is equivalent to complete mixing at this scale, at each time step, of mass, momentum, and energy. The loss of information that occurs with this mixing corresponds to an increase in entropy (Shannon 1948); the mixing of momentum is equivalent to the action of viscosity (Landau & Lifshitz 1959). In 3D flow, turbulent energy will cascade from large scales to small, at a rate set by the largest scales. We use an implicit subgrid dissipation in our large eddy simulation (ILES), which is the most computationally efficient way to deal with turbulent systems with a large range of scales (Boris 2007; Woodward 2007); the largest scales, which set the rate of cascade and contain most of the energy, are resolved on our grid and explicitly calculated, while the subgrid scales are dissipated in a way consistent with the Kolmogorov cascade.

Sytine et al. (2000) demonstrated that the piecewise parabolic method (PPM) based on the Euler equation (which has no explicit viscosity) converges to the same limit as methods based on compressible viscous equations (which do have explicit viscosity), as the grid is refined to smaller zones and smaller effective viscosity (the relevant limit for astrophysics). Porter et al. (2000) show the compatibility of mildly compressible flow with the Kolmogorov (1941) spectrum; Kritsuk et al. (2007) have pushed this to highly compressible flows as well. To represent the subgrid dissipation, which is inherent in our simulations, we explicitly introduce a volumetric dissipation rate for KE, ε_K, in our theoretical analysis. However, the turbulent cascade is a property of the whole convective flow, so connection to turbulence theory must be made through integrals of ε_K over the convection zone. Note that this is different from defining ε_K as a function of local variables (e.g., Smagorinsky 1963; Chan & Sofia 1989; Hansen & Kawaler 1994; Asida & Arnett 2000); see below.

The KE equation for convective motion was given in Meakin & Arnett (2007b). Taking the scalar product of the velocity with the equation of motion, we decompose the convective velocity u, the density ρ, and the pressure p into mean and fluctuating components (e.g., p = p₀ + p'), so the time averages are $\overline{p} = p_{0}$ and $\overline{p^{\prime }} =0$. This choice of just u, p, and ρ for this Reynolds decomposition into average and fluctuating parts gives the simplest equation for KE; the velocity u is derived from buoyancy and pressure forces (ρ and p fluctuations).

Using the hydrostatic equilibrium condition, and performing averages, gives (see Equation (A12) in Meakin & Arnett 2007b),

$$\begin{eqnarray} &&\partial _t \langle \overline{\rho E_K}\rangle + \nabla \cdot \langle \overline{\rho E_K {\bf u_0}}\rangle \nonumber \\ &&\quad = - \nabla \cdot \langle \overline{{\bf F_p + F_K }}\rangle + \langle \overline{ p^{\prime } \nabla \cdot u^{\prime } }\rangle \nonumber \\ &&\quad\ \ \ \ + \langle \overline{ \rho ^{\prime } {\bf g \cdot u^{\prime }} }\rangle - \varepsilon _K. \end{eqnarray} \tag{ 1 }$$

We use 〈p〉 to denote an average over angles at constant radius. The time average is taken over durations greater than a transit time t_transit = ℓ_CZ/(u')_rms, where ℓ_CZ is the depth of the convective region, and (u')_rms is the rms velocity across this region. This smooths the fluctuations, gives a nonlocal character to the analysis, and implies a separation of timescales into short (t ≪ t_transit) and long (t ≫ t_transit). We consider the case in which we may integrate over these short timescales and explicitly calculate the evolution on the long timescales.

Here E_K is the KE per unit mass, ρ is the mass density, u₀ is the nonfluctuating part of the fluid velocity vector and u' its fluctuating part, so

$$\begin{equation} {\bf F}_p = p^{\prime } {\bf u^{\prime }} \end{equation} \tag{ 2 }$$

is the energy flux due to pressure perturbations carried by fluctuations (this pressure–velocity correlation flux reduces to the acoustic flux when considering sound waves (Landau & Lifshitz 1959, p. 251), and contains the energy flux due to internal gravity waves),

$$\begin{equation} {\bf F}_K = \rho E_K {\bf u^{\prime }} \end{equation} \tag{ 3 }$$

is the flux of KE carried by convective turbulent motion, g is the gravitational acceleration vector, and ε_K is the volumetric rate of dissipation implied by the turbulent cascade down to small scales.

The location of convective boundaries is a long-standing problem in stellar astrophysics. It has long been known that some sort of "convective overshooting" is necessary for models to match observations. This has been conceived as mixing of material beyond the convective boundary, often by "penetrative convection." This problem exists in part because the convective boundary has been inappropriately defined (Meakin & Arnett 2007b, especially Sections 4.1 and 7). The standard definition used in MLT defines the convective boundaries using the thermodynamic nablas, which essentially mark the onset and cessation of buoyant acceleration. One of the essential properties of convection is that turbulence fills the space available. A more appropriate criterion is one in which the stellar background is stiff enough to contain the turbulent convection, and therefore must account for the relative strength of the turbulent flow and the elastic stable layers.

We will define the convective zone as that region in which the stratification of the medium is unstable to turbulent mixing. This is evaluated with the "bulk Richardson number"

$$\begin{equation} Ri_B = \frac{\Delta b \ l}{u^2}, \end{equation} \tag{ 4 }$$

where Δb = ∫_ΔrN²dr is the buoyancy jump across a layer of thickness Δr in the radial direction, N² is the Brünt–Väisälä frequency, u is the rms velocity providing shear, and l is the scale length of the turbulence (essentially the size of the largest eddies). The precise definition of the thickness Δr is a topic that deserves further study; in several cases we have noticed that the rapid change in N² near a boundary tends to make Δb insensitive to the exact value of Δr chosen for integration to obtain Δb.

In the linear limit of plane parallel flow, a region with Ri_g ≲ 0.25 has enough KE to overcome the stable stratification (Dutton 1986). Here Ri_g is the "gradient" Richardson number, and is a locally defined quantity, used in a linear stability analysis. The "bulk" Richardson number Ri_B is an inherently nonlinear quantity based on integration over an extended region. Since Ri_g,crit > 0, this allows for N² > 0 as well. This formulation recognizes that layers that are thermodynamically stable (real roots to N²) can be hydrodynamically unstable when KE is input to the zone. The bulk Richardson criterion allows more mixing than predicted by the Ledoux criterion; the Schwarzschild criterion ignores compositional effects, and also predicts more mixing than Ledoux (see Kippenhahn & Weigert 1990; Hansen & Kawaler 1994).

Stable regions adjacent to the convective zone will become Richardson unstable periodically as shear builds up from adjacent turbulent motions and waves generated by convection, leading to entrainment of material at convective boundaries (Fernando 1991). The bulk Richardson number is not only an indicator of a convective boundary; it also determines the "rate" at which the boundary migrates through the Lagrangian mesh by mixing processes (Meakin & Arnett 2007b), and thus provides a "dynamic" boundary condition for the flow.

Our analysis is based primarily upon the numerical simulations of oxygen shell burning in a 23 M_☉ star, but is found to be broadly applicable to other examples of stellar convection (such as convective envelopes, which are driven by the superadiabatic layer at the top rather than the nuclear burning shell at the bottom). Our convective zone has a depth of two pressure scale heights (see Meakin & Arnett 2007b for details).

Why do we need to do averaging? Figure 1 shows the behavior of the largest terms in the KE equation (averaged over the convection zone), throughout the duration of the simulation. The dominant terms are the integrals of the buoyancy driving $\int W_b d{\cal V}= \int g \rho ^{\prime } u^{\prime } d{\cal V}$, which is positive, and of damping, $\int \varepsilon _K d{\cal V}$, which is negative. Both show recurrent bursts (ten in 800 s). The time derivative of KE in the convection zone (labeled dE_K/dt) is constant on average. The remaining bumps are due to missing terms (divergence of F_P and F_K).

Zoom In Zoom Out Reset image size

We saved complete data every 0.5 s, which was adequate for making movies, and constructing averages of motion. A key issue is distinguishing between turbulent velocities and wave velocities with such a coarse stride in time. We have found a way to split the velocity averages into turbulent and wave components for better accuracy in the analysis (see the Appendix for a detailed discussion of how this is done).

The level of complexity shown in Figure 1 is daunting, but the near balance of buoyant driving and turbulent damping provides a clue.

The system is simplified if we integrate over time; the resulting values over four transit times are shown in Table 1. We do not calculate the damping directly, but deduce it as the remainder left from the other five terms, so that the entries sum to zero. We are comfortable with this procedure because of the good conservation properties of the numerical simulations.

Table 1. Kinetic Energy Equation Terms Averaged Over Convection Zone and Over Four Transit Times (erg s⁻¹)

Term	Value	Term	Value
$\int \langle \overline{ W_{B} } \rangle d{\cal V}$	4.576(45)	${-}\!\int \langle \overline{\varepsilon _{K}} \rangle d{\cal V}$	−4.677(45)
${-}\!\int \langle \overline{ \nabla \cdot F_{K} } \rangle d{\cal V}$	2.584(44)	${-}\!\int \langle \overline{ \nabla \cdot F_{P} } \rangle d{\cal V}$	−9.922(43)
${-}\!\int \langle \overline{ dE_{K}/dt } \rangle d{\cal V}$	−5.790(43)	$\int \langle \overline{p^{\prime } \nabla u^{\prime }} \rangle d{\cal V}$	−1.516(41)

Download table as:  ASCIITypeset image

The dominance of buoyant driving and damping is now clear; the next largest term is the divergence of the KE flux, at less than 6% of the buoyancy term.

Suppose the globally averaged damping term has the Kolmogorov form, so that

$$\begin{equation} \int _{\rm CZ} \varepsilon _K d {\cal V} = M_{\rm CZ} {v_{\rm rms}}^{3}/ \ell _{D}, \end{equation} \tag{ 5 }$$

where ℓ_D is a constant "damping length," and v_rms is the rms velocity determined over the entire convection zone. For a turbulent spectrum, ℓ_D is taken to be the largest length scale (Landau & Lifshitz 1959). In what follows we will use the term Kolmogorov damping to mean damping due to a turbulent cascade to small scales, having this cubic dependence on rms velocity.

Unlike MLT, we relate the damping length to the largest scale in the flow; this may be thought of as a "coherent structure" or a "plume which traverses the depth of the convection zone." It is not a free parameter. We introduce the notation

$$\begin{equation} \alpha _{D}= \ell _{D}/\ell _{\rm CZ}, \end{equation} \tag{ 6 }$$

so that α_D would be constant if the size of the largest eddy scales simply with the depth of the convection zone. Kritsuk et al. (2007) found the Kolmogorov damping to be very close to constant during a statistical steady state in their 3D (2048³) simulations.

The mass contained in the convection zone is M_CZ = 1.84 × 10³³g and the total KE is $\int E_{K} d{\cal V} = 8.61 \times 10^{46} \rm \ ergs$ so that the rms velocity is v_rms = 9.66 × 10⁶cm/s. This gives ℓ_D = 3.56 × 10⁸ cm ≈ 0.85 ℓ_CZ. It is natural to compare the damping length to the depth of the convection zone. Not only is ℓ_CZ the largest length available for an eddy but also if measured in pressure scale height units (ℓ_CZ/H_P), it indicates the degree of thermodynamic anisotropy across the convective region.

Table 2 gives several relevant timescales in seconds. Although these times are short in human terms, they are much closer to thermal relaxation times than are the corresponding numbers for deep simulations of the solar convection zone (our simulations are much more relaxed in this sense).

Table 2. Timescales

Timescale	Definition	Value (s)
Buoyant rise	$\int W_b d{\cal V}/KE$	14.1
Velocity damping	ℓD/vrms	27.6
KE damping	ℓD/(2vrms)	13.8
Transit time	ℓCZ/vrms	51.4
Turnover time	2ℓCZ/vrms	102.8

Download table as:  ASCIITypeset image

We may write the global damping as $M_{\rm CZ} v_{\rm rms}^{3}\big/ \ell _{D}= {1 \over 2}M_{\rm CZ} v_{\rm rms}^{2}\big/ (\ell _{D}\big/2 v_{\rm rms})$, and see that the damping time for KE in the convective zone is half the time to transit a damping length (which is approximately the depth of the convection zone). The turnover time for the convection zone is 2τ, where τ = ℓ_CZ/v_rms is the transit time, to be precise. The rise time for KE and the corresponding damping time are similar (14 s), and much shorter than a turnover time.

The term "convective efficiency" in the stellar context is usually taken to mean that the timescale for convective energy transport is much less than the timescale for radiative energy transport (Hansen & Kawaler 1994, p. 187); this insures that the actual temperature gradient is only slightly in excess of the adiabatic one. Convection can also be thought of as a thermodynamic cycle, taking a time 2ℓ_CZ/v_rms. As we saw above, the dissipation timescale for KE is about one-seventh of this (0.134), so that in this sense, convection is not thermodynamically efficient at all, but requires continual work to keep it running. The two uses of "efficiency" causes confusion, and downplays the fact that stellar convection is highly dissipative, even for slightly superadiabatic temperature gradients.

Here we present an initial, qualitative discussion of the structure of the flow, which we are investigating qualitatively for subsequent publication. We extend the idea of Cattaneo et al. (1991) that the convective velocity field has two components, a more ordered global flow and a chaotic turbulent flow. We focus on the source and sink for convective KE. Gravitational acceleration breaks the symmetry of space; we choose our z-axis parallel to this acceleration. Buoyant acceleration starts an anisotropic flow in the z-direction. The flow is unstable, and begins to break up into smaller scales, and becomes more isotropic.

Suppose that the flow occurs in narrow, vigorous down plumes and slow wide upflows for convective zones with significant anisotropy. The KE in such a downflow would dominate that in the upflow; the KE flux has two more powers of velocity than the mass flux, so that the fastest flows dominate the KE budget (Meakin & Arnett 2007b). We separate the KE into two components, u²_a/2, which corresponds to the largest eddies and u²_b/2, which represents all of the turbulent cascade. The turbulent instability begins to make the flow more isotropic, so that

$$\begin{equation} u_{z}^{2} \approx u_{a}^{2} + u_{b}^{2}, \end{equation} \tag{ 7 }$$

in the z-direction (see Batchelor 1960). If this KE is shared with the two perpendicular directions, then

$$\begin{equation} u_{x}^{2} + u_{y}^{2} \approx u_{z}^{2}. \end{equation} \tag{ 8 }$$

Further, it seems that the cascade component in the z-direction is similar to the components in the x- or the y-directions, u²_x ≈ u²_y ≈ u²_b. While the downflow matter is balanced by a return flow in x and y, the return flow is slower and broader (by our assumptions), with lower KE (which we neglect here). Thus

$$\begin{equation} u_{b}^{2} \approx u_{a}^{2}. \end{equation} \tag{ 9 }$$

One quarter of the KE is in the largest scale component u²_a/2. The ratio of vertical to horizontal rms velocities is therefore $u_{x}/u_{z} \approx 1/\sqrt{2}$.

How does this simple model compare with simulations? In their Figure 6, Meakin & Arnett (2007b) give the rms values of velocity in the vertical (v_r) and the horizontal (v_θ and v_ϕ) directions. We identify locally the r direction with the z-axis, and theta and phi with the x- and y-axes. Then from their Figure 6, we find

$$\begin{equation} 0.5 \le \frac{u_{x}}{u_{z}} \approx \frac{u_{y}}{u_{z}} \le 0.8, \end{equation} \tag{ 10 }$$

over the convection zone, away from the boundaries. The measured value of the vertical velocity contains both large-scale and cascade components. This is consistent with our estimate of $u_{x}/u_{z} \approx u_{y}/u_{z} \approx 1/\sqrt{2} \approx 0.707$. The large-scale component of KE is comparable to the cascade one, consistent with Cattaneo et al. (1991), and our simulations.

The simulations have led us to a simple, two-component model for the origin and destruction of convective KE. Anisotropy is due to buoyant acceleration in the largest eddy, which represents one component. The second is isotropic turbulence, which gives dissipation of KE after the cascade to small scales. It will be interesting to test and refine this model against a wide variety of simulations. The power spectra of the flow will be a useful tool to explore this ansatz.

In what follows we will define a turbulent rms velocity u_t from the measured horizontal velocity components, $u_{t}^{2} = {3 \over 2}\big (u_{x}^{2} +u_{y}^{2}\big) = {3 \over 4} u_{\rm rms}^{2}$. We note for future reference that $u_{\rm rms}^{3}\big/u_{t}^{3} = \big({4 \over 3}\big)^{{3 \over 2}} = 1.540$. Roughly speaking, we identify u_t with the turbulent cascade. The anisotropic component u_a is identified predominately with the largest turbulent scale, and may be sensitive to the structure of the convection zone and the location and strength of the buoyant acceleration.

We now relax one level of averaging, and consider the radial variation of all the terms in the KE equation. Figures 2–4 show six terms in six panels, each compared with the inferred damping (the solid line present in each panel). This allows us to determine what each of the five terms contribute to the inferred damping, and also shows an isotropic Kolmogorov estimate as a simple approximation (Figure 2, bottom panel). The solid line has a sharp spike at the bottom of the convection zone. This feature is due to the dominance of g-mode character in the motion near the boundary region. The Kolmogorov cascade is appropriate away from the boundaries.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In Figure 2, the top panel shows the averaged buoyancy $\int \overline{\langle W_B d{\cal V}\rangle }$, and the bottom panel shows the isotropic Kolmogorov approximation, both as dotted curves. In planar geometry (a fair approximation), we can relate this quantity to the "buoyancy flux" which is widely used in experimental fluid mechanics (Turner 1973). The buoyancy flux is

$$\begin{equation} {\bf q} = \frac{\relax \overline{\langle {\bf g \cdot u^{\prime } }\rho ^{\prime } \rangle }}{ \rho _{0}}, \end{equation} \tag{ 11 }$$

so that $\int \overline{\langle W_B d{\cal V}\rangle }\rightarrow \int \overline{\langle q \rangle }4 \pi r^{2}dr$. The enthalpy (convective heat) flux is

$$\begin{equation} {\bf F_e} = \relax \overline{\langle \rho _0 {\bf u^{\prime }} C_p T^{\prime } \rangle }, \end{equation} \tag{ 12 }$$

where C_p is the specific heat at constant pressure. For low Mach-number flows like ours⁷ the pressure fluctuations are small. The temperature and density perturbations at constant pressure are proportional,

$$\begin{equation} \frac{\rho ^{\prime }}{\rho _0} = \frac{\beta _{T}\ T^{\prime }}{T_0}, \end{equation} \tag{ 13 }$$

where β_T = −∂ln ρ/∂ln T_P, taken at constant pressure (and composition)⁸. We find, using hydrostatic equilibrium of the unperturbed star H_pg = P₀/ρ₀,

$$\begin{equation} {\bf F_e} = \frac{\rho _0 H_p {\bf q}}{\nabla }. \end{equation} \tag{ 14 }$$

where ∇ = β_TP/ρC_PT, a factor that we will see again. Except possibly for extreme cases, the enthalpy flux F_e, which is intimately related to the stellar luminosity, is itself closely related to the buoyancy flux q and hence to the convective velocity (see also Meakin & Arnett 2007b). MLT ignores this connection; we shall exploit it. The source of turbulent KE is directly proportional to the convective luminosity and thus to the radial entropy gradient ("superadiabatic gradient") of MLT.

Buoyancy is one of the two dominant terms, and with a sign opposite to that of the damping; it differs in magnitude from the inferred damping at the convective boundaries (lower) and the broad peak above the bottom of the convection zone (higher).

The lower panel shows the isotropic Kolmogorov approximation to the local damping, M_CZ u³_rms/ℓ_D = 1.54 M_CZ u³_t/ℓ_D. The condition of global balance between all driving and damping terms gives a value for ℓ_D, but the local value of u_t is used. This gives a relatively good, smoothed fit to the actual inferred damping throughout the turbulent convection zone, and departs only in the boundary layers where the velocity field is due primarily to gravity waves. The net effect of the remaining terms in the KE equation is to modify the velocity field so that the damping is more smoothly distributed over the turbulent region than is the driving.

The agreement in Figure 2, lower panel, between the actual inferred dissipation and our estimate provides strong support for our introduction of the "Kolmogorov dissipation" in both its global and local forms.

Figure 3 shows the spatial behavior of the flux divergence terms, $\bf \nabla \cdot F_{K}$ for flux of KE and $\bf \nabla \cdot F_{P}$ for the pressure correlation flux. These terms have significant positive and negative contributions, which cancel upon averaging over radius, and therefore are more important than Table 1 would suggest. The change of sign in a divergence implies that these terms move KE. In particular, they remove KE from the region where buoyancy is strong, and add it to regions where the buoyancy is weak. The divergence of the KE flux (Figure 3, upper panel) is the most significant in transporting energy. It moves KE from the region in which buoyancy driving exceeds the inferred damping, and toward the convective boundaries.

In Figure 2, top panel, we saw that there is negative buoyancy at the convective boundaries; this is due to buoyancy braking of the convective motion. The pressure correlation flux (Figure 3, lower panel) is most effective in the deficit regions right at the convective boundary, and generates elastic response (waves) in the stably stratified regions outside the convection zone. This two-step behavior, with first F_KE and then F_P carrying energy to the edge of the convection zone, is mirrored at the upper convective boundary, but at lower amplitude.

The p'∇ · u' term, shown in the top panel of Figure 4, does the same sort of thing at a much reduced level. In the bottom panel is shown the time derivative of the KE, which is small. The convection is close to a steady-state behavior.

There is a simple picture which explains these trends.

• 1.  
  The extent of turbulence is limited by energy and by boundaries. In a star, turbulence will mix even stably stratified layers so long as sufficient KE is available to supply the work necessary.
• 2.  
  Turbulence takes ordered motion on the large scales and converts it to disordered motion on small scales. It makes the flow more isotropic. Kolmogorov dissipation is derived with the assumption of homogeneity and isotropy, and so becomes a better approximation as turbulence acts.

The convective motion is driven by buoyant acceleration, parallel to the gravitational acceleration vector, and is necessarily large-scale and anisotropic. The large-scale order of this motion is destroyed by turbulence, which spreads through all accessible regions. Thus, over all the convective region, except right at the boundaries, the flow is made more isotropic, and the Kolmogorov damping becomes better approximated by the local isotropic expression

$$\begin{equation} \varepsilon = \rho (u^{\prime })_{\rm rms}^{3}\big/\ell _{D}, \end{equation} \tag{ 15 }$$

(per unit volume). After carefully distinguishing between global Kolmogorov damping (Equation 5) and the local version (Equation 15), we find that turbulence tends to drive fluxes in such a way as to make both valid.

Having explored the spatial dependence of the KE equation, we now examine its time dependence. The flow, while wildly fluctuating, has an orderly statistical behavior. To see this, we examine the behavior of KE, integrated over the volume of the convection zone, as a function of time. The integral buoyancy flux, $q_{\rm int}=\int \big\langle {\bf g \cdot u^{\prime } }{\rho ^{\prime } \over \rho } \big\rangle dr$, has dimensions of velocity cubed. It is convenient to plot KE and q^3/2_int on the same graph for detailed comparison. Figure 5 shows the time behavior of damping and driving terms. There are two types of time dependence: (1) a set of bursts (pulses), and (2) a secular increase in amplitude (related to the lack of total balance between nuclear heating and neutrino cooling over the convective zone (Arnett 1996).

Zoom In Zoom Out Reset image size

For a static steady state, these two curves would be almost identical. If for simplicity we assume a planar convective region and neglect the smaller terms in the KE equation (Equation 1) in favor of buoyancy and damping, we find the simple result ${v_{\rm rms}^{3}}\big/\ell _{D}\approx { 1 \over \ell _{\rm CZ}}\int q dr$, or ℓ_D/ℓ_CZ ≈ v³_rms/∫qdr.

The power spectra of both curves have a peak at 89 s. Both exhibit strong fluctuations; the buoyancy flux precedes the KE by roughly 20 s. We identify⁹ this lag with the time it takes for the turbulent cascade to react to changes in the flow. This is suggestively close to our previous estimate of the turbulent decay time of about 14 s.

Buoyancy is primarily a property of the largest scales, while damping is a property of the smallest. The large separation of scales means that they are not tightly coupled (except on average), a characteristic of turbulent systems.

The buoyancy flux reaches a high value before turbulence can stop its rise, so that it overshoots the steady-state condition. This leads to excessive velocities, which then cascade to excessive damping. In our simulations, the lag in dissipation behind buoyant driving aids fluctuations about the nominal steady-state condition.

This "boom and bust" cycle is reminiscent of prey–predator relations and the logistic map (May 1976), which also have chaotic behavior. We find that an iterated map of the delay in turbulent damping is quantitatively inadequate to drive the fluctuations. They are driven by nuclear burning, as our simulation (described below) of turbulent decay implies. It will be interesting to explore whether this fluctuating behavior depends upon resolution (we expect it to be dominated by the largest scale eddies, and weakly dependent on resolution).

Consider the average conditions from 200 to 800 s in the simulations. The average level of turbulent KE is about 0.7 of q^2/3_int. If we ignore the smaller terms in the KE equation, we find

$$\begin{equation} \alpha _{D}= \ell _{D}/\ell _{\rm CZ}= 1.54\ u_{t}^{3} \bigg/\!\! \int q dr \approx 0.89. \end{equation} \tag{ 16 }$$

This is roughly what we get from Table 1, so we conclude that to our accuracy, the dissipation length in our simulations is roughly the depth of the convective zone, consistent with Kolmogorov theory. The appearance of ℓ_CZ here is due to integration over the convection zone, not to any assumption about the nature of the largest eddies.

The decay of turbulence is an interesting and complex topic in its own right (Batchelor 1960). Here we will utilize it to provide independent estimates of the damping length. These estimates will differ in detail from those of driven convection because the flow details change, and the dissipation is a function of the flow properties. Nevertheless the sizes of the damping lengths are found to be comparable.

In the left panel in Figure 6, we show an independent simulation: the rate of decay of turbulent KE, after oxygen burning and neutrino cooling are artificially turned off at 439 s (in the middle of the previous simulation).

Zoom In Zoom Out Reset image size

If we neglect the small terms for energy flux escaping the convective zone, we may express the KE equation (Equation 1) more concisely as D = dK/dt − B, where D is the inferred dissipation. After 30 s, the buoyancy driving term B becomes small, and dK/dt tracks the turbulent dissipation D. Note that dK/dt is slightly below D, possibly because of entrainment of stable matter at the convective boundary which requires energy. If we globally fit the inferred damping term by ε_K = −M_CZ v³_rms/ℓ_D, we have ℓ_D = 4.16 × 10⁸cm = 0.97 ℓ_CZ.

The right panel in Figure 6 shows a fit to the KE, starting 40 s after the burning was switched off, and the fossil buoyancy had died away. The KE may be represented analytically if dK/dt = D, so that this late part of the damping curve is reasonably well approximated by a lower damping length, ℓ_D = 2.6 × 10⁸cm = 0.61 ℓ_CZ.

This independent simulation supports the identification of subgrid damping with isotropic Kolmogorov damping, and the representation of damping by the global expression (Equation 5), and ℓ_D ≈ (0.6 to 1.0)ℓ_CZ.

These and the previous estimates for ℓ_D seem to vary beyond the statistical accuracy of the analysis. Replacing the complexity of turbulent dissipation by the simple expression ε_K = −M_CZ v_rms³/ℓ_D, with constant ℓ_D is a vast change. This quantity ℓ_D is a property of the flow, and indirectly of the problem being addressed. In particular, the pure decay is more isotropic than the driven cases (Ogilvie 2003). As Table 3 indicates, the ratio ℓ_D/ℓ_CZ is different for different flows, but of order unity. Further investigation with different simulations and different physical parameters is needed to understand more precisely the general behavior of ℓ_D.

In MLT a connection is made between superadiabatic gradient Δ∇ and convective velocity. This allows the convective enthalpy flux to be expressed in terms of Δ∇ (see, Kippenhahn & Weigert 1990, Section 7, Clayton 1983, Sections 3–5, and Hansen & Kawaler 1994, Section 5). The turbulent KE equation (Equation 1) can perform the same function, if the convective velocity is identified with the rms velocity $\langle \overline{u^{\prime 2}}\rangle ^{1 \over 2}$ at a given radius, and used in the enthalpy flux in the same way. This replaces the parameterization used in MLT with a physical constraint.

We need to rewrite the total energy Equation (A6 in Meakin & Arnett 2007b) in the form of the internal energy equation used in stellar evolution. To get the internal energy equation, we subtract the KE Equation (A12) from the total energy Equation (A6), and find

$$\begin{eqnarray} \partial _{t} \langle \overline{ \rho E_{I}}\rangle &+& {\bf \nabla \cdot }\langle \overline{ {\bf u_{0}} (\rho E_{I} + p_{0})}\rangle\nonumber\\ = &-&\!\nabla \cdot [{\bf F_{I} + F_{r} }] - \langle \overline{ p^{\prime } \nabla \cdot {\bf u^{\prime }} }\rangle \nonumber \\ &+& \varepsilon _{K} + \langle \overline{ \rho _{0} {\bf u_{0} \cdot g }}\rangle +\langle \overline{ \rho \epsilon }\rangle. \end{eqnarray} \tag{ 17 }$$

The left-hand side is just ρ(dE_I/dt + p₀dV/dt) in Lagrangian (comoving) coordinates (Landau & Lifshitz 1959), defining ρ₀ = 1/V, so that

$$\begin{eqnarray} &&\hspace{-1.8pc}d\langle \overline{E_{I}}\rangle /dt + \langle \overline{p_{0} dV/dt }\rangle = - { 1 \over \rho _{0}}{\bf \nabla\ \cdot}\ [ {\bf F_{e} + F_{r} } ] +\langle \overline{\epsilon } \rangle\nonumber\\ &&\hspace{-1.8pc}\quad - \left\langle \overline{{ p^{\prime }\over \rho _{0}}{\bf \nabla \cdot u^{\prime } }}\right\rangle + \varepsilon _{K}/\rho _{0} + \left\langle \overline{ {\bf u_{0}\ \cdot}\ \left({\bf g}\,-\,\!{1\over \rho _{0} }{\bf \nabla }p_{0}\right) }\right\rangle. \end{eqnarray} \tag{ 18 }$$

The first line contains usual terms. Here F_r is the heat flux due to radiative diffusion, and in this frame the flux of internal energy F_I becomes the enthalpy flux carried by convection F_e (see, Tennekes & Lumley 1972, p. 33, and our Equation 12). The nuclear heating term includes neutrino cooling.

We will rewrite Equation (18) in a more familiar notation (see Arnett 1996, Clayton 1983, Hansen & Kawaler 1994, or Kippenhahn & Weigert 1990),

$$\begin{eqnarray} &&\hspace{-1.8pc}dE/dt + PdV/dt = -{1 \over \rho }{\bf \nabla \cdot [ F_{e} + F_{r} ] } + \epsilon - \left\langle \overline{ { p^{\prime }\over \rho _{0}}{\bf \nabla \cdot u^{\prime } } }\right\rangle \nonumber\\ &&+ \varepsilon _{K}/\rho _{0} + \left\langle \overline{ {\bf u_{0}\ \cdot}\ \left({\bf g } -{1 \over \rho _{0} }{\bf \nabla }p_{0}\right) }\right\rangle. \end{eqnarray} \tag{ 19 }$$

The first four terms contain the usual formulation. The new terms are: $\big\langle \overline{{ p^{\prime }\over \rho _{0}} {\bf \nabla \cdot u^{\prime }} }\big\rangle$ which represents the compressional work done by pressure fluctuations (which also appears in the KE equation and we have seen it to be small here), ε_K/ρ₀ which is the deposition of heat by the Kolmogorov turbulent cascade, and $\big\langle \overline{ {\bf u_{0}\, \cdot}\,\big({\bf g} -{1 \over \rho _{0} }{\bf \nabla } p_{0}\big) }\big\rangle$, which is zero in hydrostatic equilibrium without rotation or expansion, and drives meridional circulation for rotating, radiative stars (Tassoul 1978; Clayton 1983).

The ε_K/ρ₀ term allows turbulent KE to do dissipative heating, a new effect not in conventional stellar evolution, and it is not guaranteed to be negligible. Through Equation (1) it couples the divergence of the KE flux (including rotational shear) and the wave energy flux to the internal energy. As an internal energy source term, it can generate entropy and cause mixing even in radiative regions. For the Sun, this term would give rise to heating in the photosphere, chromosphere, and corona by shocks and wave motion (pressure and Alfven waves), for example.

Equations (1) and (19) represent an extension of turbulent convection theory for stellar evolution, in which (1) the algebraic relation for convective velocity and superadiabatic gradient is replaced by a differential equation (see Spiegel 1972), and (2) the Kolmogorov cascade is explicit in the formulation. Note that both space and time derivatives appear, making the system nonlocal and time dependent, unlike MLT. These derivatives allow the treatment of boundary dynamics in a physical way. We will explore specific implementations in a subsequent paper, including entrainment at convective boundaries.

Let us now consider some simplifications in order to clarify the meaning of these equations. In particular, we assume time invariance (so the time derivatives are small on average), no overall background motion (u₀ = 0), and little work done by pressure perturbations on the velocity perturbations. As shown previously, these are reasonable approximations except in extreme cases. Then Equation (1) becomes

$$\begin{equation} \overline{ \langle \rho ^{\prime } {\bf g \cdot u^{\prime }}\rangle } - {\bf \nabla \cdot (F_K + F_P) } = \varepsilon _K, \end{equation} \tag{ 20 }$$

and Equation (19) becomes

$$\begin{equation} \rho \epsilon - {\bf \nabla \cdot (F_e + F_r) } = - \varepsilon _K. \end{equation} \tag{ 21 }$$

The coupling of internal energy and turbulent KE occurs through the Kolmogorov cascade, which creates internal energy by damping KE. Eliminating ε_K, we have

$$\begin{equation} \rho \epsilon = - \overline{ \langle \rho ^{\prime } {\bf g \cdot u^{\prime }}\rangle } + {\bf \nabla \cdot (F_e + F_P + F_r + F_K) }. \end{equation} \tag{ 22 }$$

In a steady state, the nuclear energy generation must balance the divergence of fluxes out of the region and supply the work needed to maintain the convective flow. This is different from the usual formulation used in stellar evolution in that there are new terms $(\!- \overline{ \langle \rho ^{\prime } {\bf g \cdot u^{\prime }}\rangle } + {\bf \nabla \cdot (F_K + F_P) }$), which combine to equal ε_K, the dissipation rate due to the Kolmogorov turbulent cascade. If the turbulent velocity is nonzero, these corrections are also nonzero, leading to the conclusion that the standard formulation of stellar evolution is wrong in neglecting turbulent heating. The motion implied by convection will carry KE, drive pressure and gravity waves into radiative regions, and give local microscopic heating as it dissipates in a transit time, effects that should no longer be ignored in stellar evolution. Through its dependence on ζ, for a given convective enthalpy flux, the KE flux is dependent upon the equation of state.

For a low Mach-number flow,

$$\begin{equation} \overline{ \langle \rho ^{\prime } {\bf g \cdot u^{\prime }}\rangle } = \beta _T g F_e / C_P T, \end{equation} \tag{ 23 }$$

and radial coordinates so that ∇ · F → dL/dm, then Equation (22) becomes

$$\begin{equation} \epsilon = -{\nabla \over m_p}L_e + {d \over dm}(L_e + L_P+ L_r + L_K), \end{equation} \tag{ 24 }$$

where ∇ = β_TPV/C_PT, m_p = 4πr²ρH_P, and H_P is the pressure scale height. For an ideal gas, ∇ → (γ − 1)/γ, or 0.4 for γ = 5/3. For gases with a specific heat at constant pressure which is large compared to the specific heat at constant volume (such as partially ionized plasma or electron–positron plasma), ζ is considerably smaller (see below).

We now examine each of a set of important convection parameters (see Porter & Woodward 2000 and Meakin & Arnett 2007b for details). These parameters reflect the various uses of the MLT parameter α, and are a convenient and concise way to compare the simulations of compressible turbulent convection.

4.1.1. α_T

In 3D simulations, a correlation is found between rms temperature fluctuation $(T^{\prime })_{\rm rms}= \langle \overline{T^{\prime 2}}\rangle ^{1 \over 2}$ and the superadiabatic gradient Δ∇ = ∇ − ∇_a − ∇_x, using the conventional notation in astrophysics of ∇ = ∂ln T/∂ln P (Kippenhahn & Weigert 1990; Hansen & Kawaler 1994; Clayton 1983). Here, ∇_a = ∂ln T/∂ln P, taken at constant entropy and composition, and ∇_x = −∇_μ is the remaining part due to possible compositional change. Then,

$$\begin{equation} (T^{\prime })_{\rm rms}/T_{0} = \alpha _T \Delta \nabla, \end{equation} \tag{ 25 }$$

on average, over time and over the volume of the convection zone. Consider low Mach-number flow so that ρ'/ρ₀ = −β_T(T'/T₀).

In a single convective roll (e.g., Lorenz 1963; Tritton 1988), T' is the temperature perturbation amplitude at the horizontal midplane. In the vertical midplane of the roll, T' = [(dT/dr) − (dT/dr)_ad]ℓ_CZ/2 is the corresponding amplitude. Assuming the amplitudes are comparable,

$$\begin{equation} T^{\prime }/T_{0} = (\ell _{\rm CZ}/2 H_{P}) \Delta \nabla, \end{equation} \tag{ 26 }$$

implying that α_T ≈ ℓ_CZ/2H_P, and not a constant. In MLT, α_T = α/2, which is a particular choice for the assumed flow.

In the simple picture of a single convective roll, (T')_rms/T₀ can give a buoyancy torque, while Δ∇ cannot, because the gravitational acceleration vector is directed radially downward (Tritton 1988). In the Lorenz (1963) model of thermal convection, the difference between the two is intimately connected to the onset of chaotic behavior. Meakin & Arnett (2007b) find α_T ≈ 0.73 (see their Figure 17), so that

$$\begin{equation} \alpha _T = 0.73\ \ell _{\rm CZ}/2 H_{P}. \end{equation} \tag{ 27 }$$

There is considerable variation in the values of α_T shown in Table 4, with a tendency to increase for deeper (more stratified) convection zones.

Note that for two convective rolls, one atop the other, we would expect the characteristic roll size to change (ℓ ≈ ℓ_CZ → ℓ_CZ/2, so α_T → α_T/2, approximately). This explicitly shows how the convection parameters can be a function of the properties of the flow itself.

Figure 7, left panel, shows the behavior of α_T as a function of the depth of the convection zone, for the computations listed in Table 4. The two PPM calculations (Meakin & Arnett 2007b; Porter & Woodward 2000) agree with the scaling indicated in Equation (27). The calculations using the compressible viscous equations with subgrid modeling for dissipation all lie below Equation (27). In order to give some idea of the depth needed in simulations of stellar convection zones, the x-axis in Figure 7 is marked from zero to the depth of the solar convection zone (20 pressure scale heights). None of these simulations describe such an extremely anisotropic case.

Zoom In Zoom Out Reset image size

4.1.2. α_E and α_K

The enthalpy flux is

$$\begin{eqnarray} F_{c} &=& \rho _{0}C_{P} \overline{T^{\prime }u^{\prime }} \nonumber \\ &=& \alpha _{E} \rho _{0}C_{P} (T^{\prime })_{\rm rms}\ (u^{\prime })_{\rm rms}, \end{eqnarray} \tag{ 28 }$$

where Meakin & Arnett (2007b) find α_E = 0.70 ± 0.03; in Table 4, α_E is relatively constant among the simulations (the total range is about 20%). It is not ruled out that α_E might be a universal constant, or at least slowly varying. Note that α_E is just the correlation coefficient for T' and u'. It seems that α_E is not sensitive to the Prandtl number Pr. Porter & Woodward (2000) have a different Pr (ours is Pr ≈ 1) and get an α_E similar to ours.

Using Equation (25), this becomes

$$\begin{equation} F_{c} = \alpha _{E} \alpha _{T} \rho _{0}C_{P}T_{0} (u^{\prime })_{\rm rms}\Delta \nabla. \end{equation} \tag{ 29 }$$

Similarly, if the KE flux is

$$\begin{equation} F_{K} = \case{1}{2} \langle \rho u^{\prime 2} u_{z}^{\prime } \rangle, \end{equation} \tag{ 30 }$$

we may define

$$\begin{equation} \alpha _{K} = \langle \rho u^{\prime 2} u_{z}^{\prime } \rangle / \langle \rho \rangle \langle (u_{z}^{\prime })^{2} \rangle ^{3/2}, \end{equation} \tag{ 31 }$$

so that $F_{K} = { \alpha _{K} \over 2} \langle \rho \rangle \langle (u_{z}^{\prime })^{2} \rangle ^{3/2}$. Note that the sign of α_K can be negative.

4.1.3. α_v

In Meakin & Arnett (2007b), the correlation between convective velocity (squared) and superadiabatic excess Δ∇ is written as

$$\begin{equation} (u^{\prime })_{\rm rms}^2 = \big(\alpha _v^2 \big/ 4 \big) g \beta _{T} H_P \Delta \nabla. \end{equation} \tag{ 32 }$$

In MLT, a similar expression is defined, with α²/8 replacing α²_v/4. If we have local balance between buoyancy driving and turbulent damping,

$$\begin{equation} \langle \rho ^{\prime } g u^{\prime } \rangle \approx \langle \rho \rangle (u^{\prime })_{\rm rms}^{3}\big/\ell _{D}. \end{equation} \tag{ 33 }$$

Using Equations (13), (25), and (28), we have

$$\begin{equation} (u^{\prime })_{\rm rms}^{2} = (\ell _{D} \alpha _{T}\alpha _{E}) \beta _{T} g \Delta \nabla. \end{equation} \tag{ 34 }$$

Comparing this to Equation (32) we have

$$\begin{equation} \alpha _{v}^{2}\big/4 = \alpha _{T}\alpha _{E} \ell _{D} / H_{P}. \end{equation} \tag{ 35 }$$

Using Equation (27) and α_E = 0.70,

$$\begin{equation} \alpha _{v} = 1.22\ (\ell _{\rm CZ}/2H_{P})(\ell _{D}/0.9\ell _{\rm CZ})^{{1\over 2}}. \end{equation} \tag{ 36 }$$

In Table 4, we see that α_v is variable, and tends to increase with increasing depth of the convective zone. This is shown explicitly in the right panel of Figure 7. The two PPM calculations are in excellent agreement with Equation (36), taking ℓ_D = 0.9 ℓ_CZ (α_D = 0.9), while the three compressible viscous calculations lie below it.

Using Equation (29), we have

$$\begin{equation} F_{c} = (\alpha _E \alpha _T \alpha _v /2) \rho _{0} C_{P} T_{0} \sqrt{g \beta _{T} H_P }(\Delta \nabla)^{3/2}. \end{equation} \tag{ 37 }$$

Using Equation (27), Equation (36), ℓ_D = 0.9ℓ_CZ, and α_E = 0.70, we have

$$\begin{equation} \alpha _{E}\alpha _{T}\alpha _{v}/2 = 0.312 (\ell _{\rm CZ}/H_{P})^{2} \end{equation} \tag{ 38 }$$

for the factor in Equation (37).

The thermodynamic character of the electron–positron plasma in the oxygen-burning shell is significantly different from an ideal gas, an effect which must be taken into account in comparisons to simulations which use different equations of state. This is illustrated in Table 5; we use the Helmholtz equation of state of Timmes & Swesty (2000). The entropy in the oxygen burning convection zone is $S/{\cal R} \approx 4.6$ in dimensionless units, where $\cal R$ is the gas constant. The zone is about 2 pressure scale heights deep (see Column 3). For an ideal gas, the specific heat at constant pressure is $2.5\cal R$; $C_{P}/{\cal R}$ is much larger for the plasma, ranging between 10 and 16 (Column 4). Column 5 gives the value of β_T, which is unity for an ideal gas, but lies between 3 and 4 for the plasma. Column 6 gives the ratio of PV/C_PT = p₀/ρ₀C_PT₀, which is 0.4 for an ideal gas. The ratio of the buoyancy flux to the enthalpy flux (Equation (14)) is proportional to ∇ = β_TPV/C_PT. The same factor appears in the ratio of the KE flux to the enthalpy flux. For the ideal gas ζ = 0.4, but is smaller for the plasma (Column 7). For the same convective enthalpy flux, the pair-plasma has a much lower KE flux, so the velocity scale is lower.

Table 5. Thermodynamic Parameters for Entropy $S/{\cal R}=4.623$

T9	ρ6	ln Pmax/P	$C_{P}/{\cal R}$	βT	PV/CPT	∇ad
2.51	1.580	0.0	15.36	3.912	0.0609	0.2382
2.31	1.225	0.348	14.53	3.757	0.0638	0.2397
2.11	0.930	0.724	13.62	3.592	0.0673	0.2417
1.91	0.690	1.134	12.67	3.421	0.0714	0.2443
1.71	0.498	1.583	11.73	3.257	0.0762	0.2483
1.51	0.348	2.080	10.83	3.106	0.0814	0.2529

Download table as:  ASCIITypeset image

The KE flux is averaged over angle at a given radius, and averaged over two convective turnover times. As Figure 5 shows, this time span covers significant dynamic behavior. Roughly speaking, an unstable configuration forms, becomes dynamic, reforms, becomes dynamic again, and so on. Over this turnover timescale the nuclear burning provides the energy necessary to drive the turbulent KE $\big(\epsilon \tau \approx {1 \over 2}\langle \overline{u^{\prime 2}}\rangle\big)$.

The convective instability is closely related to the Rayleigh–Taylor instability (Chandrasekhar 1961a), which has been produced dramatically in high energy–density plasma experiments (Remington et al. 1999), i.e., under star-like conditions and well into the nonlinear growth regime. The characteristic behavior is the rise of mushroom-shaped higher entropy plasma and the fall of spike-shaped lower entropy plasma. If we consider a closed volume, these motions are accompanied by slower return motions which maintain mass conservation. The KE flux scales with velocity cubed, and so is dominated by the fast moving mushrooms and spikes. In a symmetrical system, we will have an upward KE flux (from the mushrooms) in the region above the horizontal midplane of the volume, and a downward KE flux in the region below. If we average over several cycles (turnover times) the KE flux will be dominated by these motions, being positive (upward) above the midplane and negative below. This is qualitatively similar to the simulation results of Meakin & Arnett (2007b; see their Figures 21 and 22).

This simple picture is complicated by an up–down asymmetry due to stratification (hydrostatic structure). The depth of the convection zone in pressure scale heights, ℓ_CZ/H_P = ln(P_bot/P_top), is a direct measure of the up–down asymmetry. Upward flows move into regions of lower pressure, and expand; downward flows move into regions of higher pressure, and are compressed. The smaller area of the downflows implies higher velocities relative to a coordinate frame containing the convection zone (a Lagrangian frame). This will favor downward (negative) KE fluxes. The neglect of KE fluxes (Böhm-Vitense 1992) corresponds to the limiting case of a shallow convection zone. For deep convection zones, there is a strong bias in favor of fast downward plumes (Nordlund & Stein 2000; Stein & Nordlund 1998), and these dominate the flow for simulations with ℓ_CZ/H_P ⩾ 4 or so. In Meakin & Arnett (2007b), which has a depth ℓ_CZ/H_P = 2, the convective zone is skewed, so that the surface, which separates the positive and the negative KE fluxes, moves nearer to the bottom of the convection zone. For Porter & Woodward (2000), where ℓ_CZ/H_P is larger, the positive KE flux is overwhelmed, and the direction of the KE flux is opposite to that of the much larger enthalpy flux. We expect this to be a general property of deep convection zones.

There is a natural limit to the depth of convection zones expected in active burning regions. Entropies in active burning regions vary relatively slowly (Arnett 1996), so that the depth of a convection zone implies a value of the temperature ratio between bottom and top. In Table 5, the convection zone extends down almost to neon burning temperatures (T ≈ 1.5 × 10⁹K). Further extension will entrain new fuel into the oxygen convective shell, which will burn at the top of the convection zone, choking the flow. For the last stages (C, Ne, O, and Si burning), the burning temperatures for different fuels are fairly close together, implying that their related convection zones will tend to be relatively shallow. They are also close together in radius, raising the issue of interactions between them (Meakin & Arnett 2006).

Here $\ell _{\rm CZ}/H_{P} =\ln P_{\rm bot}/P_{\rm top} = {\gamma \over \gamma -1 } \ln T_{\rm bot}/T_{\rm top}$. For example, for advanced burning stages or radiation dominated regions, γ ≈ 4/3, and γ/(γ − 1) = 4. For helium and hydrogen burning, T(He)/T(H) ≈ 13 so ℓ_CZ/H_P ⩽ 4ln 13 ≈ 10. A helium burning convective zone will not be deeper than about 10 scale heights unless the overlying matter is devoid of hydrogen. For carbon burning the corresponding depth is ℓ_CZ/H_P ≲ 3, unless devoid of H and He fuel. Surface convection zones may extend down to the hydrogen burning regions. Very roughly, ln T(H)/T_e ≈ 8, so ℓ_CZ/H_P ≲ 32, using the structure of an n = 3 polytrope.

Table 6 gives the ratio of kinetic energy flux to enthalpy flux for several 3D simulations. This ratio is much smaller in our simulations (by a factor of ∼20–40) than in the others, and F_K changes sign in our convection zone.

As in Meakin & Arnett (2007b), we can write the KE to enthalpy flux ratio according to mixing length relationships,

$$\begin{equation} \frac{F_{K}}{F_{c}} \sim \frac{\rho v_c^2\big/2}{\rho c_p T^{\prime }}\frac{v_c}{v_c} \sim \left(\frac{\alpha _v^2}{8 \alpha _T}\right) \nabla _{ad}, \end{equation} \tag{ 39 }$$

and then balance between buoyancy driving and turbulent damping through Equation (37) gives

$$\begin{equation} \frac{F_K}{F_c}\sim \frac{\alpha _D\alpha _E}{2}\left(\frac{l_{CZ}}{H_P}\right). \end{equation} \tag{ 40 }$$

See Equation (3.14) of Porter & Woodward (2000), which uses a gamma-law equation of state to generate the sum of kinetic and enthalpy fluxes, and implies an equivalent flux ratio.

The ratio of ∇_ad for the ideal gas to that of the electron–positron gas gives a factor of 1.6 or so. The ratios of the depth of the convection zones give another factor of ∼5. While these considerations illustrate the role played by the depth of the convection zone and go some way towards explaining the trends in KE flux it is also important to consider the geometry of the driving region, which relates to how well a local balance between buoyancy driving and turbulent damping is achieved as assumed in Equation (37). In particular, the length scale l_a over which buoyancy is imparted to the stellar plasma through either heating at the base of the convection zone or cooling at the top can affect the KE to enthalpy flux ratio dramatically. This may be understood simply: the flow depends both upon the geometry of the convective domain and upon the way in which the fluid is stirred. At present we have at least two basic patterns, a mostly negative kinetic energy flux for deep convective zones driven by surface cooling (most stellar surface convection zones) and a more complex positive-negative convective flux for shallower zones driven by nuclear burning. More simulations are underway to clarify this issue (C. Meakin & D. Arnett 2009, in preparation).

Are there limits to the linear rise in energy of convection that is implied by Figure 7? Deep convection zones (ℓ_CZ ⩾ 3H_P) have strong negative KE fluxes. For very deep zones, the extrapolated kinetic fluxes imply supersonic velocities. Such large velocities would generate shock waves, which would increase dissipation, converting KE into internal energy. The rate scales as the velocity difference cubed, which is the same scaling as turbulent dissipation in the Kolmogorov cascade (Bethe 1942; Boris 2007). This suggests that even if deeper convection zones did tend to have increasingly strong velocities, shock dissipation will become significant, and resist the increase of ℓ_D with increasing ℓ_CZ. In this sense, the increase in damping length must "saturate" with increasing depths.

Increased damping may come sooner from another effect: a change in the nature of the eddies and the size of the damping length. Physically this would occur as follows: as the deep, fast, narrow downflows drive into the convection zone they will give rise to shear instabilities at their interfaces, which will lead to mixing with the ambient fluid, and exchange momentum through a turbulent viscosity, and eventually completely dissolve into the background. This would give shorter dissipation lengths, and would begin to occur before the Mach numbers become large enough for significant shock formation.

We expect the PPM calculations for deeper convective zones to "flatten" in Figure 7 as the Yale simulations do, but at higher amplitude, as a result of higher damping rates. This hypothesis needs to be tested numerically, which will be challenging. Shock waves may cause errors at grid boundaries, deep convection zones will have longer thermal relaxation times, and maintaining sufficiently wide aspect angle implies many computational zones for adequate resolution, for example.

The simulations of the Yale group (Chan & Sofia 1989, 1996; Kim et al. 1995, 1996; Robinson et al. 2004), which use fewer zones and a strong damping (Smagorinsky 1963), appear to have a larger dissipation than the PPM simulations (Porter & Woodward 2000; Meakin & Arnett 2007b). This is qualitatively equivalent to the saturation discussed above, and may be tested if calculations using the Yale code (or equivalent) are repeated at higher resolution and/or lower dissipation.

The energy flux terms include both advective transport and waves. Here we will recall some properties of waves and their generation (Landau & Lifshitz 1959; Press 1981). The characteristic frequencies of convective motion are centered about a frequency ω_CZ ≈ v_rms/ℓ_CZ. The convective Mach number, $(u^{\prime })_{\rm rms}/{\cal C} = \rho ^{\prime }/\rho _{0} = p^{\prime }/\rho _{0} {\cal C}^{2}$, is a measure of compressibility of the flow; here ${\cal C}$ is the sound speed. If the interface between convective and radiative zones moves with the matter (is volume conserving, on average), it generates acoustic waves by dipole emission at a luminosity L_p-mode ∝ ω⁶_CZ (see, Landau & Lifshitz 1959, Section 73), or as the Mach number to the sixth power. For more vigorous motion, the perturbation may give volume changes, so that the acoustic wave generation by monopole emissivity, or Mach number to the fourth power. For subsonic flows this channel is closed to significant energy flow, but open as (u')_rms approaches the local sound speed, as it does in the surface convection zones of many stars, including the Sun, or as it may in the stage prior to core collapse in massive stars.

While the exact power dependence may depend upon specific geometries and degree of interference, a general result seems to be: the gravity wave channel dominates over the acoustic channel for low Mach-number flows (as we observe in our simulations; see Meakin & Arnett 2007a for a more detailed discussion of both the g-mode and p-mode behaviors, including mixed p- and g-modes). Both types of waves are generated by convection interacting with stably stratified bounding regions, and the luminosity of each depends upon both the convective vigor and the impedance at the boundary. Such waves can propagate into stably stratified regions (Press 1981; Press & Rybicki 1981; Young & Arnett 2005). Generation of gravity waves by convective turbulence has become an issue in questions of mixing and angular momentum transport (García López & Spruit 1991; Charbonnel & Talon 1999; Talon & Charbonnel 2004; Young & Arnett 2005).

The establishment of a robust estimate of the turbulent velocity field should improve estimates of wave generation, which has implications for mass loss, mixing in radiative regions, coronal heating, and helioseismology. In particular, the pressure correlation flux takes over the energy carried by the KE flux as the convective boundaries are approached. This gives a direct connection between the scale of turbulent velocity and g-mode wave emission (Press 1981).

A closely similar set of mean-field equations is used in the theory of the magnetic resonance instability (MRI) in accretion disks (Balbus & Hawley 1998; Pessah et al. 2006). If we had included magnetic fields in the MHD approximation, and assumed strong rotation, our procedure would have produced an equation for mechanical energy(our Equation (1) corresponds to Equation (17) of Balbus & Hawley 1998), and total energy (Equation (A6) of Meakin & Arnett 2007b to their Equation (27)). Projection onto a cylindrical coordinate system, with the rotational axis oriented parallel to the total angular momentum vector, would give an angular momentum equation (their Equation (29)). Auxiliary equations provide for magnetic field dynamics, radiation transport, and nuclear burning. This underlying similarity provides a way to write a more general set of mean-field equations, of which both stars and accretion disks are limiting cases. In turn this allows a systematic evaluation of the relative importance of different effects (rotational mixing and convective mixing, for example) which are now considered piecemeal.

Our simulations include the complete set of rotational terms, but the initial conditions imply that these terms are not exercised except as convectively induced shear. Our simulations do not include magnetic fields in the MHD approximation, but could be generalized to do so (Stone & Gardiner 2007; Pessah et al. 2006). Because they interact, rotation and magnetic fields should be included together.

Balbus & Hawley (1994) have argued that in the stellar case, the weak-field MRI dominates over merely hydrodynamic instabilities, and drives the radiative zones (but not convective zones) toward solid body rotation. Heger et al. (2000) argued the reverse (based on the Høiland instability criterion), that convective zones would tend toward rigid body rotation, and radiative zones would tend to have differential rotation. Because the Høiland instability applies to neutral fluids, not dense plasma, it is probably not relevant for stars. Helioseismology (Thompson et al. 1996; Howe et al. 2005; Brandenburg 2007) is showing that while the convective zone of the Sun shows differential rotation, the underlying radiative zone seems to be tending toward solid body rotation, as the MRI arguments suggest. Browning & Basi (2007) have shown that even deep convection zones can generate significant magnetic fields, so that the presence of a stable interface is not necessary for field generation. The rotational state of the central regions of the Sun probably depends upon the efficacy of angular momentum transport by processes related to flow of the plasma, including g-mode waves (Charbonnel & Talon 1999), as well as magnetic fields.

In general, it is appropriate to decompose the equations in temperature (T =  T₀ + T') and composition as well. The opacity and the nuclear reaction rates are often sensitive to both. Meakin & Arnett (2006) found flashing due to oxygen burning in vigorous downdrafts which were fuel-rich (the flashes were too mild to affect the flow dramatically). None of these effects are included in standard stellar evolution theory. For simplicity we suppress this complication for the moment; this means that our opacities and reaction rates are to be interpreted as averages over fluctuations in these variables as well. Further investigation of this issue, with 3D simulations, is desirable.

For this set of simulations, the heating and cooling times are at least 100 times longer than the transit times, so that we are in the regime of small Damköhler number D_m ≲ 0.01 (Zel'dovich et al. 1985; Oran & Boris 1987). The release of nuclear energy during a transit time is small relative to the internal energy in the convection zone, but comparable to the superadiabatic energy and to the turbulent KE. The burning drives the turbulent motions, but only gives moderate pulses of KE, as shown in Figure 5. This is unlike the much more complex problem of Type Ia supernova models, for which D_m is large. In our case the necessary averaging over convective cycles does not seem to be a problem.

This convenient state may not apply to the double shell flash stage for asymptotic giant branch (AGB) stars, in which wave-driven mixing and entrainment are likely to complicate the issue of mixing, and therefore figure into the question of s-process nucleosynthesis (Busso et al. 1999; Campbell & Lattanzio 2008).