A model of rotating convection in stellar and planetary interiors: I - convective penetration
Kyle C. Augustson Stéphane Mathis


 

K. C. Augustson 
kyle.augustson@cea.fr



AIM, CEA, CNRS, Université Paris-Saclay, Université Paris Diderot, Sarbonne Paris Cité, F-91191 Gif-sur-Yvette Cedex, France


  A monomodal model for stellar and planetary convection is derived for the magnitude of the rms velocity, degree of
  superadiabaticity, and characteristic length scale as a function of rotation rate as well as with thermal and viscous
  diffusivities. The convection model is used as a boundary condition for a linearization of the equations of motion in
  the transition region between convectively unstable and stably-stratified regions, yielding the depth to which
  convection penetrates into the stable region and establishing a relationship between that depth and the local
  convective Rossby number, diffusivity, and pressure scale height of those flows.  Upward and downward penetrative
  convection have a similar scaling with rotation rate and diffusivities, but they depend differently upon the pressure
  scale height due to the differing energetic processes occurring in convective cores of early-type stars versus
  convective envelopes of late-type stars.









§ INTRODUCTION


In the context of stellar and planetary physics, convection driven through buoyancy and doubly-diffusive instabilities
plays a critical role in mixing and transport processes in convectively unstable regions of stars
<cit.>.  It primarily serves to transport the energy released deep within the star
or planet through regions where radiative energy transport is inefficient.  Moreover, convection directly and nonlocally
transports heat and chemicals through advection and it acts diffusively through entrainment and dissipative processes.
In the presence of rotation, convection can also transport angular momentum through Reynolds stresses and meridional
flows to establish and maintain a differential rotation
<cit.>. However,
parameterizing the impact of rotation on the convection and its related transport properties over evolutionary time
scales remains a relatively open problem.

Convective flows cause mixing not only in regions of superadiabatic temperature gradients but in neighboring
subadiabatic regions as well, as motions from the convective region contain sufficient inertia to extend into those
regions before being buoyantly braked or turbulently eroded
<cit.>.  This convective penetration and turbulence can
thus alter the chemical composition and thermodynamic properties in those regions, softening the transition between
convectively stable and unstable regions, with the consequence being that the differential rotation, opacity, and
thermodynamic gradients are modified
<cit.>. Such processes
have an asteroseismic signature as has been observed in many kinds of stars
<cit.>. Indeed, penetrating convection
can modify the depth of the convection zone and any extant tachocline leading to a frequency shift in asteroseismically
detected modes <cit.> and advect magnetic field and angular
momentum into the stable region leading to changes in dynamo action and the differential rotation
<cit.>. In stars with a
convective core, convective overshoot and penetration can lead to a greater amount of time spent on stable burning
phases as fresh fuel is mixed into the burning region
<cit.>.  Penetrative convection induces chemical mixing
near the transition between convectively stable and unstable regions of stars, leading to changes in the spectral
characteristics of the atmosphere <cit.>. Hence, from the
standpoint of stellar evolution, an open problem is to understand how the depth of penetration and the character of
the convection in this region change with rotation, magnetism, and diffusion.

In this work, a generalized version of a heuristic model for convection for rotating systems is developed following
<cit.>. This model of convection is then employed to estimate the convective depth of penetration above and
below a convection zone. The linearized Boussinesq equations of motion that yield the growth rate used in the model are
given in <ref>.  The heat-flux maximization principle is discussed in <ref>. The model of
convection is derived for a diffusion-free setting in <ref>, when including thermal diffusion in
<ref>, and when taking into account both thermal and viscous diffusion in <ref>.  This
convection model is then used in conjunction with the linearized model of convective penetration derived in
<cit.> to estimate the depth of convective penetration as shown in <ref>.  A summary of the
results and perspectives are presented in <ref>.



§ HEAT-FLUX MAXIMIZED CONVECTION MODEL


As a first step, the heuristic model will be considered to be local such that the length scales of the flow are much
smaller than either density or pressure scale heights. This is equivalent to ignoring the global dynamics and assuming
that the convection can be approximated as local at each radius and colatitude in a star or planet.  As such, one may
consider the dynamics to be in the Boussinesq limit (Section <ref>).  The linearized dynamics is used to
construct the heat-flux maximization principle.



 §.§ Linearized Boussinesq Dynamics
 

Boussinesq devised one of the simplest convective systems, consisting of an infinite layer of a nearly incompressible
fluid with a small thermal expansion coefficient α_T=-∂lnρ/∂ T|_P that is confined between
two infinite impenetrable plates differing in temperature by Δ T=T(z_2)-T(z_c), where it is assumed for this model that
T(z_2)<T(z_c), and separated by a distance ℓ_0=z_2-z_c, as in Figure <ref>. The thermodynamic
variables are further expanded about their averages as q = ⟨ q⟩ + q̅(z) + q', with ⟨ q⟩
being the volumetric mean, q̅ being the horizontal average, and q' being the dynamical perturbation. Following
<cit.> and <cit.>, 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 α_T Δ T. This is equivalent to requiring that the maximum density
fluctuation be small in the sense that α_TΔ T≪ 1.  Under these assumptions, the equations of motion become





  
    ( - ν∇^2)   = --1/⟨ρ⟩ P' + gα_TT' - 2,
    ( - κ∇^2) T'    = β - T',




where t is the time coordinate and ∂_t is the partial derivative with respect to it,  is the
spatial coordinate gradient. To first order in α_T Δ T, the continuity equation requires that the velocity
field  be solenoidal (=0). The vertical direction  is anti-aligned to the local direction of
the constant effective gravity g, e.g.  =-𝐠/g, so that fluid with T'>0 rises when Δ T<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
<ref>). Here ⟨ρ⟩ is the mean density, P' is the pressure perturbation, T' is the
temperature perturbation,  is the vertical component of the velocity field,  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 α_T Δ T is defined as
β=dT̅/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 P dV work in the thermal energy equation.



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 <ref>. Dropping the nonlinear terms, these
operations yield the following set of linearized Boussinesq equations





  
    ( - κ∇^2) T' = β, 
       ( - ν∇^2)ζ             = 2, 
       ( - ν∇^2)∇^2 = gα_T∇_⊥^2T' - 2ζ,




