THE MASS MIXING LENGTH IN CONVECTIVE STELLAR ENVELOPES

For stellar structure and evolution calculations (e.g., Iben 1967), models for the relevant energy-transport mechanisms are needed. In the interior, radiative or conductive transfer of energy is trivial, although the atomic physics going into the opacities and conductivities are not. A satisfactory description of convection in a manner amenable to stellar structure calculations, however, has proven elusive. Stellar modelers usually use the mixing length theory (MLT; Böhm-Vitense 1958; Gough 1977), which despite its simplicity and many shortcomings has been highly successful in matching a variety of stellar observations. In this model, convection is thought of as a collection of fluid parcels that travel a distance equal to the mixing length and then thoroughly mix with the surrounding ambient medium, sharing their energy. The MLT picture is completely symmetric, so the upward flow of buoyant warm elements is accompanied by an equal number of cooler and overdense elements moving inward. In these calculations, the mixing length is a free parameter and is specified as a multiple of the pressure scale height H_P = P/ρg, where P is the pressure, ρ is the density, and g is the gravitational acceleration. In general, there are also free parameters specifying the geometry of the fluid parcel—the ratio of its horizontal to vertical size and another specifying the lateral radiative transfer between convective elements and the ambient medium. The simple mixing length model is used due to a lack of better alternatives and because it is quick to compute. Its primary function is to determine the entropy of the deep, efficient, and nearly adiabatic part of the convection zone in the star. If the magnitude of the mixing length and the other parameters could be specified a priori, simple MLT would be capable of fulfilling this function. However, there is little guidance from theory as to their values.

In addition, the mixing length model suffers from several other problems: it does not allow for overshooting motions that mix material into the surrounding stable layers, it does not allow for turbulence, it does not allow for asymmetry between upflows and downflows, especially the large temperature difference near the surface and its effect on opacities and radiative transfer. These are all phenomena that, according to realistic three-dimensional simulations, contribute to the atmospheric entropy jump and hence, the adiabat of the deep, efficient convection. Several attempts have been made to improve the mixing length model. Some have taken the mixing length to be the distance to the nearest boundary (Stothers & Chin 1997), some have introduced a two-stream model with separate properties of the upflows and downflows (Nordlund 1976), and some have introduced the idea of a spectrum of sizes of the convective cells (Canuto & Mazzitelli 1991).

Several attempts have been made to determine the appropriate value of the ratio of mixing length to pressure scale height by matching theoretical and observed evolutionary tracks in the Hertzsprung–Russell (H-R) diagram (Chieffi et al. 1995; Stassun et al. 2004), by matching observed and model properties of binary systems (Fernandes et al. 1998, 2002; Montalbán et al. 2004), and red giant temperatures (Ferraro et al. 2006). There have also been a few attempts to determine the mixing length parameter by comparing one-dimensional mixing length stellar model calculations with two- and three-dimensional numerical simulations of stellar convection by either matching entropy in the adiabatic portion of the convection zone (Ludwig et al. 1999; Freytag & Salaris 1999) or by matching temperature and density at a given pressure in the interior (Trampedach et al. 1999; R. Trampedach et al. 2011, in preparation).

MLT works reasonably well for determining the thermal structure of convection zones because real mixing of mass occurs in stratified convection. Stars are stratified. Their density decreases from their center to their surface in order to produce a pressure gradient to balance the inward pull of gravity. In order to conserve mass in such a stratified medium, most of the ascending material must turnover and head back down within about a density scale height (Stein & Nordlund 1989). Thus, mixing really does occur on the order of a density scale height.

It is possible to easily and unambiguously calculate the mass mixing length from computational fluid dynamic simulations of stellar convection. Such simulations, when they employ realistic equations of state (EOSs) to relate the pressure to the density and internal energy or temperature and employ sufficiently accurate treatment of radiative transfer to determine the heating and cooling, have as free parameters only the resolution of the computational grid used and the magnitude of the diffusion used to stabilize the calculations. These parameters, however, are not free to be adjusted to match observations. The diffusion is minimized against the requirement of a numerically stable simulation, and the resolution is chosen as a tradeoff between resolving features in the photosphere and computational cost. We would not lower the resolution or increase the diffusion to give a better match to some observations.

