A Model of Rotating Convection in Stellar and Planetary Interiors. I. Convective Penetration

Boussinesq devised one of the simplest convective systems, consisting of an infinite layer of a nearly incompressible fluid with a small thermal expansion coefficient ${\alpha }_{T}=-\partial \mathrm{ln}\rho /\partial T{| }_{P}$ that is confined between two infinite impenetrable plates differing in temperature by ${\rm{\Delta }}T=T({z}_{2})-T({z}_{c})$, where it is assumed for this model that $T({z}_{2})\lt T({z}_{c})$, and separated by a distance ${{\ell }}_{0}={z}_{2}-{z}_{c}$, as in Figure 1. The thermodynamic variables are further expanded about their averages as $q=\langle q\rangle +\bar{q}(z)+q^{\prime}$, with $\langle q\rangle$ being the volumetric mean, $\bar{q}$ being the horizontal average, and $q^{\prime}$ being the dynamical perturbation. Following Spiegel & Veronis (1960) and Chandrasekhar (1961), the dynamics of a rotating Boussinesq system in Cartesian coordinates may be ascertained from a parametric expansion of the Navier–Stokes equation under the assumption that quantities are scaled relative to and expanded in the parameter ${\alpha }_{T}{\rm{\Delta }}T$. This is equivalent to requiring that the maximum density fluctuation be small in the sense that ${\alpha }_{T}{\rm{\Delta }}T\ll 1$. Under these assumptions, the equations of motion become

$$\begin{eqnarray}&&\left({\partial }_{t}-\nu {{\rm{\nabla }}}^{2}\right){\boldsymbol{v}}=-{\boldsymbol{v}}\,\cdot \,{\boldsymbol{\nabla }}{\boldsymbol{v}},-\displaystyle \frac{1}{\langle \rho \rangle }{\boldsymbol{\nabla }}P^{\prime} +g{\alpha }_{T}T^{\prime} \hat{{\boldsymbol{z}}}-2{\boldsymbol{\Omega }}\times {\boldsymbol{v}},\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&\left({\partial }_{t}-\kappa {{\rm{\nabla }}}^{2}\right)T^{\prime} =\beta {v}_{z}-{\boldsymbol{v}}\,\cdot \,{\boldsymbol{\nabla }}T^{\prime} ,\end{eqnarray} \tag{ 2 }$$

where t is the time coordinate and ${\partial }_{t}$ is the partial derivative with respect to it, ${\boldsymbol{\nabla }}$ is the spatial coordinate gradient. To first order in ${\alpha }_{T}{\rm{\Delta }}T$, the continuity equation requires that the velocity field ${\boldsymbol{v}}$ be solenoidal (${\boldsymbol{\nabla }}\,\cdot \,{\boldsymbol{v}}=0$). The vertical direction $\hat{{\boldsymbol{z}}}$ is anti-aligned with the local direction of the constant effective gravity g, i.e., $\hat{{\boldsymbol{z}}}=-{\boldsymbol{g}}/g$, so that the fluid with $T^{\prime} \gt 0$ rises when ${\rm{\Delta }}T\lt 0$. As the Cartesian f-plane domain assumed here is meant to correspond to a small, local region in the global spherical geometry, the x direction corresponds to the azimuthal direction and the y direction to the latitudinal direction (see Figure 1). Here, $\langle \rho \rangle$ is the mean density, $P^{\prime}$ is the pressure perturbation, $T^{\prime}$ is the temperature perturbation, ${v}_{z}$ is the vertical component of the velocity field, ${\boldsymbol{\Omega }}$ is the local rotational vector that is taken to be in the y–z plane, and κ and ν are the thermal and viscous diffusivities. The vertically aligned thermal gradient to first order in the expansion parameter ${\alpha }_{T}{\rm{\Delta }}T$ is defined as $\beta =d\bar{T}/{dz}+g/{c}_{P}$, where c_P is the specific heat capacity at constant pressure. The adiabatic lapse rate g/c_p arises from the first-order expansion of the ${PdV}$ work in the thermal energy equation.

Zoom In Zoom Out Reset image size

To construct a dispersion relationship for later use in the heat-flux maximization turbulence model, one may take the z component of both the curl and the double curl of Equation (1). Dropping the nonlinear terms, these operations yield the following set of linearized Boussinesq equations:

$$\begin{eqnarray}&&\left({\partial }_{t}-\kappa {{\rm{\nabla }}}^{2}\right)T^{\prime} =\beta {v}_{z},\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&\left({\partial }_{t}-\nu {{\rm{\nabla }}}^{2}\right)\zeta =2{\boldsymbol{\Omega }}\,\cdot \,{\boldsymbol{\nabla }}{v}_{z},\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&\left({\partial }_{t}-\nu {{\rm{\nabla }}}^{2}\right){{\rm{\nabla }}}^{2}{v}_{z}=g{\alpha }_{T}{{\rm{\nabla }}}_{\perp }^{2}T^{\prime} -2{\boldsymbol{\Omega }}\,\cdot \,{\boldsymbol{\nabla }}\zeta ,\end{eqnarray} \tag{ 5 }$$

with ζ being the vertical component of the vorticity and ${{\rm{\nabla }}}_{\perp }$ being the gradient transverse to the vertical direction as in Equations (79), (84), and (85) of Chandrasekhar (1961). As seen in many papers regarding Boussinesq dynamics (e.g., Chandrasekhar 1961), this set of equations can be reduced to a single third order in time and eighth order in space equation for the vertical velocity as

$$\begin{eqnarray}&&\begin{array}{l}\left({\partial }_{t}-\kappa {{\rm{\nabla }}}^{2}\right){\left({\partial }_{t}-\nu {{\rm{\nabla }}}^{2}\right)}^{2}{{\rm{\nabla }}}^{2}{v}_{z}+g{\alpha }_{T}\beta {{\rm{\nabla }}}_{\perp }^{2}\left({\partial }_{t}-\nu {{\rm{\nabla }}}^{2}\right){v}_{z}\\ \qquad +4{\boldsymbol{\Omega }}\,\cdot \,{\boldsymbol{\nabla }}\left[{\boldsymbol{\Omega }}\,\cdot \,{\boldsymbol{\nabla }}\left({\partial }_{t}-\kappa {{\rm{\nabla }}}^{2}\right){v}_{z}\right]=0.\end{array}\end{eqnarray} \tag{ 6 }$$

For impenetrable and stress-free boundary conditions, the solutions of Equation (6) are periodic in the horizontal, sinusoidal in the vertical, and exponential in time, i.e., ${v}_{z}=v\,\sin \left[{k}_{z}\left(z-{z}_{c}\right)\right]\exp (i{{\boldsymbol{k}}}_{\perp }\,\cdot \,{\boldsymbol{r}}+{st})$, where ${{\boldsymbol{k}}}_{\perp }$ is the horizontal wavevector, s is the growth rate, ${\boldsymbol{r}}$ is the local coordinate vector, and $v$ is a to-be-determined constant velocity amplitude. To satisfy the impenetrable, stress-free, and fixed temperature boundary conditions, it is required that the vertical wavenumber be ${k}_{z}=n\pi /{{\ell }}_{0}$. This can also be considered to be equivalent to modeling only the portion of the spectrum of convective eddies or motions whose Mach numbers are small and those that are much smaller than the integral scale of the system. The introduction of this solution into Equation (6) yields the following dispersion relationship that relates s to the wavevector ${\boldsymbol{k}}$ as

$$\begin{eqnarray}&&\begin{array}{l}\left[\left(s+\kappa {k}^{2}\right){\left(s+\nu {k}^{2}\right)}^{2}{k}^{2}+g{\alpha }_{T}\beta {k}_{\perp }^{2}\left(s+\nu {k}^{2}\right)\right.\\ \qquad \left.+4{\left({\boldsymbol{\Omega }}\cdot {\boldsymbol{k}}\right)}^{2}\left(s+\kappa {k}^{2}\right)\right]{v}_{z}=0.\end{array}\end{eqnarray} \tag{ 7 }$$

This equation can be simplified and made nondimensional by dividing through by the appropriate powers of N_* and k_z, where ${N}_{* }^{2}=| g{\alpha }_{T}\beta |$ is the absolute value of the square of the Brunt–Väisälä frequency as in Barker et al. (2014), which is otherwise negative in a convective region. Utilizing the following auxiliary quantities,

$$\begin{eqnarray}\begin{array}{rcl}\hat{s} & = & \displaystyle \frac{s}{{N}_{* }},\,{z}^{3}=1+{a}^{2}=\displaystyle \frac{{k}^{2}}{{k}_{z}^{2}},\,{a}^{2}=\displaystyle \frac{{k}_{x}^{2}}{{k}_{z}^{2}}+\displaystyle \frac{{k}_{y}^{2}}{{k}_{z}^{2}}={a}_{x}^{2}+{a}_{y}^{2},\\ K & = & \displaystyle \frac{\kappa {k}_{z}^{2}}{{N}_{* }},\quad V=\displaystyle \frac{\nu {k}_{z}^{2}}{{N}_{* }},\end{array}\end{eqnarray} \tag{ 8 }$$

the dispersion relationship becomes

$$\begin{eqnarray}&&\begin{array}{l}\left(\hat{s}+{{Kz}}^{3}\right)\left({z}^{3}{\left(\hat{s}+{{Vz}}^{3}\right)}^{2}+{O}^{2}{\left(\cos \theta +{a}_{y}\sin \theta \right)}^{2}\right)\\ \qquad -\,\left({z}^{3}-1\right)\left(\hat{s}+{{Vz}}^{3}\right)=0,\end{array}\end{eqnarray} \tag{ 9 }$$