with ζ being the vertical component of the vorticity and ∇_⊥ being the gradient transverse to
the vertical direction as in Equations 79, 84, and 85 of <cit.>. As seen in many papers regarding
Boussinesq dynamics <cit.>, 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





  
    ( - κ∇^2)( - ν∇^2)^2∇^2 +
          gα_Tβ∇_⊥^2( - ν∇^2)
           + 4[( - κ∇^2) ] = 0.




For impenetrable and stress-free boundary conditions the solutions of Equation <ref> are periodic in the
horizontal, sinusoidal in the vertical, and exponential in time, e.g.
= sin[k_z(z-z_c)]exp(i_⊥ + s t), where _⊥ is the
horizontal wavevector, s is the growth rate,  is the local coordinate vector, and  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π/ℓ_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 <ref> yields the following
dispersion relationship that relates s to the wavevector  as





  
    [(s +κ k^2)(s + ν k^2)^2 k^2 + 
                gα_Tβ k_⊥^2(s + ν k^2). 
           .+ 4()^2(s + κ k^2)]  = 0.




This equation can be simplified and made non-dimensional by dividing through by the appropriate powers of N_*
and k_z, where N_*^2=|gα_Tβ| is the absolute value of the square of the Brunt-Väisälä frequency
as in <cit.>, which is otherwise negative in a convective region.  Utilizing the following auxiliary quantities





  
    ŝ   =s/N_*,    
        z^3 = 1+ a^2 = k^2/k_z^2,    
        a^2 = k_x^2/k_z^2+k_y^2/k_z^2 = a_x^2 + a_y^2, 
    
        K    = κ k_z^2/N_*,   
        V = ν k_z^2/N_*,




the dispersion relationship becomes





  
    (ŝ+ K z^3)( z^3(ŝ+ V z^3)^2 + O^2 (cosθ
          + a_ysinθ)^2) 
       -(z^3- 1)(ŝ+ V z^3)=0,




where 4(·)^2/N_*^2= k_z^2 O^2 (cosθ + a_ysinθ)^2 with





  
    O^2    = 4Ω_0^2/N_*^2,