In this paper, we present results of such mass mixing calculation for a homogeneous sample of convection simulations, covering a large part of the H-R diagram for which stars have convective envelopes. The simulations form an irregular grid, covering 4500–6900 K on the main sequence and going up to giants with log g = 2.2 and T_eff=4300–5000 K, all with solar composition. The simulations employ the Mihalas–Hummer–Däppen EOS (Hummer & Mihalas 1988), line-opacity from Kurucz (1992), and continuum opacities largely from the Marcs stellar atmosphere package (Gustafsson et al. 1975), but with many improvements. The simulations are described in more detail in Trampedach et al. (2011) and R. Trampedach et al. (2011, in preparation).

The mass mixing length, the distance over which, because of mass conservation, most of the ascending fluid turns over and starts to descend, is the inverse of the logarithmic derivative of the unidirectional mass flux (Stein et al. 2009):

$$\begin{eqnarray} \ell _m(r) & = & |d \ln F_{\uparrow m}(r) / dr |^{-1}\nonumber \\ & = & |d \ln F_{\downarrow m}(r) / dr |^{-1}. \end{eqnarray} \tag{ 1 }$$

This quantity is readily determined from three-dimensional, stellar convection simulations. The equality of the mass flux in the upflows and downflows is just another way of stating mass conservation at each depth. Figure 1 shows the mass flux in the downflows for all of the convection simulations.

Zoom In Zoom Out Reset image size

In the rest of this paper we will look at the upflow, whose mass flux is

$$\begin{equation} F_{\uparrow m}(r) = \langle \rho (r) V(r)\rangle _\uparrow A_\uparrow (r) \approx \langle \rho (r)\rangle _\uparrow \langle V(r)\rangle _\uparrow A_\uparrow (r). \end{equation} \tag{ 2 }$$

The brackets, 〈f(r)〉_↑, denote horizontal averages over the upflows of a quantity f: $\langle f(r)\rangle _\uparrow =\int _{A_\uparrow (r)}f(\vec{s},r) { d^2} \vec{s}/A_\uparrow (r)$, with A_↑(r) being the area with upflows and $\vec{s}$ is the coordinates in the horizontal directions. The approximation of the right-hand side of Equation (2) ignores correlations between density and velocity, but we only use this approximation to gain some insights into the behavior of ℓ_m(r), by re-writing it as

$$\begin{equation} \ell _m(r) \approx \left|\frac{d{\ln \langle \rho \rangle _\uparrow }}{d{r}} + \frac{d{\ln \langle V\rangle _\uparrow }}{d{r}} + \frac{d{\ln A_\uparrow }}{d{r}}\right|^{-1}. \end{equation} \tag{ 3 }$$

Below the surface region, the area of the upflows is nearly constant at about 2/3 of the total area, so its gradient has little effect on ℓ_m. From Equation (3) we see that, ignoring gradients in velocity and filling factor, the mass mixing length is simply the density scale height, H_ρ = dr/dln ρ. Another way to look at it is to consider an upflow along the decreasing density of an average stratification with gradient $\frac{d{\ln \rho }}{d{r}}$. If the upflow had retained its density, it would be overdense by $\Delta \rho /\rho =-\frac{d{\ln \rho }}{d{r}}\Delta r$, compared to its surroundings, after having traveled the distance Δr (remember that the gradient is negative). This is obviously unstable, and instead the fraction Δρ/ρ overturns into the downflows, in order to keep the hydrostatic balance. That fraction is 1 for $\Delta r= 1/\frac{d{\ln \rho }}{d{r}}=H_\rho$ in a linear stratification and 1/e in an exponential stratification. Locally, the upflow will therefore be "eroded" with an e-folding scale of H_ρ. This "erosion" of the upflows does not mean that they lose coherency and disappear over that length scale—on the contrary. We observe coherent upflows over the entire 5.3–7.3 pressure scale heights covered by the convection zones of our simulations (a similar number of pressure scale heights are included above the convection zones of the simulations). We likewise see downdrafts that reach all the way to the bottom of our simulations (in all of our simulations), with weaker ones either stopping or merging into fewer strong downdrafts, before they reach the bottom. The mass mixing length simply describes how only a very small fraction of the mass in the upflows makes it from the interior and all the way up to the photosphere, simply due to mass conservation in a stratified atmosphere.