where $4{({\boldsymbol{\Omega }}\cdot {\boldsymbol{k}})}^{2}/{N}_{* }^{2}={k}_{z}^{2}{O}^{2}{\left(\cos \theta +{a}_{y}\sin \theta \right)}^{2}$ with

$$\begin{eqnarray}&&{O}^{2}=\displaystyle \frac{4{{\rm{\Omega }}}_{0}^{2}}{{N}_{* }^{2}},\end{eqnarray} \tag{ 10 }$$

where ${{\rm{\Omega }}}_{0}$ is the bulk rotation rate of the system.

In assessing the behavior of the Boussinesq heat flux, both with and without diffusion, there are several limits of the diffusive parameters that can be considered, each with its own physical justification. In highly turbulent regimes, it is often assumed that both diffusivities can be ignored, which is the case considered in S79, and expanded upon below in Section 2.4. Within stellar interiors, the molecular value of the thermal Prandtl number ($\Pr =\nu /\kappa$) is typically very small, being of the order of 10⁻⁵ or smaller, so the limit $\nu \to 0$ may be considered as in Section 2.5. The case where thermal diffusion can be ignored ($\kappa \to 0$) can be relevant to very high Prandtl number systems such as those found in geophysical contexts, but it is not directly considered here. However, since the general case with both diffusivities can be treated here (Section 2.6), turbulent diffusivities may also be incorporated into the model, where typically it is assumed that $\Pr$ is approximately unity. Moreover, to make contact with numerical simulations and laboratory experiments, it is useful to retain diffusive effects.

The secular impacts of rotation and magnetic fields on stellar and planetary evolution are of keen interest within the astrophysical community (e.g., Maeder 2009; Mathis 2013). Accordingly, extensions to mixing length theory (MLT) and its ilk have been introduced (e.g., Stevenson 1979; Canuto & Hartke 1986; Canuto & Mazzitelli 1991; Zhou 1995; Xiong et al. 1997). As expounded upon in Stevenson (1979, hereafter S79), a surprisingly effective approach to including rotation in MLT is to hypothesize a convection model where the convective length scale, degree of superadiabaticity, and velocity are governed by the linear mode that maximizes the convective heat flux. This heuristic method was made more rigorous when incorporated into the spectrally nonlocal large-scale turbulence model of Canuto & Hartke (1986), where it is used to set the scale of energy injection of a Heisenberg–Kolmogorov spectrum of turbulent energy cascade (Kolmogorov 1941; Heisenberg 1948). Subsequently, it sets the local values of the turbulent diffusivities, although these theories have not yet widely been employed in stellar evolution computations, with one exception being Ireland & Browning (2018). One can also incorporate approximate 2D dynamics through models of the Reynolds and Maxwell-stresses, which has been significantly developed over the last several decades (e.g., Durney & Spruit 1979; Hathaway 1984; Ruediger 1989; Garaud et al. 2010).

The S79 model of rotating convection has its origins in the principle of maximum heat transport proposed by Malkus (1954). In that principle, an upper limit for a boundary-condition-dependent turbulent heat flux that depends upon the smallest Rayleigh-unstable convective eddy is established. The size of this eddy is determined with a variational technique that is similar to that developed in Chandrasekhar (1961) for the determination of the Rayleigh number, which is the ratio of the buoyancy force to the viscous force multiplied by the ratio of the thermal to viscous diffusion timescales. This technique then permits the independent computation of the rms values of the fluctuating temperature and velocity amplitudes. However, the heat-flux maximization principle is not fully rigorous as there is no a priori theoretical justification as to why the turbulence would arrange itself so as to maximize the heat flux (e.g., Howard 1963). It is instead built upon three assumptions: first, the mean-temperature gradient must be everywhere negative; second, there is a finite range of wavenumbers that are effective at transporting heat; and third, that the highest vertical wavenumber contributing to the heat transport is marginally stable with respect to a given mean-temperature profile. As Howard (1963) has shown, however, the first assumption can be relaxed if the power integrals of the Boussinesq equations are considered in place of the Boussinesq equations. It was also shown that there will be a hierarchy of such constraint integrals that close the theory for each order of expansion. So, it may be more practical to appeal to a more mathematically complete means to compute the most probable equilibrium states, which is the variational technique established in Prigogine & Glansdorff (1965). Indeed, as Spiegel (1962), Townsend (1962), and Howard (1963) all point out, these optimal solutions are formed from solutions to the linear equations, and hence they are not true solutions to the full nonlinear Boussinesq equations. However, this flaw can be mollified if, as in S79, it is conjectured that the convective modes are amplitude-limited by instabilities in a nonlinearly saturated state.

The Rayleigh–Bénard experiments carried out in Townsend (1959) aimed to directly test the heat-flux maximization principle of Malkus (1954), where the temperature and its gradient are measured to compute the heat flux and the moments of the temperature distribution. These values are then compared against the mean-temperature gradient predicted in Malkus (1954). Townsend (1959) found that indeed, the measured mean-temperature gradient is well described by those expected from the heat-flux maximization principle, being inversely proportional to the vertical coordinate outside of the boundary layer. This is in contrast to the self-similarity solution of Priestley (1954), which predicts a scaling proportional to the inverse cube root of the vertical coordinate. Howard (1963) reexamined those experimental data with an emphasis on testing whether or not both the mean fields and the fluctuating fields predicted in Malkus (1954) correspond to the measured flow structures; they are both found to match the theoretical values to within a factor of order unity. Additionally, Howard (1963) constructed a variational technique to assess the Rayleigh number ($\mathrm{Ra}=-g{\alpha }_{T}{\rm{\Delta }}T{{\ell }}_{0}^{3}/\nu \kappa$) and the Nusselt number (which measures the ratio of the total heat flux to the rate of thermal diffusion) of the heat-flux maximizing flows. It was also shown that, when the continuity equation is imposed in addition to the power integrals, the experimentally measured Nusselt number is asymptotically well approximated by the variationally obtained value, although its magnitude is overestimated by about a factor of 2. This approach can also be taken to compare Malkus's theory to more recent laboratory and numerical experiments of the Rayleigh–Bénard convection (e.g., Funfschilling et al. 2005; Ahlers et al. 2006; Anders & Brown 2017). Furthermore, similar experiments of the Rayleigh–Bénard convection in rotating cylinders have measured the Nusselt number scaling with respect to the Rayleigh number and with the Ekman number $[\mathrm{Ek}=\nu /(2{\rm{\Omega }}{{\ell }}_{0}^{2})]$ in the rotating cases (King et al. 2009, 2012, 2013; Zhong et al. 2009; King & Aurnou 2012; Aurnou et al. 2015; Cheng et al. 2015; Weiss et al. 2016). Thus, a Malkus-like turbulence theory that includes rotation, as in S79 and as shown here, can in principle also be compared to those experiments, which will be a focus of later work.

Laboratory experiments (Coates & Ivey 1997; Weiss et al. 2016) and numerical simulations have lent some credence to this convection model (Käpylä et al. 2005; Barker et al. 2014; Sondak et al. 2015). In particular, those simulations indicate that the low convective Rossby number scaling regime established in S79 appears to hold up well for three decades in convective Rossby number and for about one decade in Nusselt number. Moreover, a recent implementation of Stevenson's asymptotic scaling regimes in convective Rossby number has been used to provide a local model of a rotationally dependent mixing-length theory (Ireland & Browning 2018). There, it was shown that the entropy gradient of fully convective stars can be significantly modified by rapid rotation in the bulk of their interiors, excluding the surface region where the convective Rossby number is large. Therefore, the adiabat of the star is modified, leading to structural and evolutionary changes. What remains to be shown is how such a model of convection can impact the mixing for more modest rotation rates where the convective Rossby numbers are closer to unity, the depth of convective penetration, as well as the amplitude of wave-driven and shear-induced transport mechanisms.

Following S79, the nonlinear saturation conjecture used in conjunction with the maximum heat transport principle requires that the primary contribution to each scale of the temperature and velocity fluctuation is from the convective term at the corresponding scale. This implies that its growth rate is $s=\left|{\boldsymbol{v}}\right|\left|{\boldsymbol{k}}\right|$, which also provides a way to estimate the heat flux $F=\rho {c}_{P}{\mathfrak{R}}[{v}_{z}T^{\prime} ]$ with linear dynamics. Moreover, in the theory of Malkus and Howard, the mode that carries the largest fraction of the heat flux is the fundamental vertical mode, i.e., n = 1, so that ${k}_{z}=\pi /{{\ell }}_{0}$, but whose horizontal structure is undetermined. It is this mode that Stevenson's turbulence model is built around. Nominally, however, the maximization should be carried out scalewise for each of the modes that contribute to the heat flux. This task will be left for later work that examines non-Boussinesq turbulence using Malkus's theory as a point of comparison. Thus, Equation (1), with the assumptions used above to formulate the dispersion relationship, yields

$$\begin{eqnarray}&&s{\boldsymbol{v}}=-i{\boldsymbol{k}}P^{\prime} /\langle \rho \rangle +g{\alpha }_{T}T^{\prime} \hat{{\boldsymbol{z}}}-2{\boldsymbol{\Omega }}\times {\boldsymbol{v}}-\nu {k}^{2}{\boldsymbol{v}},\end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}&&i{\boldsymbol{k}}\,\cdot \,{\boldsymbol{v}}=0.\end{eqnarray} \tag{ 12 }$$

Taking the dot product of ${\boldsymbol{v}}$ and Equation (11), using Equation (12) to eliminate the pressure, and taking a volume average, it can be seen that