where Ω_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 neglected, which is the case considered in , and
expanded upon below in <ref>. Within stellar interiors, the molecular value of the thermal Prandtl number
(Pr=ν/κ) is typically very small, being of the order of 10^-5 or smaller, so the limit
ν→ 0 may be considered as in <ref>.  The case where thermal diffusion can be ignored
(κ→ 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
(<ref>), 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.



 §.§ Maximum Heat Transport


The secular impacts of rotation and magnetic fields on stellar and planetary evolution are of keen interest within the
astrophysical community <cit.>. Accordingly, extensions to MLT and its ilk have been
introduced <cit.>. As expounded upon in <cit.>, 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 <cit.>, where it is used to set the scale of energy injection of a
Heisenberg-Kolmogorov spectrum of turbulent energy cascade <cit.>. Subsequently, it sets the
local values of the turbulent diffusivities. Though, these theories have not yet widely been employed in stellar
evolution computations, with one exception being <cit.>.  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 <cit.>.

The  model of rotating convection has its origins in the principle of maximum heat transport
proposed by <cit.>.  In that principle, an upper limit for a boundary condition dependent turbulent heat flux
is established that depends upon the smallest Rayleigh unstable convective eddy. The size of this eddy is determined
with a variational technique that is similar to that developed in <cit.> 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. Yet the heat-flux maximization principle is not fully rigorous as there
is no apriori theoretical justification as to why the turbulence would arrange itself so as to maximize the heat flux
<cit.>. 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 <cit.> 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 <cit.>. Indeed, as <cit.> and
<cit.> 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
, 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 <cit.> aim to directly test the heat-flux maximization
principle of <cit.>, 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 <cit.>.  <cit.> finds 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 <cit.>, which
predicts a scaling proportional to the inverse cube root of the vertical coordinate. <cit.> reexamines those
experimental data with an emphasis on testing wether or not both the mean-fields and the fluctuating fields predicted in
<cit.> correspond to the measured flow structures, which are both found to match the theoretical values to
within a factor of order unity. Additionally, <cit.> constructs a variational technique to assess the Rayleigh
number (Ra = -gα_TΔ T ℓ_0^3/νκ) 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 is 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
two. This approach can also be taken to compare Malkus's theory to more recent laboratory and numerical experiments of
Rayleigh-Bénard convection <cit.>. Furthermore, similar experiments of
Rayleigh-Bénard convection in rotating cylinders have measured the Nusselt number scaling with respect to the Rayleigh
number and with the Ekman number [Ek=ν/(2Ωℓ_0^2)] in the rotating cases
<cit.>. Thus, a Malkus-like turbulence theory that
includes rotation, as in  and as shown here, can in principle also be compared to those
experiments, which will be a focus of later work.

Laboratory experiments <cit.> and numerical simulations have lent some credence to this convection model
<cit.>. In particular, those simulations indicate that the low convective Rossby number scaling
regime established in  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 have been used to provide a local model of a rotationally-dependent mixing-length theory
<cit.>. There, it is 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.



 §.§ General Framework


Following , 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 = || ||, which also provides a way to estimate the heat flux
F = ρ c_P [ T'] 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, e.g. n=1 so that k_z = π/ℓ_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 scale-wise 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 <ref>, with the assumptions used above to formulate the dispersion relationship, yields





  
    s  = -i P'/⟨ρ⟩ +gα_TT' -2 - ν k^2, 
       i = 0.




Taking the dot product of  and Equation <ref>, using Equation <ref> to eliminate
the pressure, and taking a volume average it can be seen that





  
    (s+ν k^2) ⟨^2⟩ = gα_T⟨ T'⟩.




Therefore,





  
    F = ρ c_P[(s+ν k^2) ⟨^2⟩]/gα_T.




Employing the hypothesis that the velocity in the quasi-stationary turbulent state scales as
⟨^2⟩ = s s^*/k^2, one has the heat flux





  
    F    = ⟨ρ⟩ c_P/gα_T[(s+ν k^2) ]s s^*/k^2
       = ⟨ρ⟩ c_P N_*^3/gα_T k_z^2[[ŝ]/z^3+V]|ŝ|^2,




where the definitions given in Equations <ref> 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





  
    ŝ^2 = a^2 - O^2(cosθ+sinθa_y)^2/1+a^2,




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>Ocosθ, which is shown to
be satisfied by the heat-flux maximization in <cit.> and in the next subsection.  Likewise, in the
nonrotating but diffusive case, one has that





  
    [ŝ+V z^3][z^3(ŝ+K z^3)(ŝ+V z^3)-a^2]=0.




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





  
    ŝ   = z^-3√((V+K)^2z^12+4 z^3(a^2-K V z^9))
       -1/2(V+K)z^3.




When linearized with respect to V and K, this yields





  
    ŝ   = a z^-3/2-1/2(V+K)z^3,




which is maximized when a_y=0. So, in all the 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
non-rotating Rayleigh-Bénard convection, but nevertheless the 2D modes can approximate many of the behaviors of the
3D modes <cit.>.

It is also useful to recast some of the parameters so that the system can be more readily compared to the
 results. The characteristic velocity _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α_Tβ_0|=3/5N^2_*,0, with
β_0 being the thermal gradient and g_0 being the effective gravity in the nonrotating case, and
k_0^2=5/2 k_z^2 as in <cit.> Equation 36, leading to





  
    _0 = s_0/k_0 = √(6)/5N_*,0/k_z =  √(6)/5π N_*,0ℓ_0.




This permits the definition of the convective Rossby number , which for this Boussinesq system is





  
    =_0/2Ω_0 ℓ_0=√(6) N_*,0/10πΩ_0,




which implies that





  
    O    =2Ω_0/N_* = _0/N_* ℓ_0 = √(6)N_*,0/5π N_*.




Next consider the variation of the superadiabaticity, which for this system is given by ϵ = H_P β/T,
meaning that N_*^2 = |g α_T T ϵ/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





  
    q    = N_*,0/N_*.




Therefore, all parametric quantities have the following equivalencies





  
    O    = q √(6)/5π= q O_0,
    
        K    = q κ k_z^2/N_*,0 = q K_0,
    
        V    = q ν k_z^2/N_*,0 = q V_0.




So, the dispersion relationship (Equation <ref>) and the heat flux (Equation <ref>) may be
written as





  
    (ŝ+ K_0 q z^3)( z^3(ŝ+ V_0 q z^3)^2 + O_0^2q^2cos^2θ) 
           -(z^3- 1)(ŝ+ V_0 q z^3)=0, 
       F = F_0/q^3[ŝ^3/z^3+V_0 qŝ^2],




where F_0 = ⟨ρ⟩ c_P N_*,0^3/(gα_T k_z^2). When maximized for a positive real value of
ŝ, Equation <ref> implies that





  
    dŝ/dz = 3ŝ^2/z(3ŝ+2V_0 q z^3).




Depending upon the method of finding ŝ, 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 <ref>)
to obtain dŝ/dz, which yields





  
    dŝ/dz=-3z^2[V_0q(1-2z^3)+K_0 q^3(O_0^2cos^2θ+4V_0^2z^9) . 
       .+ ŝ(ŝ^2+2q(K_0+2V_0) z^3 ŝ + 3 V_0 q^2(2K_0+V_0)z^6-1)]
       [z^3(3ŝ^2+ 2q(K_0+ 2V_0)z^3ŝ+ V_0q^2(2K_0+ V_0)z^6- 1).
       .+ q^2 O_0^2cos^2θ+1]^-1.




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[F]=max[F]_0 so the heat flux is equal to the maximum value max[F]_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 <cit.>.  Furthermore, the flux will be latitude dependent due to the
rotational variation of the local gravity <cit.>. Indeed, an analysis of 3D nonlinear
global-scale convection simulations of rapidly rotating oblate stars has been carried out in <cit.>.  These
simulations indicate that the Von Zeipel theorem <cit.> 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 <cit.> and the von Zeipel theorem, one may show that
the heat flux F_V on an isobar with P=⟨ P⟩ varies as





  
    F_V = L(⟨ P⟩)/4π G M_*(⟨ P⟩) g,




where L(⟨ P⟩) is the total luminosity on the isobar and
M_*(⟨ P⟩) is the total mass integrated to the isobar. Approximating the effective gravity and the
shift in the mass as in <cit.> (Equation 2.3), this becomes





  
    F_V ≈L(⟨ P⟩) (g_0-Ω_0^2rsinθ)/4π G M[1-Ω_0^2/(2π G ⟨ρ⟩)],




where r and θ correspond to the origin of the local coordinate system as in Figure
<ref>. Taking the ratio of the rotating and non-rotating fluxes, one has that





  
    F_V/F_V,0≈L(⟨ P⟩)/L(r)1-Ω_0^2rsinθ/g_0/1-Ω_0^2/(2π G⟨ρ⟩),




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





  
    F_V/F_V,0≈L(⟨ P⟩)/L(r)1-3N_*,0^2rsinθ/(50π^2 Ro_c^2 g_0)/1-3 N_*,0^2/(100π^3 Ro_c^2 G⟨ρ⟩).




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





  
    F_V/F_V,0≈L(⟨ P⟩)/L(r)1-3|α_TΔ T| sinθ/(50π^2 Ro_c^2 )/1-|α_TΔ T|g_0
        r^2/(25π^2Ro_c^2 G M)
       ≈L(⟨ P⟩)/L(r)[1+|α_T Δ
          T|/25π^2Ro_c^2(g_0r^2/GM-3/2sinθ)] +⋯.




Therefore, since the linearized change in flux is proportional to the Boussinesq expansion parameter
|α_T Δ T| and the inverse square of the convective Rossby number, this requires
Ro_c > 5^-1π^-1 |α_T Δ 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, 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 ŝ 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 (e.g., V_0≪ 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





  
    max[F]/max[F]_0   =25/6√(5/3)[ŝ^3/q^3z^3+V_0ŝ^2/q^2]_max
       ≈.25/6√(5/3)ŝ^3/q^3z^3|_max =1,




implying that





  
    ŝ=s̃ q z + 𝒪(V_0),




where s̃=2^1/3 3^1/2 5^-5/6 and max[F]_0=6/25√(3/5) F_0 follows from the
definition of the flux and the maximizing wavevector used to define _0 above in Equation <ref>.

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






  
    /_0   = k_0/s_0s/k = 5/√(6)N_*/N_*,0ŝ/z^3/2 =
                              5/√(6)ŝ/q z^3/2 = (5/2)^1/6 z^-1/2.




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 <ref>.





 §.§ Diffusion-free Models


As considered in , the simplest case to consider is one in which all diffusive mechanisms are
neglected. Given Equation <ref>, the maximization condition requires





  
    dŝ/dz   = ŝ/z,




where the derivative of the growth rate can be determined from the nondiffusive limit of Equation <ref> as 





  
    dŝ/dz   = -3z^2ŝ(ŝ^2-1)/1+O_0^2 q^2cos^2θ+z^3(3ŝ^2-1);




when equated this yields the following equation for the growth rate





  
    ŝ^2 = 4z^3-O_0^2q^2cos^2θ-1/6z^3.




However, the dispersion relationship (Equation <ref>) for the nondiffusive case requires






  
    ŝ^2 = z^3-O_0^2q^2cos^2θ-1/z^3.



    
Using the heat flux invariance condition (Equation <ref>) so that
ŝ=s̃ q z, 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





  
    2z^5-5z^2-3/s̃^2 O_0^2cos^2θ = 0,
       q^2 - z^3-1/s̃^2z^5+O_0^2cos^2θ=0.




Together with Equations <ref>, these two equations yield the following scalings with
convective Rossby number and maximizing wavevector as





  
    ϵ/ϵ_0   = 25π^2_c^2s̃^2z^5+6cos^2θ/25π^2_c^2(z^3-1),
    
        2   z^5-5z^2-18cos^2θ/25π^2_c^2s̃^2=0.



    
With the quintic form of Equation <ref>, 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 _c≪1 it
is clear that z ∝_c^-2/5 since z^2≪ z^5. It then follows that
ϵ/ϵ_0 ∝_c^-4/5 and /_0 ∝_c^1/5, all of which
match the scalings given in , the direct numerical simulations of <cit.> and
<cit.>, and can be seen in Figure <ref>(a).  In contrast, at very large convective Rossby
number, these quantities converge to unity as expected.



  



At a first glance, one can see from Equations <ref> and <ref> 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θ, 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
minimum at the equator.  Hence, the heat carrying motions are of a larger scale at larger colatitudes than at the pole
as expected in theoretical considerations <cit.> and as seen in numerical simulations of
spherical rotating convection in both fully convective domains and in spherical shells
<cit.>. These behaviors are illustrated in Figure
<ref>(b).



 §.§ Adding Thermal Diffusion


The procedure developed above can be applied again in the case with thermal diffusion as the heat flux maximization
condition remains unchanged from Equation <ref>. Taking the limit that V→ 0 in Equation
<ref>, one finds that 





  
    dŝ/dz   =-3z^2(q^3 K_0 O_0^2cos^2θ +ŝ(ŝ^2+2 K_0 q z^3ŝ-1))/1+q^2 O_0^2cos^2θ+z^3(3ŝ^2+2K_0 q z^3ŝ-1),




which implies the following heat flux maximization condition





  
    ŝ(1+q^2O_0^2cos^2θ+2z^3(3ŝ^2+4q K_0z^3ŝ-2)) 
           +3q^3 K_0 O_0^2 z^3cos^2θ =0.




The dispersion relationship (Equation <ref>) with V_0→ 0 is





  
    (ŝ+ K_0 q z^3)(ŝ^2 z^3+ O_0^2q^2cos^2θ) -ŝ(z^3- 1)=0.




Again utilizing the maximized heat flux constraint ŝ=s̃qz and the definitions
in Equations <ref>, the two previous equations can be shown to be equivalent to





  
    q^2 - s̃(z^3-1)/(s̃ + K_0 z^2)( s̃^2z^5 +O_0^2cos^2θ)= 0,
       s̃z(2z^5- 5z^2-3/s̃^2 O_0^2cos^2θ)
           + K_0[4z^8- 7z^5-O_0^2/s̃^2cos^2θ(z^3+ 2)]=0.




So the connection to the nondiffusive case is clear and can be recovered directly when taking the limit as
K_0→ 0.  Employing the relationship between O and the convective Rossby number given in Equation
<ref>, the superadiabaticity and wavevector can be defined as





  
    ϵ/ϵ_0 = (s̃ + K_0 z^2)( 25 π^2 _c^2 s̃^2z^5 +6 cos^2θ)/25 π^2 _c^2 s̃(z^3-1),
       s̃ z[2z^5- 5z^2- 18cos^2θ/25π^2_c^2s̃^2]
           +K_0[4z^8-7z^5- 6cos^2θ(z^3+ 2)/25π^2_c^2s̃^2]= 0.




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 <ref> indicates that the wavevector must increase with decreasing
convective Rossby number, whereas its behavior with varying the thermal diffusion rate (K_0) is not obvious and is illustrated in
Figure <ref>(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_0.  For K_0>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
<ref>(a). The velocity decreases with lower convective Rossby number given that the wavevector increases,
with the same nontrivial behavior with the thermal diffusivity.





 §.§ Including Viscosity


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 ŝ and from





  
    F/F_0    = s̃^3[ŝ/q^3z^3+V_0 q]ŝ^2=1.




In particular, the above equation, Equation <ref>, and Equation <ref> give





  
    ŝ^3+(V_0 qŝ^2-s̃^3) q^3 z^3=0,
       (ŝ+K_0 q z^3)[z^3(ŝ+V_0 q z^3)^2+6cos^2θ/25π^2_c^2]
           +(1-z^3)(ŝ+V_0 q z^3)=0, 
       6q^2cos^2θ/25π^2_c^2[ŝ^2+K_0 q z^3(3ŝ+2 q V_0 z^3)]
           +(ŝ+V_0 q z^3)[6ŝ^3z^3 + 4 q ŝ^2 z^6(3V_0 + 2K_0).
            +ŝ(1-4z^3+2V_0 q^2 z^9(8K_0+3V_0))
           .+2V_0 q z^3(1-2z^3+4K_0V_0 q^2 z^9)]=0.




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 neglected if it is again assumed that V_0≪ N_*,0), so that as above
ŝ≈s̃ q z + 𝒪(V_0). Subsequently, following the method shown for the case
of thermal diffusion, one may find the implicit wavevector derivative of the growth rate ŝ (Equation
<ref>), which implies that the constraining dispersion relationship and the equation resulting from the heat
flux maximization are





  
    (s̃+V_0 z^2)[q^2z^5(s̃+V_0z^2)(s̃+K_0z^2)-z^3+1]
       +6 cos^2θ q^2 (K_0 z^2+s̃)/25 π ^2 _c^2=0,
       s̃(1-4z^3)+3V_0 z^2(1-2z^3)+6cos^2θ(s̃+3K_0 z^2)/25π^2 _c^2
       +2z^5(s̃+V_0z^2)[3s̃^2+s̃(4K_0+5V_0)z^2+6K_0V_0z^4]




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





  
    z^3(V_0 z^2+s̃)^2[3V_0 K_0 z^4(2 z^3- 3).
               .+s̃ z^2(V_0+ K_0)(4z^3- 7)+s̃^2(2 z^3- 5)]
       -6 cos^2θ/25 π^2 _c^2[2s̃(K_0-V_0)+3s̃^2z.
               .+s̃(K_0+5V_0)z^3+3K_0V_0z^5]=0.




The horizontal wavevector is then recovered from the roots of the fourteenth-order Equation
<ref>, whereas the superadiabaticity is defined as





  
    ϵ/ϵ_0=(s̃+K_0 z^2)(25 π^2 _c^2
          s̃^2z^5(s̃+ V_0 z^2)^2+ 6 cos^2θ)/25π^2 _c^2 s̃(z^3-
          1)(s̃+ V_0 z^2).




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 <ref>,
the wavevector must increase with decreasing convective Rossby number, whereas again its behavior with varying K_0 is not
obvious, but is relatively benign as is shown in Figure <ref>(a). In contrast, viscosity has a much larger
impact on the amplitude of the change in the wavevector as exhibited in Figure <ref>(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.  

However, for values of {K_0, V_0}<N_*,0, one may find a heuristic fit to the horizontal wavevector for the diffusive case
given by





  
    z ≈ 2^-11/6 5^1/3[3_c^2√(2)/5]^√(2)/2tanh[17/10_c^23/124 V_0^-37/160]
       [2+(√(2)-1)(1+tanh[√(2)/2(3/2+log_10_c/(V_0+
          K_0)^7/6)])],




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 <ref>, <ref>, and
<ref>, where it is well-described by hyperbolic tangents to within a few percent even in the transitional
region of _c≈ 1/10. The accuracy of the fit is illustrated in Figure <ref>(b).





 §.§ Diffusion Approximation for Turbulent Transport


In stellar models, a diffusive treatment of transport processes can be adopted for thermodynamics, chemicals, angular
momentum, and magnetic field <cit.>. Within a convection zone, the turbulent
mixing of those quantities can be approximated through a parameterized vertical diffusion
D_v = 1/3  α H_P v_MLT. Since this depends upon 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 α_T T ϵ/H_P. Thus, allowing H_P to vary, one has that
N^2/N_0^2 = x H_P,0/H_P = x/h. To maintain the constancy of the heat flux, this then requires
x=ϵ/ϵ_0 to depend linearly on h, with x→ xh.  

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





  
    D_/D_0= h/_0 = (5/2)^1/6h/√(z),




where D_0 is the local diffusion coefficient in the absence of rotation. The scaling of D_ is shown in
Figure <ref> and described in detail in the following section, as it suffices to note that the diffusion
will increase with 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.



§ CONVECTIVE PENETRATION




 §.§ Models Excluding Rotation


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
<cit.>. There the convective heat flux is 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 ℓ_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 upon the energy flux in the star <cit.>. Other models have also been developed that
variously extend MLT to include nonlinear processes through finite amplitude analyses
<cit.>. Subsequently, modal expansions of the Boussinesq and anelastic equations provided
solutions that established that nonlinear penetration was substantially deeper than what linear theory predicted
<cit.>. In anelastic convection, the density stratification induces asymmetries between
upflows and downflows that further increases the penetration depth <cit.>.  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 <cit.>.

Three pioneering studies analytically estimate 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 <cit.>, 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 <cit.>, <cit.> uses analytical
models of turbulent plume ensembles and examines their statistical properties, finding that the penetration depth scales
as _p^3/2 f^1/2, where _p is the plume velocity and f is the fractional area filled by the
plumes. Likewise, following <cit.> and <cit.> and by directly parameterizing the filling factor
f, <cit.> found a similar scaling result through a spatial linearization of the equations of
motion and the energy transport. Both of these studies predict penetration depths as fairly large fractions of a
pressure scale height.  also confirmed the modal results of <cit.> 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
 linearized model for convective penetration. This model neglects 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  model derived
above provides a means of estimating the latter alteration. In particular, one can harness the linearized
framework of the  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.



 §.§ Convective Penetration with Rotation




  §.§.§ Theory


The  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
(x,y,z) = (z) h(x,y) such that f=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 <cit.>.  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
<cit.>. Yet, 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 α_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





  
    F_conv   = ρ c_P  T' = fρ c_P δ T 
       = F_0[dlnM/dlnr-dlnL/dlnr+dlnχ/dlnr]_r_0z/r_0,




where M(r) is the mass, L(r) is the luminosity, and χ(r) is the radiative conductivity such that
χ = ρ c_P κ as seen in  (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 ρ' = α_T T'. To follow  as closely as
possible, the system is considered only at the pole so that the direct effects of the local Coriolis acceleration
2Ω_0sinθv_x may be neglected, 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  (Equation 3.8), the force balance becomes





  
    c/2d^2/dz = g δ T/T,




where c=h^3/h^2.  This equation can then be linked to the convective energy flux
through Equation <ref> 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 , 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 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 takes place near the convective cores of intermediate mass and high mass
stars.  The first of these two regimes neglects the variation in the total luminosity and mass and so only retains the
term related to changes in the radiative conductivity. It scales as





  
    L_P/H_P=[2/3(1-f) f _z^3/g α_T κχ_P∇_ad]^1/2,




where ∇_ad is the adiabatic temperature gradient and
χ_P=∂lnκ/∂lnP|_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





  
    L_P/H_P=[2(1-f) f _z^3 r_0/9g α_T κ∇_ad H_P(ρ P/ρ P - ρ Q/ρ Q-r_0χ_P/3H_P)]^1/2,




where r_0 is the radius at the edge of the core. Note that ∇_ad = ∇ + ϵ, 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





  
    L_P/L_P,0=(/_0)^3/2 = (5/2)^1/4z^-3/4.






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 <ref>). 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 a larger thermal conductivity.  In the case of a 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 <ref>.



  §.§.§ Comparison with Penetrative Convection Experiments


Penetrative convection has been studied extensively in fully compressible and anelastic numerical experiments in 2D
within the context of the solar convection zone <cit.>, for the convective
cores of intermediate mass stars <cit.>, and in the context of modeling laboratory experiments for
temperature stratified water <cit.>. These studies all tend to
find that penetration into stable layers below the convection zone is more extensive than into overlying stable layers
and that the depth of penetration depends sensitively on the stiffness of the interface between them and upon the
strength of the convective driving.  Moreover, in rotating convection models of equatorial planes of intermediate mass
stars, it is directly seen that rotation reduces the overshooting depth as the Coriolis force leads to an azimuthal
deflection of convective plumes and to less available radial kinetic energy <cit.>. While 2D simulations are a
useful first step to understanding convective penetration, they are limited in assessing the convective properties of
stellar interiors given the 3D nature of the convection there, which leads to much different convective structures and
spectral energy transfer properties.

In constrast to the behavior of penetrative convection seen in local simulations and from scaling arguments such as
those described above, the depth of penetration in fully nonlinear global-scale 3D simulations tends to increase with
decreasing convective Rossby number and in spherical geometry
<cit.>.  However,
these simulations are typically in a low Péclet number regime Pe ≈𝒪(10) and thus are still
influenced by thermal diffusion. The stiffness of the interface often has been a priori softened so that the simulation
is more computationally feasible, meaning that the degree of convective overshoot relative to penetration is larger than
in a stellar context and that the excitation mechanisms are also more dominated by Reynolds stresses rather than
buoyantly driven.  So, it may be that the stiffness of the transition region is too weak in current global-scale
modelling approaches, permitting large-scale entraining and overshooting flows rather than restricting them to
plume-like dynamics.  Moreover, the sweeping motions of the relatively laminar globally connected flows induce mixing on
a scale of the correlation length of those flows <cit.>, which enhances the depth of penetration and is
something that will be assessed in later work.

Three-dimensional simulations have been conducted that can simultaneously resolve gravity waves as well as the
mechanisms that excite them, namely shear and convective penetration in both nonrotating Cartesian settings
<cit.> and f-planes <cit.>.  As in
2D simulations, these 3D simulations also find the following: penetration into stable layers below the convection zone
is more extensive than into overlying stable layers, the depth of penetration depends sensitively on the stiffness of
the interface between them, and the penetration depth is latitudinally dependent.

In the 3D f-plane simulations of rotating convection described in <cit.> and <cit.>, it is found
that the penetration depth into a stable layer below a convective region scales as
L_P ∝_c^0.15, that is it decreases with increasing rotational influence, due primarily to a
reduction in the flow amplitude.  Likewise, in a similar suite of f-plane simulations examined in <cit.>, it is
found that there is a decrease in the penetration depth with increasing rotation rate that scales as
L_P ∝_c^0.2 at the pole and to L_P ∝_c^0.4 at mid-latitude.  For
penetration into a stable layer above a convection zone, on the other hand, f-plane simulations from <cit.>
indicate that there is a much weaker rotational scaling, being statistically consistent with no scaling. However, the
range of parameters examined is quite restricted due to the computational requirements to resolve structures as well as
maintain highly supercritical flows at much lower convective Rossby numbers. Nevertheless, the depth of convective
penetration as assessed in those numerical simulations appears to be roughly consistent with the heuristic model derived
above, where L_P/L_P,0∝_c^3/10, which follows from
/_0 ∝_c^1/5 in the nondiffusive and low convective Rossby number limit of the
convection model.

<cit.> have examined the influence of the stiffness of the convective-radiative transition in numerical
analogs of a laboratory experiments involving the buoyancy transition of water that occurs near 4^∘ C, which is
also a decent analog for convective stellar cores.  In these 2D simulations, there is a smooth transition between plume
dominated and entrainment dominated dynamical regimes that is independent of the Peclet number.  Instead, it is
sensitive to the stiffness of the convective-radiative transition region and the Rayleigh number.  According to those
local domain numerical simulations, and to the theory developed in <cit.>, the interface stiffness provides
a further distinction drawn between the regimes of penetrative convection. Particularly, the penetration depth depends
upon the relative stability S=(2 n_s-3)/(3-2 n_c) between the stably-stratified and convective regions, where n_s is
the polytropic index in the stable region and n_c is the index in the convective region.  Indeed, <cit.>
employs a linearized convection model similar to the one derived in , which has a continuous thermal
conductivity, but instead the conductivity is taken to be piecewise constant. Using this model, it can be shown that the
total penetration depth depends upon the depth of the adiabatic layer and the depth of the thermal adjustment layer.
However, the depth of the two regions scale differently with the relative stability parameter S.  In particular, the
depth of the adiabatic layer scales as L_a ∝ S^-1 and the thermal adjustment layer scales as
L_t ∝ S^-1/4.  Therefore, as S increases the depth of the adiabatic layer decreases significantly, leaving
the thermal adjustment layer as the primary means of defining the depth of penetration.  Note that as stated in
<cit.> (Equation 3.18), the penetration model developed there and expanded on above scales proportionally to the
inverse of the interface stiffness. For the simulations of penetration into a stably-stratified layer below a convection
zone, <cit.> and <cit.> find that the depth of penetration scales primarily as
L_P ∝ S^-1/4, for values of S between about unity and around 30.  Yet for simulations of convective
penetration into a stably-stratified region above a convection zone, the opposite behavior is found, namely that
L_P ∝ S^-1 for S between unity and ten. Given the similarity in the parameters used in those simulations,
one might conclude that the penetration into a stably-stratified layer below a convection zone has physics more akin to
a thermal adjustment layer, whereas penetration into a layer above is more akin to an adiabatic penetrative layer.





 §.§ A Diffusive Approach


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
<cit.>. A complementary model has been established
through an extreme-value statistical analysis of 3D penetrative convection simulations <cit.>.  A frequent
assumption made of turbulent flows is to assume that the underlying distribution of a flow quantity is normal. However,
it is shown through the simulations and analysis of <cit.> and <cit.> that the more rare events
occurring in the tails of the actual distribution function, e.g. 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 a turbulent diffusion based upon the Gumbel distribution <cit.>. This model
for diffusion and overshoot has been applied to the problem of lithium depletion occurring in rotating low-mass stars
<cit.>, 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
<cit.> may be improved. Using the above extension of the  model, one can estimate both the
penetration depth and the level of diffusion. Taking the parameters of the Gumbel distribution as in <cit.>,
yields the following description of the radial dependence of the diffusion coefficient





  
    D_(r) =(5/2)^1/6α H_P_0 h/3√(z)
           {1-exp[-exp((r-r_c)/λ L_P+μ/λ)]},




where r_c is the base of the convection zone and where μ and λ are the empirically determined
parameters from <cit.>.  An illustration of the scaling behavior due both to the variation of D_ with
the mixing length velocity  and with the depth of penetration L_P are shown in Figure <ref>, where
the parameters are taken to make contact with Sun-like stars where the transition region begins around
r≈ 0.7 R_. Note, however, the convective velocity will vary in the convection zone. So Figure
<ref> is only meant to illustrate the rough behaviors of this diffusive model near the convectively
stable-to-unstable transition point. From Equation <ref> and as seen in Figure <ref>, 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.

In <cit.>, 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 is evolved that permits the
measurement and fitting of a height-dependent turbulent chemical diffusivity that mimics the effect of the overshooting
convection.  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 <cit.> have also addressed the
shape of the overshooting region in a Boussinesq system and its associated potential mixing properties in 3D
simulations.  As in <cit.>, 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
modelled with a Gaussian. Moreover, the radial correlation length scale downflows is larger than the length scale
associated with the averaged kinetic energy profile. This length scale may be interpreted as 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 <cit.> and <cit.>, although it is in contrast to the models proposed in
<cit.> and <cit.>. For the time being, the examinations in <cit.> and <cit.>
will have to suffice as evidence that the overshooting is reduced with rotation rate, while the shape of the mixing
region likely depend upon the model setup and equations solved in similar simulations.



§ SUMMARY AND DISCUSSION


A model of rotating convection originating with <cit.> has been extended to include thermal and viscous
diffusion for any convective Rossby number.  Moreover, a systematic means of developing such models for an arbitrary
dispersion relationship have also been shown.  An explicit expression is given for the scaling of the horizontal wavenumber
in terms of the convective Rossby number and diffusion coefficients under the constraint that the values of the
diffusive time scale are less than the convective time scale (Equations <ref>, <ref>, and
<ref>).  The scalings of the velocity and superadiabaticity in terms of that wavenumber are also given.
Asymptotically at low convective Rossby number and without diffusion, these match the expressions given in
<cit.> (see Figure <ref>), as well as the numerical results found in the 3D simulations
of <cit.> and <cit.>.

In <ref>, the model of rotating convection is employed to assess the convective Rossby number scaling
of the depth of convective penetration, utilizing the linearized model of <cit.>. Due to the reduced velocity
and increased superadiabaticity at lower convective Rossby number, the penetration depth decreases proportionally to
_c^3/10 when diffusive processes are neglected. This estimate of the penetration depth is then
employed to construct an estimate for the diffusive mixing coefficient in the region of penetration based upon the
numerical experiments of <cit.> and <cit.>. In the 3D f-plane simulations of rotating convection of
<cit.> and <cit.> it is found that L_P∝_c^0.15. Similarly, in the
simulations of <cit.> and <cit.> it is found that L_P varies in latitude, with
L_P∝_c^0.2-0.4 from the pole to mid-latitudes. Hence, the depth of convective penetration as
assessed in those numerical simulations is roughly consistent with those that follow from the coupling of the convection
model with the results of <cit.>.

This convection model will soon be exploited to begin to study the impact of reduced convective penetration in stellar
and planetary evolution models when rotation is included. Yet its impact on the structure of the convection can already
be anticipated. From Figure (<ref>), since this convection model is local, a qualitative picture
of those aforementioned impacts can be constructed as follows: in regions where the local convective Rossby number is
less than unity, the mixing-length velocity will be reduced and the superadiabatic temperature gradient will be
enhanced.  Therefore, depending upon the value of the superadiabatic gradient in the nonrotating case, the full
temperature gradient may increase sufficiently to modify the location of the ionization zones as well as regions of
large opacity changes.  However, most of those thermodynamic changes will be felt most keenly deeper within the
convection zone where the convective Rossby number can be smaller than near the photosphere for stars, where the
convective Rossby number is typically of order unity or larger.  Near the boundary of convective regions the typical
mixing-length velocity is small due to the increasing importance of the radiative transport of energy.  Thus, in those
regions, the local convective Rossby number can be quite small and yet the thermal diffusivity can become larger
relative to the dynamical time scale.  With that in mind, and again appealing to Figure (<ref>),
it is clear that the velocity will be even further reduced and the superadiabatic gradient further increased. This
implies that the transition to the region of convective penetration and the radiative zone itself will be shallower (or
deeper for a convective core) when compared to the nonrotating case.  Hence, in addition to the reduced depth of
convective penetration with decreasing convective Rossby number seen in <ref>, the convection zone
depth will be reduced when considering convection in the presence of rotation with this model. Due to decreased
diffusion and transport, sharper thermodynamic and chemical gradients may be present in convectively stable regions.

Being a first step, there are several improvements that can be made to this convection model. One clear omission in this
work is the impact of the magnetic field. However, magnetism has been considered in both <cit.> and
<cit.>. An alternative model for the impact of the magnetic field on the superadiabaticity, derived from the
variational principle of <cit.>, can be found in <cit.>. Thus, convection and magnetic field will
be the focus of a paper in this series.  Second, the conjecture of <cit.> gives the scaling of only the
dominant heat-carrying mode.  To build a theory of the turbulent spectrum and to assess the changing cascade of energy
and helicity, one could extend the theory in a manner similar to that shown in <cit.>, with a mode by mode
analysis of the heat flux to construct its accompanying spectrum.  One attempt at generalizing this turbulence model has
already been undertaken utilizing the Heisenberg-Kolmogorov turbulence model, as shown in <cit.>. It would
also be prudent to further numerically assess the validity of these scaling laws with 3D simulations with a larger range
of Nusselt number and to examine the impact of diffusion. The theory can be further improved by considering more
sophisticated models of the structure of the flows, such as applying the results obtained for rotating plumes considered
in <cit.> or <cit.>.  Additionally, one relatively simple improvement to the treatment of convective
penetration will be to include the effects of rotation in the linearized equations of motion in an f-plane, expanding
further upon the <cit.> model.

These waves, which are internal gravity waves modified by rotation through the Coriolis acceleration
<cit.>, propagate in stably stratified stellar radiation zones. They provide a means for
transporting angular momentum and mixing chemical species in these regions
<cit.>. They are one of the mechanisms that may explain the strong
extraction of angular momentum in the Sun <cit.>, in subgiant and red giant stars
<cit.>, and in intermediate-mass and massive stars
<cit.> as revealed by the weak differential rotation observed thanks to helio- and asteroseismology
<cit.>. However, their amplitude
and frequency spectrum strongly depend on the properties of the turbulent convective flows that excite them
stochastically at the radiation/convection interface and in the bulk of convective regions
<cit.>. Therefore, one should
consider the impact of the modification of their convective excitation source by rotation. Although some first attempts
have been made to do this <cit.>, a systematic study should be done. This will be adressed
in the second article of this series.

Another possible use of the convection model developed here is for tidal dissipation. Turbulent friction is one of the
crucial physical mechanisms driving the dissipation of tidal flows in stellar and planetary convective regions
<cit.>. This friction acts both on the equilibrium tide and on tidal inertial waves in
these layers. Recently, in <cit.>, the asymptotic version of Stevenson's theoretical scaling laws was used to
construct a corresponding local model of tidal waves in order to understand the consequences of rotating convective
turbulence for linear tidal dissipation. With this convective model, the turbulent friction acting on the tides was
significantly reduced in rapidly rotating objects. Thus, to better capture the behavior of the tidal friction in the
convective regions in stars and planets at modest convective Rossby number, it will be useful to apply the model of
convection developed here.



§ ACKNOWLEDGMENTS
 The authors thank the anonymous referee for their constructive comments, which have improved
  the description of the convection model and its implications. K. C. Augustson and S. Mathis acknowledge support from
  the ERC SPIRE 647383 grant and PLATO CNES grant at CEA/DAp-AIM.  The authors also thank Q. André, A. Astoul,
  M. Browning, A. S. Brun, V. Prat, R. Raynaud, and J. Toomre for fruitful conversations.