In stars, the velocity decreases inward from the top of the convection zone, so its derivative is of the opposite sign to that of the density. This decreases the variation of the unidirectional mass flux and increases its scale length compared to the density scale height. In a deep solar simulation, we find a logarithmic derivative of the upward mass flux of ≃4/(5H_ρ) (Figure 3), and the vertical velocity must therefore vary approximately as

$$\begin{equation} V_{\rm vert} \propto \rho ^{-1/5}. \end{equation} \tag{ 4 }$$

The various terms of Equation (3) are shown for the deep solar simulation in Figure 2. We clearly see that the density gradient is the main term and that the velocity gradient produces the deviation from ℓ_m = H_ρ. The gradient of the filling factor is only significant in the surface layers. The difference between the approximation of the right-hand side of Equation (3) and the full expression can be seen by comparison with Figure 3, which is based on the latter.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Near the surface, the mass mixing length is not a constant multiple of either the pressure or density scale heights (Figures 2 and 3). This is due to the abrupt changes in all three mass-flux ingredients in the surface layers. The density gradient becomes very small (for all but the coolest dwarfs the density in the upflows has a sub-photospheric inversion) resulting in a spike in the density's contribution to the mixing length (see the right-hand panel of Figure 2). The small density gradient is however accompanied by a large velocity gradient due to buoyancy braking of the upflows where the density gradient becomes small or inverted, which reduces the density induced peak in the mixing length. In the interior, the density and velocity gradients have opposite signs and a nearly constant ratio. However, as the surface is approached the buoyancy driving increases, which steepens the velocity gradient and gives rise to the wide sup-photospheric bump in the mixing length. The photospheric side of this bump is considerably smoothed by the contribution from the filling factor. Division by the pressure scale height to get the ℓ/H_P presented in Figures 4 and 5 decreases the difference between the atmospheric and the interior values.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

To determine the mass mixing length, ℓ_m, for a simulation, we first calculate the horizontal and time averages of the vertical upward (and downward) momentum (which is the unidirectional mass flux). The unidirectional mass flux increases with increasing depth as the density increases and the velocity decreases, but more slowly than ρ^−1/2. We then calculate the inverse of the logarithmic derivative of this total upward or downward unidirectional mass flux. We have computed the mass mixing length for the set of realistic, compressible, convection simulations by R. Trampedach et al. (2011, in preparation). Figure 1 shows the unidirectional (downward) mass flux for all the stars in our sample. Next, we determine the multiple of the pressure and density scale heights that best match the inverse logarithmic derivative of the unidirectional mass flux between its local maximum at the steep temperature gradient near the surface and where the bottom boundary conditions begin to affect it, due to boundary effects on the velocity, seen as the upturn at the bottom in Figures 3–5. This increase in the mass mixing length near the bottom boundary is an artifact of the bottom boundary condition that the upflows are forced toward vertical at the boundary.

For the Sun, we have deep (20 Mm) convection calculations (Stein et al. 2009) and the determination of the mass mixing length is straightforward, as shown in Figure 3. From this we see that the ratio of the mass mixing length to the pressure scale height, α_P, is nearly constant with depth in the interior of the convection zone until the bottom boundary is approached. The ratio to the density scale height, α_ρ, is closer to 1, but varies slightly and linearly with depth. The mixing length is greater than the density scale height because buoyancy is continually accelerating the flow resulting in a significant contribution from the second term in Equation (3). The mass mixing length in units of the distance to the top of the convection zone, $\alpha _{r-r_{\rm CZ}}$, is also nearly constant with depth in the interior, but has a value of only 0.45 rather than one. In all three cases, the constant proportionality is destroyed at the same point near the bottom boundary due to boundary effects (see Figure 1). At the top of the convection zone, in the highly superadiabatic layer, the mass mixing length becomes large compared to all three measures. Hence, larger values of the mixing length parameters, α_X, are needed near the surface than in the interior, with $\alpha _{r-r_{\rm CZ}}$ becoming largest of all, a factor of five, compared to only a factor of 1.5 for the two other methods. Canuto & Mazzitelli (1991) suggest that a mixing length equal to the distance to the boundary, $\alpha _{r-r_{\rm CZ}}=1$, should work for the entire Sun. This does not match the actual mass mixing in the simulations which is half of the expected value in the interior, adiabatic part of the convection zone, and a factor of five larger near the surface. From a relative logarithmic pressure of 0.5 and down to 4.7 where the bottom boundary effects are felt, this normalization of the mass mixing length results in a concave curve, and it is not obvious which range should be considered close to constant. If we allow larger fluctuations around the mean, that mean value will increase slightly as will the range over which the ratio is within these fluctuations. To obtain a range similar to that of α_P would require accepting fluctuations as large as ±0.27, apart from them not being random fluctuations, but rather a systematic concaveness.