$$\begin{eqnarray}&&\left(s+\nu {k}^{2}\right)\langle {v}^{2}\rangle =g{\alpha }_{T}\langle T^{\prime} {v}_{z}\rangle .\end{eqnarray} \tag{ 13 }$$

Therefore,

$$\begin{eqnarray}&&F=\displaystyle \frac{\rho {c}_{P}{\mathfrak{R}}[\left(s+\nu {k}^{2}\right)\langle {v}^{2}\rangle ]}{g{\alpha }_{T}}.\end{eqnarray} \tag{ 14 }$$

Employing the hypothesis that the velocity in the quasi-stationary turbulent state scales as $\langle {v}^{2}\rangle ={{ss}}^{* }/{k}^{2}$, one has the heat flux

$$\begin{eqnarray}&&F=\displaystyle \frac{\langle \rho \rangle {c}_{P}}{g{\alpha }_{T}}\displaystyle \frac{{\mathfrak{R}}[\left(s+\nu {k}^{2}\right)]{{ss}}^{* }}{{k}^{2}}\end{eqnarray} \tag{ 15 }$$

$$\begin{eqnarray}&&=\,\displaystyle \frac{\langle \rho \rangle {c}_{P}{N}_{* }^{3}}{g{\alpha }_{T}{k}_{z}^{2}}\left[\displaystyle \frac{{\mathfrak{R}}[\hat{s}]}{{z}^{3}}+V\right]{\left|\hat{s}\right|}^{2},\end{eqnarray} \tag{ 16 }$$

where the definitions given in Equations (8) have been used to simplify the latter expression.

It is useful to eliminate some of the possible solutions to the maximization of this heat flux. First, consider the nondiffusive case, where

$$\begin{eqnarray}&&{\hat{s}}^{2}=\displaystyle \frac{{a}^{2}-{O}^{2}{\left(\cos \theta +\sin \theta {a}_{y}\right)}^{2}}{1+{a}^{2}},\end{eqnarray} \tag{ 17 }$$

from which it is clear that any value of a_y will reduce the growth rate. This also places a lower limit on the value of a for which the solutions are real with a growing branch, with ${a}_{x}\gt O\cos \theta$, which is shown to be satisfied by the heat-flux maximization in Stevenson (1979) and in the next subsection. Likewise, in the nonrotating but diffusive case, one has

$$\begin{eqnarray}&&\left[\hat{s}+{{Vz}}^{3}\right]\left[{z}^{3}\left(\hat{s}+{{Kz}}^{3}\right)\left(\hat{s}+{{Vz}}^{3}\right)-{a}^{2}\right]=0.\end{eqnarray} \tag{ 18 }$$

This dispersion relationship has a growing real solution for sufficiently small V and K as

$$\begin{eqnarray}\begin{array}{rcl}\hat{s} & = & {z}^{-3}\sqrt{{\left(V+K\right)}^{2}{z}^{12}+4{z}^{3}\left({a}^{2}-{{KVz}}^{9}\right)}\\ & & -\,\displaystyle \frac{1}{2}\left(V+K\right){z}^{3}.\end{array}\end{eqnarray} \tag{ 19 }$$

When linearized with respect to V and K, this yields

$$\begin{eqnarray}&&\hat{s}={{az}}^{-\tfrac{3}{2}}-\displaystyle \frac{1}{2}\left(V+K\right){z}^{3},\end{eqnarray} \tag{ 20 }$$

which is maximized when a_y = 0. So, in all instances of parameter regimes considered below, the maximization will be carried out over ${z}^{3}=1+{a}_{x}^{2}$ with a_y = 0. Note that this is equivalent to considering only the two-dimensional rolls aligned with the rotation axis, which in the rotating case is due largely to the Taylor–Proudman constraint. These modes also have a larger growth rate than three-dimensional modes for nonrotating Rayleigh–Bénard convection, but nevertheless, the 2D modes can approximate many of the behaviors of the 3D modes (e.g., van der Poel et al. 2013).

It is also useful to recast some of the parameters so that the system can be more readily compared to the S79 results. The characteristic velocity ${v}_{0}$ is derived from the growth rate and maximizing wavevector in the nonrotating and nondiffusive case, which are ${s}_{0}^{2}=3/5| {g}_{0}{\alpha }_{T}{\beta }_{0}| =3/5{N}_{* ,0}^{2}$, with ${\beta }_{0}$ being the thermal gradient and g₀ being the effective gravity in the nonrotating case, and ${k}_{0}^{2}=5/2{k}_{z}^{2}$ as in Equation (36) in Stevenson (1979), leading to

$$\begin{eqnarray}&&{v}_{0}=\displaystyle \frac{{s}_{0}}{{k}_{0}}=\displaystyle \frac{\sqrt{6}}{5}\displaystyle \frac{{N}_{* ,0}}{{k}_{z}}=\displaystyle \frac{\sqrt{6}}{5\pi }{N}_{* ,0}{{\ell }}_{0}.\end{eqnarray} \tag{ 21 }$$

This permits the definition of the convective Rossby number ${\mathrm{Ro}}_{{\rm{c}}}$, which for this Boussinesq system is

$$\begin{eqnarray}&&{\mathrm{Ro}}_{{\rm{c}}}=\displaystyle \frac{{v}_{0}}{2{{\rm{\Omega }}}_{0}{{\ell }}_{0}}=\displaystyle \frac{\sqrt{6}{N}_{* ,0}}{10\pi {{\rm{\Omega }}}_{0}},\end{eqnarray} \tag{ 22 }$$

which implies that

$$\begin{eqnarray}&&O=\displaystyle \frac{2{{\rm{\Omega }}}_{0}}{{N}_{* }}=\displaystyle \frac{{v}_{0}}{{N}_{* }{\mathrm{Ro}}_{{\rm{c}}}{{\ell }}_{0}}=\displaystyle \frac{\sqrt{6}{N}_{* ,0}}{5\pi {N}_{* }{\mathrm{Ro}}_{{\rm{c}}}}.\end{eqnarray} \tag{ 23 }$$

Next consider the variation of the superadiabaticity, which for this system is given by $\epsilon ={H}_{P}\beta /T$, meaning that ${N}_{* }^{2}=| g{\alpha }_{T}T\epsilon /{H}_{P}|$, where H_P is the pressure scale height. The comparison made throughout this paper is between cases with additional included physical effects and the base case that is nondiffusive and nonrotating. The potential temperature gradient in this latter case is ascertained from the Malkus–Howard turbulence model, which yields a value of ${N}_{* ,0}$. So far, all quantities have been normalized with respect to N_*. Instead, it is useful to compare them relative to ${N}_{* ,0}$ and to introduce the ratio of superadiabaticities as an additional unknown as

$$\begin{eqnarray}&&q={N}_{* ,0}/{N}_{* }.\end{eqnarray} \tag{ 24 }$$

Therefore, all parametric quantities have the following equivalences:

$$\begin{eqnarray}\begin{array}{rcl}O & = & q\displaystyle \frac{\sqrt{6}}{5\pi {\mathrm{Ro}}_{{\rm{c}}}}={{qO}}_{0},\\ K & = & q\displaystyle \frac{\kappa {k}_{z}^{2}}{{N}_{* ,0}}={{qK}}_{0},\\ V & = & q\displaystyle \frac{\nu {k}_{z}^{2}}{{N}_{* ,0}}={{qV}}_{0}.\end{array}\end{eqnarray} \tag{ 25 }$$

So, the dispersion relationship (Equation (9)) and the heat flux (Equation (16)) may be written as

$$\begin{eqnarray}&&\begin{array}{l}\left(\hat{s}+{K}_{0}{{qz}}^{3}\right)\left({z}^{3}{\left(\hat{s}+{V}_{0}{{qz}}^{3}\right)}^{2}+{O}_{0}^{2}{q}^{2}{\cos }^{2}\theta \right)\\ \qquad -\left({z}^{3}-1\right)\left(\hat{s}+{V}_{0}{{qz}}^{3}\right)=0,\end{array}\end{eqnarray} \tag{ 26 }$$

$$\begin{eqnarray}&&F=\displaystyle \frac{{F}_{0}}{{q}^{3}}\left[\displaystyle \frac{{\hat{s}}^{3}}{{z}^{3}}+{V}_{0}q{\hat{s}}^{2}\right],\end{eqnarray} \tag{ 27 }$$

where ${F}_{0}=\langle \rho \rangle {c}_{P}{N}_{* ,0}^{3}/\left(g{\alpha }_{T}{k}_{z}^{2}\right)$. When maximized for a positive real value of $\hat{s}$, Equation (27) implies that

$$\begin{eqnarray}&&\displaystyle \frac{d\hat{s}}{{dz}}=\displaystyle \frac{3{\hat{s}}^{2}}{z\left(3\hat{s}+2{V}_{0}{{qz}}^{3}\right)}.\end{eqnarray} \tag{ 28 }$$

Depending upon the method of finding $\hat{s}$, it can be quite complicated to compute directly as a function of z. Instead, it is useful to compute the implicit derivative of the dispersion relationship (Equation (26)) to obtain $d\hat{s}/{dz}$, which yields

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{d\hat{s}}{{dz}}=-3{z}^{2}\left[{V}_{0}q\left(1-2{z}^{3}\right)+{K}_{0}{q}^{3}\left({O}_{0}^{2}{\cos }^{2}\theta +4{V}_{0}^{2}{z}^{9}\right)\right.\\ \left.+\hat{s}\left({\hat{s}}^{2}+2q\left({K}_{0}+2{V}_{0}\right){z}^{3}\hat{s}+3{V}_{0}{q}^{2}\left(2{K}_{0}+{V}_{0}\right){z}^{6}-1\right)\right]\\ \left[{z}^{3}\left(3{\hat{s}}^{2}+2q\left({K}_{0}+2{V}_{0}\right){z}^{3}\hat{s}+{V}_{0}{q}^{2}\left(2{K}_{0}+{V}_{0}\right){z}^{6}-1\right)\right.\\ {\left.+{q}^{2}{O}_{0}^{2}{\cos }^{2}\theta +1\right]}^{-1}.\end{array}\end{eqnarray} \tag{ 29 }$$

Finally, to assess the scaling of the superadiabaticity, the velocity, and the horizontal wavevector, a further assumption must be made in which the maximum heat flux is invariant to any parameters, namely that $\max \left[F\right]=\max {\left[F\right]}_{0}$, so the heat flux is equal to the maximum value $\max {\left[F\right]}_{0}$ obtained in the Malkus–Howard turbulence model for the nonrotating case. This is based on the assumption that the total heat flux should not change with rotation, which is not true in general as the centripetal acceleration of the star can lead to a variation in the central temperature and density of the star, and thus to a variation of the total luminosity of the star at a fixed mass (e.g., Rieutord et al. 2016). Furthermore, the flux will be latitude-dependent due to the rotational variation of the local gravity (e.g., Maeder 1999; Wang et al. 2016). Indeed, an analysis of 3D nonlinear global-scale convection simulations of rapidly rotating oblate stars has been carried out in Wang et al. (2016). These simulations indicate that the Von Zeipel theorem (von Zeipel 1924) holds well even for the emergent flux at the outer boundary, with the exception of a region near the equator where there is a flux enhancement. Thus, while acknowledging that the flux derived from the Von Zeipel theorem is not rigorously representative of the heat flux in a convection zone, it provides a very good approximation. Following Maeder (1999) and the von Zeipel theorem, one can show that the heat flux F_V on an isobar with $P=\langle P\rangle$ varies as

$$\begin{eqnarray}&&{F}_{V}=\displaystyle \frac{L\left(\langle P\rangle \right)}{4\pi {{GM}}_{* }\left(\langle P\rangle \right)}g,\end{eqnarray} \tag{ 30 }$$

where $L\left(\langle P\rangle \right)$ is the total luminosity on the isobar and ${M}_{* }\left(\langle P\rangle \right)$ is the total mass integrated to the isobar. Approximating the effective gravity and the shift in the mass as in Maeder (1999, Equation (2.3)), this becomes

$$\begin{eqnarray}&&{F}_{V}\approx \displaystyle \frac{L\left(\langle P\rangle \right)\left({g}_{0}-{{\rm{\Omega }}}_{0}^{2}r\sin \theta \right)}{4\pi {GM}\left[1-{{\rm{\Omega }}}_{0}^{2}/\left(2\pi G\langle \rho \rangle \right)\right]},\end{eqnarray} \tag{ 31 }$$

where r and θ correspond to the origin of the local coordinate system as in Figure 1. Taking the ratio of the rotating and nonrotating fluxes, one has

$$\begin{eqnarray}&&\displaystyle \frac{{F}_{V}}{{F}_{V,0}}\approx \displaystyle \frac{L\left(\langle P\rangle \right)}{L\left(r\right)}\displaystyle \frac{1-{{\rm{\Omega }}}_{0}^{2}r\sin \theta /{g}_{0}}{1-{{\rm{\Omega }}}_{0}^{2}/\left(2\pi G\langle \rho \rangle \right)},\end{eqnarray} \tag{ 32 }$$

which in terms of the convective Rossby number and the buoyancy frequency is

$$\begin{eqnarray}&&\displaystyle \frac{{F}_{V}}{{F}_{V,0}}\approx \displaystyle \frac{L\left(\langle P\rangle \right)}{L\left(r\right)}\displaystyle \frac{1-3{N}_{* ,0}^{2}r\sin \theta /\left(50{\pi }^{2}{\mathrm{Ro}}_{{\rm{c}}}^{2}{g}_{0}\right)}{1-3{N}_{* ,0}^{2}/\left(100{\pi }^{3}{\mathrm{Ro}}_{{\rm{c}}}^{2}G\langle \rho \rangle \right)}.\end{eqnarray} \tag{ 33 }$$

Noting the definition of the buoyancy frequency amplitude and the local density, one can see that

$$\begin{eqnarray}&&\displaystyle \frac{{F}_{V}}{{F}_{V,0}}\approx \displaystyle \frac{L\left(\langle P\rangle \right)}{L\left(r\right)}\displaystyle \frac{1-3| {\alpha }_{T}{\rm{\Delta }}T| \,\sin \theta /\left(50{\pi }^{2}{\mathrm{Ro}}_{{\rm{c}}}^{2}\right)}{1-| {\alpha }_{T}{\rm{\Delta }}T| {g}_{0}{r}^{2}/\left(25{\pi }^{2}{\mathrm{Ro}}_{{\rm{c}}}^{2}{GM}\right)}\end{eqnarray} \tag{ 34 }$$

$$\begin{eqnarray}&&\approx \displaystyle \frac{L\left(\langle P\rangle \right)}{L\left(r\right)}\left[1+\displaystyle \frac{| {\alpha }_{T}{\rm{\Delta }}T| }{25{\pi }^{2}{\mathrm{Ro}}_{{\rm{c}}}^{2}}\left(\displaystyle \frac{{g}_{0}{r}^{2}}{{GM}}-\displaystyle \frac{3}{2}\sin \theta \right)\right]+\cdots .\end{eqnarray} \tag{ 35 }$$

Therefore, because the linearized change in flux is proportional to the Boussinesq expansion parameter $| {\alpha }_{T}{\rm{\Delta }}T|$ and the inverse square of the convective Rossby number, this requires ${\mathrm{Ro}}_{{\rm{c}}}\gt {5}^{-1}{\pi }^{-1}| {\alpha }_{T}{\rm{\Delta }}T{| }^{1/2}$ for the approximations of the model to hold. Thus, for now, these effects will be ignored as they are sensitive to the global-scale dynamics, the microphysics of the nuclear energy generation rate, and the equation of state. Therefore, building this convection model consists of three steps: defining a dispersion relationship that links $\hat{s}$ to q and z, maximizing the heat flux with respect to z, and assuming an invariant maximum heat flux that then closes this three-variable system.

In the case of planetary and stellar interiors, the viscous damping timescale is generally longer than the convective-overturning timescale (i.e., ${V}_{0}\ll {N}_{* ,0}$). Thus, the maximized heat-flux invariance is much simpler to treat. In particular, the heat-flux invariance condition under this assumption is then

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{\max \left[F\right]}{\max {\left[F\right]}_{0}} & = & \displaystyle \frac{25}{6}\sqrt{\displaystyle \frac{5}{3}}{\left[\displaystyle \frac{{\hat{s}}^{3}}{{q}^{3}{z}^{3}}+\displaystyle \frac{{V}_{0}{\hat{s}}^{2}}{{q}^{2}}\right]}_{\max }\\ & \approx & {\left.\displaystyle \frac{25}{6}\sqrt{\displaystyle \frac{5}{3}}\displaystyle \frac{{\hat{s}}^{3}}{{q}^{3}{z}^{3}}\right|}_{\max }=1,\end{array}\end{eqnarray} \tag{ 36 }$$

implying that

$$\begin{eqnarray}&&\hat{s}=\tilde{s}{qz}+{ \mathcal O }({V}_{0}),\end{eqnarray} \tag{ 37 }$$

where $\tilde{s}={2}^{1/3}{3}^{1/2}{5}^{-5/6}$ and $\max {\left[F\right]}_{0}=6/25\sqrt{3/5}{F}_{0}$ follows from the definition of the flux and the maximizing wavevector used to define ${v}_{0}$ above in Equation (21).

The assumption of this convection model is that the magnitude of the velocity is defined as the ratio of the maximizing growth rate and wavevector. With the above approximation, the velocity amplitude can be defined relative to the nondiffusive and nonrotating case scales without a loss of generality as

$$\begin{eqnarray}&&\displaystyle \frac{v}{{v}_{0}}=\displaystyle \frac{{k}_{0}}{{s}_{0}}\displaystyle \frac{s}{k}=\displaystyle \frac{5}{\sqrt{6}}\displaystyle \frac{{N}_{* }}{{N}_{* ,0}}\displaystyle \frac{\hat{s}}{{z}^{3/2}}=\displaystyle \frac{5}{\sqrt{6}}\displaystyle \frac{\hat{s}}{{{qz}}^{3/2}}={\left(\displaystyle \frac{5}{2}\right)}^{\tfrac{1}{6}}{z}^{-\tfrac{1}{2}}.\end{eqnarray} \tag{ 38 }$$

So, only the maximizing wavevector needs to be found in order to ascertain the relative velocity amplitude. For reference, the symbols that will be frequently used from this section are listed in Table 1.

Table 1.  Frequently Used Symbols in the Convection Model

$a={k}_{x}/{k}_{z}$	Maximizing horizontal wavevector
${k}_{z}=\pi /{{\ell }}_{0}$	Maximizing vertical wavevector
${K}_{0}=\kappa {k}_{z}^{2}/{N}_{* ,0}$	Normalized thermal diffusivity
${O}_{0}=\sqrt{6}/\left(5\pi {\mathrm{Ro}}_{{\rm{c}}}\right)$	Normalized Coriolis coefficient
${\mathrm{Ro}}_{{\rm{c}}}=\sqrt{6}{N}_{* ,0}/\left(10\pi {{\rm{\Omega }}}_{0}\right)$	Convective Rossby number
$\hat{s}=s/{N}_{* }$	Normalized growth rate
${v}_{0}=\sqrt{6}{N}_{* ,0}{{\ell }}_{0}/\left(5\pi \right)$	Velocity of the nonrotating case
${V}_{0}=\nu {k}_{z}^{2}/{N}_{* ,0}$	Normalized viscosity
$q={N}_{* ,0}/{N}_{* }$	Ratio of buoyancy timescales
${z}^{3}={k}^{2}/{k}_{z}^{2}$	Normalized wavevector