The reasoning behind choosing ℓ = r − r_CZ is based on the MLT concept of convective elements having size ℓ, and that they cannot be larger than their distance from the top of the convection zone. Convection simulations show, however, that the convective topology is one of upward and downward moving streams with continual entrainment of upward moving fluid into the downflows. Convection is a continuous not a discrete process. There is therefore no conceptual problem with an ℓ>r − r_CZ.

As α_P is the most constant measure of the mass mixing length in the interior of the deep solar simulation (Figure 3), and we have similar experiences with other simulations, we advocate a mixing length formulation using α_P. We note that with an adiabatic stratification in hydrostatic equilibrium, H_P = H_ρ/γ₁. So H_P < H_ρ due to the simultaneous variation of ρ and T, and consequently α_P ≃ γ₁α_ρ. In the deep solar simulation γ₁ reaches the completely ionized value of 5/3, but for the other 26 simulations, listed in Table 1, only the six simulations along the high-T_eff/low-g edge get close to fully ionized, and the rest do not get closer than γ₁ ∼ 1.22–1.32 (their minima are in the 1.16–1.22 range). The radial variation of γ₁ therefore also affects the differences between mixing length normalizations with H_ρ or H_P, although the small gradient in the interior makes this a minor effect.

For other stars the available simulations are shallower (about 3 dex less in pressure), so some of the simulations have a smaller range between the surface effects and the bottom boundary effects, and the mixing length for those simulations are, for the moment, less precisely determined.

Table 1 presents the mass mixing length in units of H_P for the 26 stars for which we have performed convection simulations and for which the simulations have a deep enough domain to determine the mass mixing length. The table is ordered in order of increasing log g, and for similar gravity, in order of increasing T_eff. Figures 4 and 5 show the individual α_P determinations and Figure 6 shows how α_P varies as a function of effective temperature and gravity.

Zoom In Zoom Out Reset image size

The basic structure of stellar convection zones is determined by the need to conserve mass. Stars are highly stratified so most of the mass ascending from depth must turn over and go back down within the order of a scale height. The high-density mass from the interior cannot be supported in the low-density outer layers. This produces a true mixing by entrainment of the almost laminar upflows into the turbulent downflows. There is therefore a real physical basis for the MLT of convection. From computational fluid dynamic simulations of stellar convection zones, one can calculate this real mixing of mass to determine the length scale on which it occurs. Naively, one might expect it to be some multiple of the density scale height. However, as we observed for the Sun in Figure 3, it is actually the ratio of the mass mixing length to the pressure scale height which is nearly constant over many orders of magnitude change in the pressure. All of our values of mass mixing lengths are therefore presented in terms of the pressure scale height, H_P.