Download table as:  ASCIITypeset image

As considered in S79, the simplest case to consider is one in which all diffusive mechanisms are ignored. Given Equation (28), the maximization condition requires

$$\begin{eqnarray}&&\displaystyle \frac{d\hat{s}}{{dz}}=\displaystyle \frac{\hat{s}}{z},\end{eqnarray} \tag{ 39 }$$

where the derivative of the growth rate can be determined from the nondiffusive limit of Equation (29) as

$$\begin{eqnarray}&&\displaystyle \frac{d\hat{s}}{{dz}}=-\displaystyle \frac{3{z}^{2}\hat{s}\left({\hat{s}}^{2}-1\right)}{1+{O}_{0}^{2}{q}^{2}{\cos }^{2}\theta +{z}^{3}\left(3{\hat{s}}^{2}-1\right)};\end{eqnarray} \tag{ 40 }$$

when equated, this yields the following equation for the growth rate

$$\begin{eqnarray}&&{\hat{s}}^{2}=\displaystyle \frac{4{z}^{3}-{O}_{0}^{2}{q}^{2}{\cos }^{2}\theta -1}{6{z}^{3}}.\end{eqnarray} \tag{ 41 }$$

However, the dispersion relationship (Equation (26)) for the nondiffusive case requires

$$\begin{eqnarray}&&{\hat{s}}^{2}=\displaystyle \frac{{z}^{3}-{O}_{0}^{2}{q}^{2}{\cos }^{2}\theta -1}{{z}^{3}}.\end{eqnarray} \tag{ 42 }$$

Using the heat-flux invariance condition (Equation (37)) so that $\hat{s}=\tilde{s}{qz}$, the two above equations can be manipulated to isolate a single equation for the maximizing value of z and a dependent equation for q as

$$\begin{eqnarray}&&2{z}^{5}-5{z}^{2}-\displaystyle \frac{3}{{\tilde{s}}^{2}}{O}_{0}^{2}{\cos }^{2}\theta =0,\end{eqnarray} \tag{ 43 }$$

$$\begin{eqnarray}&&{q}^{2}-\displaystyle \frac{{z}^{3}-1}{{\tilde{s}}^{2}{z}^{5}+{O}_{0}^{2}{\cos }^{2}\theta }=0.\end{eqnarray} \tag{ 44 }$$

Together with Equations (25), these two equations yield the following scalings with convective Rossby number and maximizing wavevector as

$$\begin{eqnarray}&&\displaystyle \frac{\epsilon }{{\epsilon }_{0}}=\displaystyle \frac{25{\pi }^{2}{\mathrm{Ro}}_{{\rm{c}}}^{2}{\tilde{s}}^{2}{z}^{5}+6{\cos }^{2}\theta }{25{\pi }^{2}{\mathrm{Ro}}_{{\rm{c}}}^{2}\left({z}^{3}-1\right)},\end{eqnarray} \tag{ 45 }$$

$$\begin{eqnarray}&&2{z}^{5}-5{z}^{2}-\displaystyle \frac{18{\cos }^{2}\theta }{25{\pi }^{2}{\mathrm{Ro}}_{{\rm{c}}}^{2}{\tilde{s}}^{2}}=0.\end{eqnarray} \tag{ 46 }$$

With the quintic form of Equation (46), its solutions are not representable as radicals or other functions. So, it is left in its current form to be solved numerically for a given colatitude and convective Rossby number. However, at low convective Rossby number, these equations can be asymptotically approximated: for ${\mathrm{Ro}}_{{\rm{c}}}\ll 1$, it is clear that $z\propto {\mathrm{Ro}}_{{\rm{c}}}^{-2/5}$ because ${z}^{2}\ll {z}^{5}$. It then follows that $\epsilon /{\epsilon }_{0}\propto {\mathrm{Ro}}_{{\rm{c}}}^{-4/5}$ and $v/{v}_{0}\propto {\mathrm{Ro}}_{{\rm{c}}}^{1/5}$, both of which match the scalings given in S79 and the direct numerical simulations of Käpylä et al. (2005) and Barker et al. (2014), and can be seen in Figure 2(a). In contrast, at very large convective Rossby number, these quantities converge to unity as expected.

Zoom In Zoom Out Reset image size

At first glance, one can see from Equations (45) and (46) that at the equator, there is no rotational scaling of any of these quantities. However, there is a rapid transition from this behavior just a few degrees away from the equator. At a fixed convective Rossby number, the superadiabaticity has a latitudinal dependence of ${\cos }^{2}\theta$, meaning that its value at a larger colatitude is reduced with respect to its maximum value at the pole. Therefore, in order to maintain a constant heat flux, the velocity must increase toward the equator. In contrast, the horizontal wavevector is also a monotonically decreasing function of colatitude, with a maximum at the pole and a minimum at the equator. Hence, the heat-carrying motions are of a larger scale at larger colatitudes than at the pole as expected from theoretical considerations (e.g., Busse 2002; Dormy et al. 2004) and as seen in numerical simulations of spherical rotating convection in both fully convective domains and in spherical shells (e.g., Glatzmaier 1985; Brun et al. 2005; Miesch 2005; Browning 2008; Augustson et al. 2016). These behaviors are illustrated in Figure 2(b).

The procedure developed above can be applied again in the case with thermal diffusion as the heat-flux maximization condition remains unchanged from Equation (39). Taking the limit that $V\to 0$ in Equation (29), one finds that

$$\begin{eqnarray}&&\displaystyle \frac{d\hat{s}}{{dz}}=-\displaystyle \frac{3{z}^{2}\left({q}^{3}{K}_{0}{O}_{0}^{2}{\cos }^{2}\theta +\hat{s}\left({\hat{s}}^{2}+2{K}_{0}{{qz}}^{3}\hat{s}-1\right)\right)}{1+{q}^{2}{O}_{0}^{2}{\cos }^{2}\theta +{z}^{3}\left(3{\hat{s}}^{2}+2{K}_{0}{{qz}}^{3}\hat{s}-1\right)},\end{eqnarray} \tag{ 47 }$$

which implies the following heat-flux maximization condition:

$$\begin{eqnarray}&&\begin{array}{l}\hat{s}\left(1+{q}^{2}{O}_{0}^{2}{\cos }^{2}\theta +2{z}^{3}\left(3{\hat{s}}^{2}+4{{qK}}_{0}{z}^{3}\hat{s}-2\right)\right)\\ \qquad +3{q}^{3}{K}_{0}{O}_{0}^{2}{z}^{3}{\cos }^{2}\theta =0.\end{array}\end{eqnarray} \tag{ 48 }$$

The dispersion relationship (Equation (26)) with ${V}_{0}\to 0$ is

$$\begin{eqnarray}&&\left(\hat{s}+{K}_{0}{{qz}}^{3}\right)\left({\hat{s}}^{2}{z}^{3}+{O}_{0}^{2}{q}^{2}{\cos }^{2}\theta \right)-\hat{s}\left({z}^{3}-1\right)=0.\end{eqnarray} \tag{ 49 }$$

Again utilizing the maximized heat-flux constraint $\hat{s}=\tilde{s}{qz}$ and the definitions in Equations (25), the two previous equations can be shown to be equivalent to

$$\begin{eqnarray}&&\begin{array}{l}{q}^{2}-\displaystyle \frac{\tilde{s}\left({z}^{3}-1\right)}{\left(\tilde{s}+{K}_{0}{z}^{2}\right)\left({\tilde{s}}^{2}{z}^{5}+{O}_{0}^{2}{\cos }^{2}\theta \right)}=0,\\ \quad \tilde{s}z\left(2{z}^{5}-5{z}^{2}-\displaystyle \frac{3}{{\tilde{s}}^{2}}{O}_{0}^{2}{\cos }^{2}\theta \right)\end{array}\end{eqnarray} \tag{ 50 }$$

$$\begin{eqnarray}&&+{K}_{0}\left[4{z}^{8}-7{z}^{5}-\displaystyle \frac{{O}_{0}^{2}}{{\tilde{s}}^{2}}{\cos }^{2}\theta \left({z}^{3}+2\right)\right]=0.\end{eqnarray} \tag{ 51 }$$

So the connection to the nondiffusive case is clear and can be recovered directly when taking the limit as ${K}_{0}\to 0$. Employing the relationship between O and the convective Rossby number given in Equation (23), the superadiabaticity and wavevector can be defined as

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{\epsilon }{{\epsilon }_{0}}=\displaystyle \frac{\left(\tilde{s}+{K}_{0}{z}^{2}\right)\left(25{\pi }^{2}{\mathrm{Ro}}_{{\rm{c}}}^{2}{\tilde{s}}^{2}{z}^{5}+6{\cos }^{2}\theta \right)}{25{\pi }^{2}{\mathrm{Ro}}_{{\rm{c}}}^{2}\tilde{s}\left({z}^{3}-1\right)},\\ \tilde{s}z\left[2{z}^{5}-5{z}^{2}-\displaystyle \frac{18{\cos }^{2}\theta }{25{\pi }^{2}{\mathrm{Ro}}_{{\rm{c}}}^{2}{\tilde{s}}^{2}}\right]\end{array}\end{eqnarray} \tag{ 52 }$$

$$\begin{eqnarray}&&+{K}_{0}\left[4{z}^{8}-7{z}^{5}-\displaystyle \frac{6{\cos }^{2}\theta \left({z}^{3}+2\right)}{25{\pi }^{2}{\mathrm{Ro}}_{{\rm{c}}}^{2}{\tilde{s}}^{2}}\right]=0.\end{eqnarray} \tag{ 53 }$$

Immediately, one can see that the superadiabaticity increases as the thermal diffusion is increased because the thermal gradient must increase to drive convection. The nondiffusive component is unchanged and so it also grows with lower convective Rossby numbers. The form of Equation (53) indicates that the wavevector must increase with decreasing convective Rossby number, whereas its behavior with varying thermal diffusion rate (K₀) is not obvious and is illustrated in Figure 4(a) for the fully diffusive case. There, it is seen that thermal diffusion induces a shift in the value of the wavevector ratio that is a factor of ${(7/10)}^{1/3}$, or approximately 13%, lower for convective Rossby numbers below the value of K₀. For ${K}_{0}\gt 1$, the asymptotic value of z at large convective Rossby number is reduced by the same factor, from ${(5/2)}^{1/3}$ to ${(7/4)}^{1/3}$, indicating that the horizontal scale of the flows becomes larger. The velocity scales inversely with the horizontal wavevector; so, as before, the velocity must decrease with increasing rotation rate. These behaviors are shown for a typical value of the thermal diffusion in Figure 3(a). The velocity decreases with lower convective Rossby number given that the wavevector increases, with the same nontrivial behavior with the thermal diffusivity.

Zoom In Zoom Out Reset image size

When including viscosity within this heuristic framework, the maximum heat flux should formally be derived from the algebraic system formed by the dispersion relationship for $\hat{s}$ and from

$$\begin{eqnarray}&&F/{F}_{0}={\tilde{s}}^{3}\left[\displaystyle \frac{\hat{s}}{{q}^{3}{z}^{3}}+{V}_{0}q\right]{\hat{s}}^{2}=1.\end{eqnarray} \tag{ 54 }$$

In particular, the above equation and Equations (9) and (28) give

$$\begin{eqnarray}&&{\hat{s}}^{3}+\left({V}_{0}q{\hat{s}}^{2}-{\tilde{s}}^{3}\right){q}^{3}{z}^{3}=0,\end{eqnarray} \tag{ 55 }$$

$$\begin{eqnarray}&&\begin{array}{l}\left(\hat{s}+{K}_{0}{{qz}}^{3}\right)\left[{z}^{3}{\left(\hat{s}+{V}_{0}{{qz}}^{3}\right)}^{2}+\displaystyle \frac{6{\cos }^{2}\theta }{25{\pi }^{2}{\mathrm{Ro}}_{{\rm{c}}}^{2}}\right]\\ \qquad +\left(1-{z}^{3}\right)\left(\hat{s}+{V}_{0}{{qz}}^{3}\right)=0,\end{array}\end{eqnarray} \tag{ 56 }$$

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{6{q}^{2}{\cos }^{2}\theta }{25{\pi }^{2}{\mathrm{Ro}}_{{\rm{c}}}^{2}}\left[{\hat{s}}^{2}+{K}_{0}{{qz}}^{3}\left(3\hat{s}+2{{qV}}_{0}{z}^{3}\right)\right]\\ \qquad +\left(\hat{s}+{V}_{0}{{qz}}^{3}\right)\left[6{\hat{s}}^{3}{z}^{3}+4q{\hat{s}}^{2}{z}^{6}\left(3{V}_{0}+2{K}_{0}\right)\right.\\ \qquad +\hat{s}\left(1-4{z}^{3}+2{V}_{0}{q}^{2}{z}^{9}\left(8{K}_{0}+3{V}_{0}\right)\right)\\ \qquad \left.+2{V}_{0}{{qz}}^{3}\left(1-2{z}^{3}+4{K}_{0}{V}_{0}{q}^{2}{z}^{9}\right)\right]=0.\end{array}\end{eqnarray} \tag{ 57 }$$

While this set of equations can be solved numerically, they can be simplified in the case of planetary and stellar interiors to a relatively simple scaling behavior with viscosity. In particular, the viscous component of the heat flux may be ignored if it is again assumed that ${V}_{0}\ll {N}_{* ,0}$), so that as above $\hat{s}\approx \tilde{s}{qz}+{ \mathcal O }({V}_{0})$. Subsequently, following the method shown for the case of thermal diffusion, one may find the implicit wavevector derivative of the growth rate $\hat{s}$ (Equation (29)), which implies that the constraining dispersion relationship and the equation resulting from the heat-flux maximization are

$$\begin{eqnarray}&&\begin{array}{l}\left(\tilde{s}+{V}_{0}{z}^{2}\right)\left[{q}^{2}{z}^{5}\left(\tilde{s}+{V}_{0}{z}^{2}\right)\left(\tilde{s}+{K}_{0}{z}^{2}\right)-{z}^{3}+1\right]\\ \qquad +\,\displaystyle \frac{6{\cos }^{2}\theta {q}^{2}\left({K}_{0}{z}^{2}+\tilde{s}\right)}{25{\pi }^{2}{\mathrm{Ro}}_{{\rm{c}}}^{2}}=0,\end{array}\end{eqnarray} \tag{ 58 }$$

$$\begin{eqnarray}&&\begin{array}{l}\tilde{s}\left(1-4{z}^{3}\right)+3{V}_{0}{z}^{2}\left(1-2{z}^{3}\right)+\displaystyle \frac{6{\cos }^{2}\theta \left(\tilde{s}+3{K}_{0}{z}^{2}\right)}{25{\pi }^{2}{\mathrm{Ro}}_{{\rm{c}}}^{2}}\\ +2{z}^{5}\left(\tilde{s}+{V}_{0}{z}^{2}\right)\left[3{\tilde{s}}^{2}+\tilde{s}\left(4{K}_{0}+5{V}_{0}\right){z}^{2}+6{K}_{0}{V}_{0}{z}^{4}\right].\end{array}\end{eqnarray} \tag{ 59 }$$

So one may solve for q from the former equation, the solution of which upon substitution into the latter equation yields an equation solely for the wavevector z:

$$\begin{eqnarray}&&\begin{array}{l}{z}^{3}{\left({V}_{0}{z}^{2}+\tilde{s}\right)}^{2}\left[3{V}_{0}{K}_{0}{z}^{4}\left(2{z}^{3}-3\right)\right.\\ \qquad \left.+\tilde{s}{z}^{2}\left({V}_{0}+{K}_{0}\right)\left(4{z}^{3}-7\right)+{\tilde{s}}^{2}\left(2{z}^{3}-5\right)\right]\\ \qquad -\displaystyle \frac{6{\cos }^{2}\theta }{25{\pi }^{2}{\mathrm{Ro}}_{{\rm{c}}}^{2}}\left[2\tilde{s}\left({K}_{0}-{V}_{0}\right)+3{\tilde{s}}^{2}z\right.\\ \qquad \left.+\tilde{s}\left({K}_{0}+5{V}_{0}\right){z}^{3}+3{K}_{0}{V}_{0}{z}^{5}\right]=0.\end{array}\end{eqnarray} \tag{ 60 }$$

The horizontal wavevector is then recovered from the roots of the fourteenth-order Equation (60), whereas the superadiabaticity is defined as

$$\begin{eqnarray}&&\displaystyle \frac{\epsilon }{{\epsilon }_{0}}=\displaystyle \frac{\left(\tilde{s}+{K}_{0}{z}^{2}\right)\left(25{\pi }^{2}{\mathrm{Ro}}_{{\rm{c}}}^{2}{\tilde{s}}^{2}{z}^{5}{\left(\tilde{s}+{V}_{0}{z}^{2}\right)}^{2}+6{\cos }^{2}\theta \right)}{25{\pi }^{2}{\mathrm{Ro}}_{{\rm{c}}}^{2}\tilde{s}\left({z}^{3}-1\right)\left(\tilde{s}+{V}_{0}{z}^{2}\right)}.\end{eqnarray} \tag{ 61 }$$

Clearly, the inclusion of viscosity increases the degree of superadiabaticity beyond that due to the thermal diffusion. Moreover, as before, it also grows with lower convective Rossby numbers. Given the form of Equation (60), the wavevector must increase with decreasing convective Rossby number, whereas again its behavior with varying K₀ is not obvious, but is relatively benign as is shown in Figure 4(a). In contrast, the viscosity has a much larger impact on the amplitude of the change in the wavevector as exhibited in Figure 4(b). The velocity decreases with lower convective Rossby number, given that the wavevector increases, with the same nontrivial behavior seen earlier in the case with thermal diffusion.

Zoom In Zoom Out Reset image size

However, for values of $\left\{{K}_{0},{V}_{0}\right\}\lt {N}_{* ,0}$, one may find a heuristic fit to the horizontal wavevector for the diffusive case given by

$$\begin{eqnarray}&&\begin{array}{l}z\approx {2}^{-\tfrac{11}{6}}{5}^{\tfrac{1}{3}}\coth {\left[3{\mathrm{Ro}}_{{\rm{c}}}^{\tfrac{2\sqrt{2}}{5}}\right]}^{\tfrac{\sqrt{2}}{2}}\tanh \left[\displaystyle \frac{17}{10}{\mathrm{Ro}}_{{\rm{c}}}^{\tfrac{23}{124}}{V}_{0}^{-\tfrac{37}{160}}\right]\\ \left[2+\left(\sqrt{2}-1\right)\left(1+\tanh \left[\displaystyle \frac{\sqrt{2}}{2}\left(\displaystyle \frac{3}{2}\right.\right.\right.\right.\\ +\left.\left.\left.\left.{\mathrm{log}}_{10}\displaystyle \frac{{\mathrm{Ro}}_{{\rm{c}}}}{{\left({V}_{0}+{K}_{0}\right)}^{\tfrac{7}{6}}}\right)\right]\right)\right],\end{array}\end{eqnarray} \tag{ 62 }$$