In Figure 6, we present the variation of α = ℓ/H_P with the atmospheric parameters T_eff and log g. First, comparing the range of α with the uncertainties in Table 1, we see that the variation over the H-R diagram is small, but significant. The range of values found here (1.6–2.2) is consistent with those used in stellar evolution calculations (Salaris & Cassisi 2008; Montalbán & D'Antona 2006; Montalbán et al. 2004; Fernandes et al. 1998; Chieffi et al. 1995; Pedersen et al. 1990), although such calculations assume a constant value for α as the star evolves. We have also compared the mass mixing length values found here with those determined by matching the adiabat in the convection zone from the simulations with that obtained from MLT stellar models (R. Trampedach 2011, in preparation). We find excellent agreement except for the coolest, low mass, dwarf stars where the mass mixing length values are larger. It should also be remembered that in addition to the other mixing length parameters specifying geometry and radiative energy exchange, other input physics, atmospheric opacities, and the T–τ relation also influence the appropriate value of the mixing length to use in stellar structure and evolution calculations.

We find α is smallest toward the high-T_eff and low-g side of the diagram. These are also the simulations that have the least efficient convection, where high speeds and large superadiabatic gradients are necessary for carrying the flux to the photosphere. The least efficient convection has the smallest mass mixing lengths, as expected. The most efficient mixing occurs in the coolest dwarf with the most sedate convection. The solar value from our grid of simulations of α_☉ = 1.76 ± 0.08 agrees well with several determinations from the literature, lending credence to the solar models and to the common practice of calibrating α_MLT against such models. We also see from Figure 6, however, that the solar value is less applicable for lower mass stars, which also have the largest change of α during their evolution. The evolution of higher-mass stars and the ascent up the Hayashi-track seems well described by a single value of α.

The surface layers cannot be properly accounted for with an MLT-like formulation. Crucial components like overshooting, turbulent pressure, radiative transfer in the presence of large temperature inhomogeneities, etc., are not included in this simple model, but greatly affect the surface layers. The mass mixing length found here is the appropriate mixing length to use in a formulation that can realistically include all these phenomena—an ideal MLT formulation. Still lacking such a formulation, the solar mass mixing length found in the present work does not necessarily agree with that found from matching solar evolution models to the Sun's present luminosity and radius. The latter furthermore depends on all the other ingredients going into the atmosphere of the evolution model, such as atmospheric opacities and how the photospheric transition from optically thick to optically thin is treated, e.g., through some match to an atmosphere model or through T–τ relations. The mixing length of MLT is strongly connected to the actual length scale of mixing, as found in this paper, but this connection is weakened by all these confounding factors entering into the MLT version.

The increase in the ratio of mass mixing length to pressure scale height in the superadiabatic layers near the surface raises the question of why stellar evolution calculations with the interior values of α work well in matching observations of global stellar properties. In this regard, note that to match line profiles (e.g., Balmer lines) a small value of α that produces a steep temperature gradient is needed in the atmosphere. In this case, a large value of α must be used in the interior (Montalbán et al. 2004). Our result suggests such a two-value set of mixing length parameters may indeed be appropriate. However, the main conclusion should be that MLT is not well defined and different variations are possible that will yield similar results.

Calculations, such as those presented here, will reduce the uncertainty in the appropriate choice of mixing length in stellar structure and evolution calculations. In particular, we have used the behavior with depth to discriminate between various mixing length formulations, and find that the mass mixing length, ℓ, is best described with a constant α_P = ℓ/H_P. In a future paper (R. Trampedach et al. 2011, in preparation), we present calibrations of α_P(MLT) against the simulations, based on one-dimensional MLT envelope models that employ the exact same EOS and opacities as the simulations and use a new self-consistent treatment of T–τ relations, computed from the simulations (Trampedach et al. 2011). That calibration will be aimed at reproducing the asymptotic adiabat with a standard MLT formulation, and will contain our recommended values for α_P(MLT). Our present calculation of α is model independent in that we compute the actual physical mixing length of the simulations, without connecting it to a particular MLT formulation. Our present result, that α_P seems convergent with depth, is a confirmation that MLT is a physically reasonable model for stellar convection, despite its many shortcomings, particularly near the surface. The present work also confirms the soundness of the scheme employed in the R. Trampedach et al. (2011, in preparation) α_P(MLT) calibration, since, for most of the simulations, it results in α-values very similar to what we find here. This agreement between the two independent methods suggests that α_P(MLT) has been approximately disentangled from the various uncertainties in atmospheric physics that have historically been absorbed into α_MLT.

With the results from Trampedach et al. (2011) and R. Trampedach et al. (2011, in preparation), supported by our present results, one-dimensional MLT models can be made more self-consistent and with a reduced set of free parameters. From an MLT point of view, there are of course still the undetermined geometric form factors that could conceivably be adjusted to match the adiabat in conjunction with the mass mixing length found here.

The calculations reported here were performed at the Australian Partnership for Advanced Computations (APAC) and at the NASA Advanced Supercomputing (NAS) Division of the Ames Research Center. R.T. was supported by the Australian Research Council (grants DP 0342613 and DP 0558836) and NASA grant NNX08AI57G. R.F.S. was supported by NASA grants NNX07AO71G, NNX07AH79G and NNX08AH44G, and NSF grant AST0605738. This support is greatly appreciated.