which can be used more directly to estimate the scaling of the relative superadiabaticity and velocity amplitudes. The form of the fit follows from the observation of the qualitative features of the scaling of the horizontal wavevector given successively in Figures 2–4, where it is well described by hyperbolic tangents to within a few percent even in the transitional region of ${\mathrm{Ro}}_{{\rm{c}}}\approx 1/10$. The accuracy of the fit is illustrated in Figure 4(b).

In stellar models, a diffusive treatment of transport processes can be adopted for thermodynamics, chemicals, angular momentum, and magnetic field (e.g., Heger et al. 2000; Maeder 2009; Mathis 2013). Within a convection zone, the turbulent mixing of those quantities can be approximated through a parameterized vertical diffusion ${D}_{{\rm{v}}}=1/3\,\alpha {H}_{P}{v}_{\mathrm{MLT}}$. Because this depends on the pressure scale height H_P, it is worth noting that it impacts the superadiabaticity as well. The definition of the Brunt–Väisälä frequency is ${N}^{2}=g{\alpha }_{T}T\epsilon /{H}_{P}$. Thus, allowing H_P to vary, one has ${N}^{2}/{N}_{0}^{2}={{xH}}_{P,0}/{H}_{P}=x/h$. To maintain the constancy of the heat flux, this then requires $x=\epsilon /{\epsilon }_{0}$ to depend linearly on h, with $x\to {xh}$.

Within the context of this convection model and for a constant mixing-length parameter α, the diffusion coefficient ratio scales as

$$\begin{eqnarray}&&\displaystyle \frac{{D}_{v}}{{D}_{0}}=h\displaystyle \frac{v}{{v}_{0}}={\left(\displaystyle \frac{5}{2}\right)}^{\tfrac{1}{6}}\displaystyle \frac{h}{\sqrt{z}},\end{eqnarray} \tag{ 63 }$$

where D₀ is the local diffusion coefficient in the absence of rotation. The scaling of ${D}_{v}$ is shown in Figure 6 and described in detail in the following section, as it suffices to note that the diffusion will increase with the scale-height ratio h and decrease as the rotation rate is increased. This convection model does, however, exclude any horizontal diffusion associated with the increased horizontal shear typically present in rotating convection.

Standard MLT is not able to address convective penetration directly; rather, it must be modified to account for this nonlocal effect. Such nonlocalities arise naturally in the generalizations of MLT as they directly model nonlinear interactions. One of the first and simplest such nonlocal generalizations of MLT can be found in Spiegel (1963). There, the convective heat flux was treated as an unknown distribution in position and velocity space, and as such it follows the dynamics described by a collisional Boltzmann equation. The collisional processes are a stochastic source arising from nonlinear processes and a damping term. The latter is the connection to MLT as it limits the path length of a convective eddy to the mixing length ${{\ell }}_{\mathrm{MLT}}$, but it is not required to be small in some sense as it needs to be in standard mixing-length theory. Another approach is to include the kinetic energy flux, which leads to a convection theory that depends upon the entropy profile in the star and thus upon global integral constraints on the energy flux in the star (e.g., Roxburgh 1978). Other models that variously extend MLT to include nonlinear processes through finite amplitude analyses have also been developed (e.g., Veronis 1963; Sparrow et al. 1964). Subsequently, modal expansions of the Boussinesq and anelastic equations provided solutions that established that nonlinear penetration was substantially deeper than what linear theory predicted (e.g., Musman 1968; Moore & Weiss 1973). In anelastic convection, the density stratification induces asymmetries between upflows and downflows that further increase the penetration depth (e.g., Zahn et al. 1982; Massaguer et al. 1984). Yet as is often the case, the penetration depth determined in these models depends sensitively on their assumptions, with a wide disparity in the calculated penetration depths. The problems and inconsistencies with many of these models were reviewed in Renzini (1987).

Three pioneering studies analytically estimated the convective penetration depth by prescribing the convective flows as a boundary condition and examining the impact on the region of penetration. One of the first models to attempt to estimate the importance of overshooting above a convective core is from Roxburgh (1965), in which ergodicity and isotropy are assumed, leaving a monodimensional and time-independent equation of motion. A linearization of the thermodynamics with respect to the vertical displacement then permits the equation of motion to be integrated to yield an estimate of the depth of penetration of a fluid element. Building upon the work of Roxburgh (1978), Schmitt et al. (1984) used analytical models of turbulent plume ensembles and examined their statistical properties, finding that the penetration depth scales as ${v}_{p}^{3/2}{f}^{1/2}$, where ${v}_{p}$ is the plume velocity and f is the fractional area filled by the plumes. Likewise, following Roxburgh (1965, 1978) and by directly parameterizing the filling factor f, Zahn (1991, hereafter Z91) found a similar scaling result through a spatial linearization of the equations of motion and the energy transport. Both of these studies predicted penetration depths as fairly large fractions of the pressure scale height. Z91 also confirmed the modal results of Massaguer et al. (1984), where the flow asymmetry due to the density stratification causes downflows to penetrate farther into an underlying stable layer than upflows penetrate into an overlying layer. As a first estimate of the convective penetration depth, one may consider the Z91 linearized model for convective penetration. This model ignores the effects of rotation in both the physical description of the model as well as in its parameterization of the convective dynamics. The former modification will be saved for later work. However, the extended version of the S79 model derived above provides a means of estimating the latter alteration. In particular, one can harness the linearized framework of the Z91 model to obtain an order-of-magnitude estimate of the depth of penetration while allowing for the effects of rotation to be included through a modified value of the superadiabaticity and mixing-length velocity within the convective region. It is also a formalism that may be easily implemented in stellar evolution codes.

3.2.1. Theory

The Z91 model has four phenomenologically distinct zones: a convective zone that ends when the radiative energy flux becomes equal to the convective energy flux, a subadiabatic penetrative zone where flows render the region nearly adiabatic, a transitional thermal boundary layer with overshooting convection where the Péclet number is reduced to below unity, and a radiative zone where waves and conduction carry the total energy flux. There are a few basic assumptions underlying the model. The first assumption is that the velocity and the temperature share the same planform, namely that they are perfectly correlated. The second assumption is that only the downflows for downward penetrating flows and upflows for upward penetrating flows are effective at carrying enthalpy, which is parameterized through their filling factor f, which results from a horizontal average over their planform ${v}_{z}(x,y,z)=v(z)h(x,y)$ such that $f=\overline{{h}^{2}(x,y)}$. From simulations, the value of f appears to be about 1/3, with slight variations of order of 10% that are sensitive to the level of turbulence, rotational, and magnetic effects (Brummell et al. 2002; Stein et al. 2009a, 2009b). However, rotation does impact the structure of the convection as a function of colatitude, whose structure in turn depends upon the degree of supercriticality of the system. For rotating systems, the convective structures should possess a modicum of alignment with the rotational axis, which implies that the mechanism of penetration will vary with colatitude. For instance, consider that columnar flow structures will mostly conform to the classical paradigm of radially aligned penetrating flows near the pole, but they will increasingly become horizontally shearing, yet still localized, flows at lower colatitudes. This will be combined with any shear-induced mixing by large-scale flows, such as differential rotation (e.g., Zahn 1992; Maeder 2003; Mathis et al. 2004; Mathis 2013). However, for the purposes of this model, the filling factor f shall be considered to be a fixed parameter of the theory. Additionally, it shall be assumed that the local magnitude of the gravity g, the local conductivity κ, the thermal expansion coefficient ${\alpha }_{T}$, the nuclear energy generation rate Q, and the specific heat capacity at constant pressure c_P are unaffected by rotation, namely that the centripetal acceleration may be ignored to first order.

A further assumption is that the convective enthalpy flux can be linearized in the region of penetration such that

$$\begin{eqnarray}&&{F}_{\mathrm{conv}}=\rho {c}_{P}\overline{{v}_{z}T^{\prime} }=f\rho {c}_{P}v\delta T\end{eqnarray} \tag{ 64 }$$

$$\begin{eqnarray}&&=\,{F}_{0}{\left[\displaystyle \frac{d\mathrm{ln}M}{d\mathrm{ln}r}-\displaystyle \frac{d\mathrm{ln}L}{d\mathrm{ln}r}+\displaystyle \frac{d\mathrm{ln}\chi }{d\mathrm{ln}r}\right]}_{{r}_{0}}\displaystyle \frac{z}{{r}_{0}},\end{eqnarray} \tag{ 65 }$$

where M(r) is the mass, L(r) is the luminosity, and $\chi (r)$ is the radiative conductivity such that $\chi =\rho {c}_{P}\kappa$ as seen in Z91 (Equation (4.4)). Finally, another assumption is that the vertical inertial force balances the buoyancy force on average such that the pressure contribution vanishes and so that the density perturbations can be linearized with $\rho ^{\prime} ={\alpha }_{T}T^{\prime}$. To follow Z91 as closely as possible, the system is considered only at the pole so that the direct effects of the local Coriolis acceleration $2{{\rm{\Omega }}}_{0}\sin \theta {v}_{x}$ can be ignored, which is equivalent to having a buoyant braking timescale that is shorter than the rotational timescale. Instead, the Coriolis effect implicitly influences the penetration depth by modifying the upper boundary value of the velocity, which is taken from the convection zone. Thus, when horizontally averaged and a pressure equilibrium is assumed as in Z91 (Equation (3.8)), the force balance becomes

$$\begin{eqnarray}&&\displaystyle \frac{c}{2}\displaystyle \frac{{{dv}}^{2}}{{dz}}=g\displaystyle \frac{\delta T}{T},\end{eqnarray} \tag{ 66 }$$

where $c=\overline{{h}^{3}}/\overline{{h}^{2}}$. This equation can then be linked to the convective energy flux through Equation (64) and integrated across the penetrative region to yield an estimate for the penetration depth L_P relative to the pressure scale height H_P.

In Z91, there are two cases for the scaling of the convective penetration depth with velocity and spherically symmetric thermodynamic quantities: one for penetrative convection into a stable region below a convection zone, such as what takes place near the base of the solar convection zone, and one for penetrative convection into a stable region above the convection zone, such as what takes place near the convective cores of intermediate-mass and high-mass stars. The first of these two regimes ignores the variation in the total luminosity and mass and so only retains the term related to changes in the radiative conductivity. It scales as

$$\begin{eqnarray}&&\displaystyle \frac{{L}_{P}}{{H}_{P}}={\left[\displaystyle \frac{2}{3}\displaystyle \frac{\left(1-f\right){{fv}}_{z}^{3}}{g{\alpha }_{T}\kappa {\chi }_{P}{{\rm{\nabla }}}_{\mathrm{ad}}}\right]}^{\tfrac{1}{2}},\end{eqnarray} \tag{ 67 }$$

where ${{\rm{\nabla }}}_{\mathrm{ad}}$ is the adiabatic temperature gradient and ${\chi }_{P}\,=\partial \mathrm{ln}\kappa /\partial \mathrm{ln}P{| }_{S}$ is the adiabatic logarithmic derivative of the radiative conductivity with respect to pressure. The opposite case of convective penetration into a stable region above a convection zone scales as

$$\begin{eqnarray}&&\displaystyle \frac{{L}_{P}}{{H}_{P}}={\left[\displaystyle \frac{2\left(1-f\right){{fv}}_{z}^{3}{r}_{0}}{9g{\alpha }_{T}\kappa {{\rm{\nabla }}}_{\mathrm{ad}}{H}_{P}\left(\rho P/\overline{\rho P}-\rho Q/\overline{\rho Q}-{r}_{0}{\chi }_{P}/3{H}_{P}\right)}\right]}^{\tfrac{1}{2}},\end{eqnarray} \tag{ 68 }$$

where r₀ is the radius at the edge of the core. Note that ${{\rm{\nabla }}}_{\mathrm{ad}}={\rm{\nabla }}+\epsilon$, with ∇ being the temperature gradient and being the superadiabatic gradient as above. However, the basic assumption of the model is that does not grow large enough to modify the background temperature gradient in a steady state, and so its variation with rotation rate does not strongly influence the depth of penetration. Thus, the ratio of the penetration depth with rotation and diffusion to the nonrotating inviscid value for convective penetration into a stable layer either above or below a convection therefore scales as

$$\begin{eqnarray}&&\displaystyle \frac{{L}_{P}}{{L}_{P,0}}={\left(\displaystyle \frac{v}{{v}_{0}}\right)}^{3/2}={\left(\displaystyle \frac{5}{2}\right)}^{\tfrac{1}{4}}{z}^{-\tfrac{3}{4}}.\end{eqnarray} \tag{ 69 }$$

As seen in the previous section, the velocity amplitude of the mode that maximizes the heat flux decreases with lower diffusivities and lower convective Rossby numbers. Therefore, the penetration depth necessarily must decrease when the convective Rossby number is decreased. This behavior follows intuitively given that the reduced vertical momentum of the flows implies that the temperature perturbations are also reduced (via Equation (13)). Thus, due to the decreased buoyant thermal equilibration time and the reduced inertia of the flow, the penetration depth must decrease. In contrast, the velocity and the horizontal scale of the flow increase with greater diffusivities in order to offset the reduced temperature perturbations in the case of larger thermal conductivity. In the case of larger viscosity, the horizontal scale of the velocity field is increased, whereas for a fixed thermal conductivity, the thermal perturbations are of a smaller scale. Thus, to maintain the heat flux, the amplitude of the velocity must increase in order to compensate for the reduced correlations between the two fields. The scaling behaviors of the penetration depth are illustrated as a function of diffusivities and convective Rossby number in Figure 5.

Zoom In Zoom Out Reset image size

3.2.2. Comparison with Penetrative Convection Experiments

This model of penetrative convection can be introduced into stellar models utilizing a diffusive approach. Such a parameterization of mixing processes has been extensively examined (e.g., Freytag et al. 1996; Herwig et al. 1999; Denissenkov et al. 2013; Viallet et al. 2015; Lecoanet et al. 2016). A complementary model has been established through an extreme-value statistical analysis of 3D penetrative convection simulations (Pratt et al. 2017b). A frequent assumption made of turbulent flows is to assume that the underlying distribution of a flow quantity is normal. However, it was shown through the simulations and analysis of Pratt et al. (2017a, 2017b) that the rarer events occurring in the tails of the actual distribution function, i.e., the extremal penetrative flows with a higher velocity and a greater entropy deficit than the average flow, have a larger impact than is expected in Gaussian turbulence models on the turbulent mixing coefficients. This analysis was carried out with a novel method to quantify this non-Gaussianity of the dynamics. Their method permits a statistical examination of the turbulent particle dispersion and subsequently the construction of a model for turbulent diffusion based upon the Gumbel distribution (Pratt et al. 2017a). This model for diffusion and overshoot has been applied to the problem of lithium depletion occurring in rotating low-mass stars (Baraffe et al. 2017), where it was found to mimic the observed trends well when the effects of rotation on the diffusion treatment were introduced. There, the rotational effects are included as an empirically defined threshold rotation rate above which the diffusion is quenched ostensibly, due to the action of rotation on the convective flows. Also, the penetrative convective depth was a parameter and was not dynamically estimated. However, the initial estimate of Baraffe et al. (2017) may be improved. Using the above extension of the Z91 model, one can estimate both the penetration depth and the level of diffusion. Taking the parameters of the Gumbel distribution as in Pratt et al. (2017a) yields the following description of the radial dependence of the diffusion coefficient

$$\begin{eqnarray}\begin{array}{c}\begin{array}{rcl}{D}_{v}\left(r\right) & = & {\left(\displaystyle \frac{5}{2}\right)}^{\tfrac{1}{6}}\displaystyle \frac{\alpha {H}_{P,0}{v}_{0}h}{3\sqrt{z}}\\ & & \times \,\left\{1-\exp \left[-\exp \left(\left(r-{r}_{c}\right)/\lambda {L}_{P}+\mu /\lambda \right)\right]\right\},\end{array}\end{array}\end{eqnarray} \tag{ 70 }$$

where r_c is the base of the convection zone and where μ and λ are the empirically determined parameters from Baraffe et al. (2017). An illustration of the scaling behavior due both to the variation of ${D}_{v}$ with the mixing-length velocity $v$ and with the depth of penetration L_P is shown in Figure 6, where the parameters are taken to make contact with Sun-like stars where the transition region begins around $r\approx 0.7\,{R}_{\odot }$. Note, however, that the convective velocity will vary in the convection zone. So, Figure 6 is only meant to illustrate the rough behaviors of this diffusive model near the convectively stable-to-unstable transition point. From Equation (70) and as seen in Figure 6, the radial structure of the diffusion coefficient follows from the scaling of the velocity, namely the diffusion will globally decrease with decreasing convective Rossby number and will further vary if the local angular velocity is taken into account. The depth of penetration is perhaps most notable, in that its strong rotational dependence can lead to severe restrictions on the region in which the diffusion acts. This is evidenced through the transition between a long diffusive tail at high convective Rossby number to a step-like function at low convective Rossby number. Thus, as discussed earlier, the amount of and the depth over which mixing induced above or below a convective region may be severely reduced in stars whose convection possesses a low convective Rossby number.

Zoom In Zoom Out Reset image size

In Lecoanet et al. (2016), the mixing induced below a convection zone is characterized as a means of determining whether or not carbon flames can be disrupted by overshooting convection in evolved massive stars of between about 7–11 solar masses. To estimate the amount of mixing, the results of an idealized model of these flames are assessed with 3D hydrodynamic simulations of the Boussinesq equations. In particular, a passive scalar field that permits the measurement and fitting of a height-dependent turbulent chemical diffusivity that mimics the effect of the overshooting convection is evolved. This diffusion coefficient is modeled as a composition of two error functions, which has a super-exponential tail in the mixing region somewhat similar to the model constructed above.

What is currently lacking, however, are extensive parameter studies using 3D simulations of convective overshoot that also include rotation. Such simulations would at least provide numerical evidence regarding the hypothesized dependence of the mixing depth and its shape. For instance, the simulations carried out in Korre et al. (2018) have also addressed the shape of the overshooting region in a Boussinesq system and its associated potential mixing properties in 3D simulations. As in Lecoanet et al. (2016), it is found that the horizontally averaged kinetic energy is super-exponentially reduced with depth below the convection zone. Specifically, in those simulations, the average kinetic energy is best modeled with a Gaussian. Moreover, the radial correlation length scale of downflows is larger than the length scale associated with the averaged kinetic energy profile. This length scale may be interpreted as the depth to which the strongest downflows travel before stopping. This behavior is at least roughly consistent with the model proposed here and with the results of Pratt et al. (2017a, 2017b), although it is in contrast to the models proposed in Freytag et al. (1996) and Herwig (2000). For the time being, the examinations in Julien et al. (1996) and Brummell et al. (2002) will have to suffice as evidence that the overshooting is reduced with rotation rate, while the shape of the mixing region likely depends upon the model setup and equations solved in similar simulations.