THE MECHANICAL GREENHOUSE: BURIAL OF HEAT BY TURBULENCE IN HOT JUPITER ATMOSPHERES

Our goal is to understand energy balance. We focus on the deep atmosphere which is optically thick both to incoming stellar irradiation and the planet's emitted flux. Here, the equation of radiative diffusion

$$\begin{equation} {dT \over dP} = { F_{\rm rad}\over k_{\rm rad}} \end{equation} \tag{ 1 }$$

relates the outgoing radiative flux F_rad to the variation of temperature T with pressure P via the radiative diffusion coefficient

$$\begin{equation} k_{\rm rad} = {16 \sigma T^3 g \over 3 \kappa }, \end{equation} \tag{ 2 }$$

where σ is the Stephan–Boltzmann constant. Hydrostatic balance, dP/dz = −ρg, allows pressure to replace height z as the vertical coordinate, with ρ the atmospheric density. We hold gravity g constant, invoking the plane-parallel approximation for thin atmospheres.

For the Rosseland mean opacity, our calculations will use a power law,

$$\begin{equation} \kappa = \kappa _1 P^\alpha T^{\beta } \equiv \kappa _o \left(P \over P_{\rm kb}\right)^\alpha \left(T \over T_{\rm 2k}\right)^{\beta } \,. \end{equation} \tag{ 3 }$$

The two forms are equivalent, with the constant κ₁ being more compact, while κ_o has units of opacity and is normalized to P_kb = 1 kbar and T_2k = 2000 K. Unless stated otherwise, calculations will use κ ∝ P, i.e., α = 1, β = 0 as a rough approximation to collision induced molecular opacity. Our normalization choice of κ_o = 0.18 cm² g^-1 will be justified below (Section 2.2). We discuss alternate opacity laws in Section 4.4. While realistic opacities are only well approximated by power laws over a limited range, this considerable simplification is useful for developing intuition.

At the top of our atmosphere, we set the temperature to T_deep, an approach used in AB06 and advocated by Iro et al. (2005). This approach is valid when the incident stellar flux exceeds the emitted radiation, resulting in a deep isothermal region at the top of the optically thick atmosphere. In this physical situation, the precise location of the upper boundary is not important, as we explain further in Section 2.2. We point the reader to Hansen (2008) and Guillot (2010) for sophisticated analytic treatments of RE, which are very useful for interpreting computational models. Note that some (Baraffe et al. 2003) but not all (Seager & Sasselov 2000) detailed radiative transfer solutions show an extended isotherm in hot Jupiter atmospheres. One difference between models is the choice of non-gray opacities, which affects the depths at which starlight in different frequency intervals is absorbed.

The incident stellar flux averaged over the full planetary surface, F_irr ≡ σT⁴_*, gives a characteristic temperature

$$\begin{equation} T_\ast \approx 2000\, {\rm K} {L_{\ast,\odot }^{1/4} \over M_{\ast,\odot }^{1/6} P_{{\rm day}}^{1/3}} \end{equation} \tag{ 4 }$$

where stellar mass and luminosity, M_*,☉ and L_*,☉ are normalized to solar values, and the orbital period P_day is normalized to a (24 hr) day. Horizontal temperature gradients are important for driving winds in the weather layer. However, these winds efficiently smooth temperature gradients at pressures ≳ bar, where timescales for advection are shorter than for radiative losses (Showman et al. 2009a). Thus, one-dimensional models are appropriate for basic considerations of energy balance.

Because the planet is not a perfect blackbody, T_deep may not match T_*. Greenhouse or anti-greenhouse effects depend on the relative transparency of the atmosphere to stellar and emitted longwave radiation. (To be clear, we are now referring to standard radiative effects, not the mechanical greenhouse.) If incoming starlight penetrates below the infrared photosphere, then the greenhouse effect gives T_deep > T_*. If, however, significant incoming radiation is absorbed above the photosphere, a stratospheric thermal inversion gives T_deep < T_*. See Hubeny et al. (2003) for a more quantitative analysis. In most of our examples, we adopt T_deep = 1500 K, as appropriate for short (∼1 day) orbital periods with a thermal inversion, or for longer periods with no thermal inversion, a greenhouse effect, and/or a more luminous host star.

Giant planets, including hot Jupiters, become unstable to convection at depth. The lapse rate of the atmosphere

$$\begin{equation} \nabla \equiv {d \ln T \over d \ln P} = {3 \kappa P \over 16 g} {F_{\rm rad}\over \sigma T^4} \end{equation} \tag{ 5 }$$

characterizes its stability, and the final equality follows from Equation (1). Convection occurs where ∇>∇_ad. We set ∇_ad = 2/7, the adiabatic index of an ideal diatomic gas. In reality, non-ideal interactions lower ∇_ad at the high pressures of exoplanet atmospheres and promote convection. At the top the atmosphere, where the optical depth τ = κP/g ≈ 1 and F_rad ≪ σT⁴_deep, Equation (5) shows that ∇ ≪ 1 and the atmosphere is indeed stable and nearly isothermal. For reasonable opacity choices, ∇ increases with depth and gives a transition to convection (even under the ideal gas approximation).

In convective regions, we set ∇ = ∇_ad, i.e., an adiabatic profile with $T \propto P^{\nabla _{\rm ad}}$. The efficiency of convective energy transport makes the modest super-adiabaticity negligible. The level of the adiabat is determined by the internal entropy. A global calculation of entropy is beyond our illustrative scope. Instead, motivated by Hubbard (1977), we label our adiabats by T₁, the temperature it would have at P₁ = 1 bar pressure, even though the adiabat likely does not extend to such a low pressure. We define a reference entropy (per unit mass) S_ref, corresponding to T₁ = 250 K. Relative entropy values for different T₁ are then computed as

$$\begin{equation} \Delta S \equiv S-S_{\rm ref} = C_P \ln (T_1/250\, {\rm K}), \end{equation} \tag{ 6 }$$

where the specific heat (at constant pressure) C_P = R/∇_ad is assumed constant. For the gas constant R = k_B/(μm_p), we use a mean molecular weight μ = 2.34 times the proton mass m_p.

The stable atmosphere matches smoothly onto the convective adiabat at the RCB. Since the temperature T_c and pressure P_c at the RCB lie on the interior adiabat we require

$$\begin{equation} T_1 = T_{\rm c} \left(P_1\over P_{\rm c}\right)^{\nabla _{\rm ad}}\,. \end{equation} \tag{ 7 }$$

The location of the RCB is crucial for global evolution. The secular cooling of the convective interior is determined by the radiative flux, F_c, leaving the RCB. Combining Equations (3), (5), and then (7) at the RCB gives

$$\begin{eqnarray} F_{\rm c} &=&{\nabla _{\rm ad}16 g \sigma T_{\rm c}^{4-\beta } \over 3 \kappa _1 P_{\rm c}^{1+\alpha }}\\ &=& {\nabla _{\rm ad}16 g \sigma \over 3 \kappa _1 } \left({T_1 \over P_1^{\nabla _{\rm ad}} P_{\rm c}^{\nabla _\infty - \nabla _{\rm ad}}}\right)^{4-\beta } \propto {T_1^4 \over P_{\rm c}^{6/7}}, \end{eqnarray} \tag{ 8 }$$

with

$$\begin{equation} \nabla _\infty \equiv (1+\alpha)/(4-\beta) = 1/2\,. \end{equation} \tag{ 10 }$$

Core flux increases with the interior entropy.⁵ Pushing the RCB to higher pressures decreases the core flux if ∇_∞ > ∇_ad. This condition is satisfied for our opacity choice and is generally required for a transition to convection (as we will show shortly). We emphasize that the dependence of F_c on P_c is independent of the mechanism that changes P_c, though previous works have mostly considered irradiation. These basic considerations are useful in interpreting numerical studies of planetary cooling histories and radii evolution (Burrows et al. 2000; Baraffe et al. 2003; Chabrier et al. 2004).

We now apply the two approximations of RE to the stable layer. First, radiation is the only relevant energy transport mechanism. Thus, F_rad in Equation (1) is the total flux of energy. Second, the flux is constant through the radiative layer with F_rad = F_c, the flux from the convective interior. This assumes that local (thermal) energy release is negligible.

Figure 2 plots radiative equilibrium atmospheres for κ ∝ P, with two values of T_deep matched on to interior adiabats labeled by T₁. We obtain analytic RE solutions by integrating Equation (1) with F_rad = F_c and T = T_deep at P = 0 to get

$$\begin{equation} T = T_{\rm deep}\left[1 + {\nabla _{\rm ad}\over \nabla _\infty - \nabla _{\rm ad}} \left({P \over P_{\rm c}}\right)^{1+\alpha }\right] ^{1 /(4-\beta)}\,. \end{equation} \tag{ 11 }$$

Zoom In Zoom Out Reset image size

This solution uses Equation (8) and we have imposed the requirement T(P_c) = T_c to find

$$\begin{equation} T_{\rm c} = T_{\rm deep} \left(\nabla _\infty \over \nabla _\infty - \nabla _{\rm ad}\right)^{1/(4-\beta)}. \end{equation} \tag{ 12 }$$

A valid solution—one that transitions to convection—thus requires ∇_∞ > ∇_ad and α> − 1 (which together assure β < 4). As Figure 2 shows, T_c increases with T_deep but is independent of the interior entropy.

The RCB sinks to larger pressure as entropy decreases or as T_deep increases,

$$\begin{equation} P_{\rm c} = k_{\nabla } P_1 \left(T_{\rm deep} \over T_1\right)^{1/\nabla _{\rm ad}}, \end{equation} \tag{ 13 }$$

which follows from Equations (7) and (12) with the constant

$$\begin{equation*} k_{\nabla } \equiv \left(\nabla _\infty \over \nabla _\infty - \nabla _{\rm ad}\right)^{\nabla _\infty /[\nabla _{\rm ad}(1+\alpha)]} \approx 2.1\, \nonumber \,. \end{equation*}$$

The core flux for RE atmospheres follows from Equations (8), (12), and (13) as

$$\begin{equation} F_{\rm c} = k_F {g \over \kappa _1} \left(T_1 \over T_{\rm deep}^{1- \nabla _{\rm ad}/\nabla _\infty } \right)^{(1 +\alpha) / \nabla _{\rm ad}} \propto {T_1^{7} \over T_{\rm deep}^3}\,. \end{equation} \tag{ 14 }$$

This gives the well-known result that increased irradiation reduces the cooling of the planet, while higher entropy planets are more luminous. The constant

$$\begin{equation*} k_F \equiv {16 \sigma \nabla _{\rm ad}\over 3 P_1^{1+\alpha }}\left(1 - {\nabla _{\rm ad}\over \nabla _\infty }\right)^{\nabla _\infty /\nabla _{\rm ad}-1}\,. \nonumber \end{equation*}$$

Equation (14) is consistent with, but more specific than, Equation (9) in assuming that RE sets the location of P_c.

We chose the parameters for Figure 2—used throughout this work—by roughly matching the analytic solutions to more detailed hot Jupiter models, as in AB06.

We constrain the entropy parameter T₁ by appealing to the typical P_c ≈ 1 kbar location of the RCB in hot Jupiters with modest, i.e., Jovian, entropies. With T₁ = 260 K, we reproduce a 1 kbar RCB for T_deep = 1500 K. We also consider larger values of T₁ to describe more inflated planets, but keep P_c ≫ 1 bar.

The normalization of the opacity determines the core flux.⁶ Requiring F_c = σ(100 K)⁴ for the standard parameters and P_c = 1 kbar gives κ_o = 0.18 cm² g^-1 for g = 10³ cm² s⁻¹. We emphasize that this is not a realistic opacity law (in particular it is too low at small pressures). We are merely choosing parameters that allow the simple analytic model to mimic properties of more detailed hot Jupiter models.

The lapse rate for RE solutions is (from Equation (11))

$$\begin{equation} \nabla = \nabla _{\rm ad}{ \left(P/P_{\rm c}\right) ^{1+\alpha } \over \left(1-{\nabla _{\rm ad}\over \nabla _\infty }\right) + {\nabla _{\rm ad}\over \nabla _\infty }\left(P/ P_{\rm c}\right)^{1+\alpha }}, \end{equation} \tag{ 15 }$$

demonstrating that ∇ = ∇_ad at P = P_c and that the solution becomes isothermal ∇ → 0 at low pressures. Smooth opacity laws give a monotonic increase in ∇ with P. Opacity windows give more complicated profiles of ∇, including multiple zones of convection (see Section 4.4).

We define P_deep, the effective depth of the isothermal layer, as the location where ∇ = ∇_ad/2. This occurs at

$$\begin{equation} P_{\rm deep} = \left({\nabla _\infty -\nabla _{\rm ad}\over 2 \nabla _\infty -\nabla _{\rm ad}}\right)^{1/(1+\alpha)}P_{\rm c} \approx 0.55 P_{\rm c} \,. \end{equation} \tag{ 16 }$$

Our definition of P_deep differs from AB06, who define P_deep as a characteristic scale that might exceed P_c.

We now revisit the validity of applying the boundary condition T = T_deep at P = 0. Due to the isothermal layer at low pressures, applying the boundary condition at any P ≪ P_deep gives indistinguishable solutions. However, solutions are only physically valid in optically thick regions, for P ≫ P_thick ∼ g/κ_min ∼ 10 mbar[κ_min/(0.1 cm² g^-1)]⁻¹. The relevant opacity κ_min is the smaller of the opacities to starlight and emitted radiation near P_thick. Indeed, the penetration of starlight below the infrared photosphere can push the top of the isotherm to ∼1 bar (Guillot 2010). As long as P_thick ≪ P_deep, solutions are physically consistent below P_thick.

We derive F_eddy using basic elements of mixing length theory. This theory is usually applied to convectively unstable regions, but instead we apply it to forced turbulence in convectively stable regions. We will show that in this case the energy flux is inward. We leave the forcing mechanism of turbulent motions unspecified and their strength a free parameter.

We consider parcels of gas which conserve entropy and maintain pressure equilibrium with their surroundings as they exchange positions over a vertical distance ℓ and then dissolve. These parcels contain excess heat δq = ρC_PδT with

$$\begin{eqnarray} \delta T &=& \left(\left. {dT \over dz}\right|_{\rm ad} - {dT \over dz}\right)\ell \nonumber \\ &=& {-}{\ell T \over C_P} {dS \over dz}, \end{eqnarray} \tag{ 19 }$$

where δf gives the difference of any quantity f between the parcel and its surroundings. For stable stratification (dS/dz > 0), rising parcels (ℓ>0) cool down and sinking parcels heat up. We express the heat flux, F_eddy = wδq with w the characteristic eddy speed, in terms of the turbulent diffusion, K_zz = wℓ, as

$$\begin{eqnarray} F_{\rm eddy}&=& {-} K_{zz} \rho T {dS \over dz} \nonumber \\ &=& {-} K_{zz} \rho g \left(1 - {\nabla \over \nabla _{\rm ad}}\right) \,. \end{eqnarray} \tag{ 20 }$$

The flux is always negative for stable stratification. It vanishes at the RCB where dS/dz = 0 (and ∇ = ∇_ad). We do not model overshoot, which could allow energy exchange (in either direction) between convectively stable and unstable zones. In the upper isothermal regions F_eddy ∝ −K_zzP, declines in magnitude with height, unless K_zz increases with height to compensate, as we will consider in Section 4.2.

The above assumption that parcels conserve entropy assumes a radiative cooling time longer than the eddy turnover time. We are ignoring radiative losses during eddy motion in this work, because we consider optically thick regions with long cooling times. This assumption will eventually break down in optically thin regions, notably the inversion layer itself, which is strongly stratified. Whether eddy fluxes affect the structure of the inversion layer is left to future work.

To understand the energetics of F_eddy, we analyze its divergence, which describes cooling and (when negative) heating,

$$\begin{equation} {d F_{\rm eddy}\over dz} = K_{zz} {\rho g \over R} {dS \over dz} - K_{zz} \rho T {d^2S \over dz^2} - {d K_{zz} \over dz} \rho T {dS \over dz}\,. \end{equation} \tag{ 21 }$$

The terms on the right-hand side represent cooling rates, $-\delta \dot{q}$, which we relate to rates of work, $\delta \dot{w}$, using the first law, δq = δe − δw. Thus, $-\delta \dot{q} = \delta \dot{w}/\nabla _{\rm ad}$ since the internal energy, δe = ρC_VδT = (1 − ∇_ad)δq with C_V = C_P − R. The first term in Equation (21) arises from buoyant work $\delta \dot{w}_{\rm B} = \delta \rho g w = K_{zz} \rho g (dS/dz)/C_P$, where δρ/ρ = −δT/T by pressure equilibrium. The first term on the right-hand side of Equation (21) is thus $-\delta \dot{q}_{\rm B} = \delta \dot{w}_{\rm B}/\nabla _{\rm ad}$. This buoyant cooling will be evident in the stratified regions (P < P_deep) of Figure 6.

The second term represents the tendency of mixing to heat by filling in entropy minima (for d²S/dz² > 0). More specifically, it arises from compressional work, which vanishes for a constant entropy gradient because the work done on rising and sinking parcels cancels out. For varying dS/dz consider two parcels that arrive at z, one from above and one from below. Compressional work is done on the parcels at a rate

$$\begin{equation} P \nabla \cdot \mbox{\boldmath {$v$}}_{\pm } = -{P \over \rho } {\delta \rho _\pm \over \delta t} = \mp {P \ell \over C_P} \left.{dS \over dz}\right|_{ z \pm \ell } {w \over \ell }, \end{equation} \tag{ 22 }$$

where the top (bottom) sign refers to sinking (rising) parcels, ∇ · v is a velocity divergence, and δt = ℓ/w gives the expansion rate. The net work is the sum of these terms, $\delta \dot{w}_{\rm C} = -\nabla _{\rm ad}K_{zz} \rho T d^2S/dz^2$. We identify the second right-hand side term in Equation (21) as $-\delta \dot{q}_{\rm C} = \delta \dot{w}_{\rm C}/\nabla _{\rm ad}$. This compressional heating dominates in deeper regions (P > P_deep) of Figure 6. The third and final term represents a flux imbalance that arises from non-uniform eddy diffusion as in Section 4.2.

Our atmospheric model is described by Equations (1), (17), (18), and (20) which reduce to the following pair of coupled ordinary differential equations (ODEs)

$$\begin{equation} {dT \over dP} = {F + F_{\rm iso} \over k_{\rm rad} + {F_{\rm iso} P /(\nabla _{\rm ad}T)}} \end{equation} \tag{ 23 }$$

$$\begin{equation} \hspace{-67pt}{d F \over dP} = -{\epsilon \over g} \end{equation} \tag{ 24 }$$

where

$$\begin{equation} \hspace{-60pt} F_{\rm iso} \equiv K_{zz} \rho g \end{equation} \tag{ 25 }$$

is (minus one times) the isothermal limit of the eddy flux. The thermal profile in Equation (23) describes the combined effects of radiative and eddy fluxes. Ignoring mixing (F_iso → 0) recovers standard radiative diffusion of Equation (1). Strong mixing (F_iso → ∞) creates an isentropic profile (∇ → ∇_ad).

Solving the coupled ODEs for T and F requires prescriptions for K_zz and and boundary conditions. As with the RE solutions, we fix T_deep at the top of the atmosphere and match onto an adiabat with a given T₁ at the bottom. This matching generally requires iterative techniques. We avoid this complication by finding self-similar solutions. This technique is only possible because of the idealizations—notably power-law opacities and ideal gas EOS—described in Section 2.1.

To find self-similar solutions, we normalize all quantities to their value at the RCB. As with the analytic RE solutions, we do not know where the RCB is located until we find the solution. We integrate outward from the RCB with trivial boundary conditions: the dimensionless T and F are unity.

We set the strength of diffusion not with a physical value for K_zz, but by the parameter

$$\begin{equation} \psi _{\rm c} = {K_{zz} \rho _{\rm c} g\over F_{\rm c} }, \end{equation} \tag{ 26 }$$

the ratio of F_iso to F = F_rad at the RCB. For dissipation, we must similarly specify P/(gF) at the RCB. See Appendix B for details.

To get a physical solution, we scale a dimensionless solution to any desired value of T_deep and T₁. Setting T = T_deep at P = 0, the solution gives T_c at the RCB. Specification of T₁ then fixes P_c via Equation (7). We then determine F_c and K_zz via Equations (8) and (26).

Turbulence that gives rise to diffusion, K_zz, will also dissipate at some rate . We now consider a prescription that sets a lower bound on from turbulence. We also allow for stronger heating, perhaps from non-turbulent sources.

In a Kolmogorov cascade, the dissipation rate = w³/ℓ and the diffusion K_zz = wℓ give a simple relation = K_zz/t²_o, with t_o = ℓ/w the turnover time of the dominant eddies. Unfortunately, we lack a reliable model for eddy timescales. Moreover, turbulence in a stratified atmosphere is likely anisotropic and not well described by Kolmogorov scalings. Fortunately, our results are not very sensitive to this limitation, as we will show in subsequent sections. Eddies with long turnover times, t_o > 1/N, organize into horizontally extended pancakes, where the squared buoyancy frequency is

$$\begin{equation} N^2 = {g^2\over R T}\left[\nabla _{\rm ad}-\nabla \right]\,. \end{equation} \tag{ 27 }$$

Assuming that the buoyancy frequency sets the relevant timescale, t_o = 1/N, gives a dissipation rate

$$\begin{equation} \epsilon _{\rm buoy} \approx K_{zz}N^2\,. \end{equation} \tag{ 28 }$$

In strongly stratified, isothermal regions the buoyancy is

$$\begin{equation} N_{\rm deep} = g/\sqrt{C_P T_{\rm deep}} \sim c_{\rm s}/H \end{equation} \tag{ 29 }$$

with c_s the sound speed and H the scale height. For sub-sonic turbulence with w ≪  c_s, our prescription gives ℓ ∼ w/N ≪ H. This is consistent with the expectation that stratification limits the vertical extent of turbulent structures to less than a scale height. Forced turbulence with t_o < 1/N is also possible. Since this would give even stronger dissipation, we consider _buoy a reasonable lower bound on dissipation for stratified turbulence.

Near the RCB, as N → 0 it is unreasonable to expect the dissipation to vanish entirely. Thus, we also include a floor to the dissipation

$$\begin{equation} \epsilon _o = f_\epsilon K_{zz} N_{\rm deep}^2 \end{equation} \tag{ 30 }$$

with the dimensionless normalization f giving the ratio of _o to _buoy in isothermal regions. If we use rotation as the other relevant timescale, then the floor would be quite low with f ∼ Ω²/N²_deep ∼ 10⁻³P^−7/3_day. Our full prescription considers both terms,

$$\begin{equation} \epsilon = \epsilon _o + \epsilon _{\rm buoy}, \end{equation} \tag{ 31 }$$

as discussed in Section 4.3.

We describe solutions to the atmospheric model of Section 3. We start by considering turbulent mixing with a constant diffusivity K_zz and no dissipation, i.e., = 0. This requires integration of Equation (23). While the decay of turbulence always gives some dissipation, the effect on energetics can be small (as we will show in Section 4.3).

Our main goal is to constrain the amount of turbulent diffusion that can be maintained in convectively stable regions. The simplest and most generous—i.e., allowing the largest levels of turbulent diffusion—constraint comes when we ignore dissipation.

An appeal to energetic balance shows why this should be the case. Recalling that turbulence drives a downward eddy flux, F_eddy < 0, in convectively stable regions, we rearrange flux balance as −F_eddy = F_rad − F. For a large downward eddy flux, we need the total flux F to be small compared to the outgoing radiative flux. Ignoring dissipation helps in this regard by preventing F from increasing through the layer. Pushing the RCB to high pressure also helps by lowering the constant F = F_c from the RCB. The solutions below show that strong mixing does indeed push the RCB to high pressure.

4.1.1. An Upper Limit to K_zz

Figure 3 shows the effect of varying K_zz while holding irradiation (T_deep = 1500 K) and internal entropy (T₁ = 250 K) fixed. As K_zz increases, the downwelling eddy flux pushes the RCB to higher pressures, P_c, and lowers the flux from the interior, F_c. These effects are coupled since T_c ∝ P^−6/7_c along a fixed adiabat (Equation (9)). Turbulent mixing—like strong stellar irradiation—reduces the planet's cooling rate.

Zoom In Zoom Out Reset image size

At a critical value of K_zz = K_zz,crit, P_c diverges to infinity while F_c drops to zero. This upper limit to diffusion is K_zz,crit ≈ 1660 cm² s^-1 for the adiabat and the T_deep chosen in Figure 3. Of course a real planet cannot extend to infinite pressure, to say nothing of the plane-parallel approximation. The point is that our steady state model cannot energetically support K_zz > K_zz,crit.

Stronger turbulence could in principle exist, since we are saying nothing about what forces turbulence. In a non-equilibrium state with K_zz > K_zz,crit, the downwelling eddy flux would then increase the internal entropy and inflate the planet.

Figure 4 shows that higher entropy planets have a higher steady state K_zz,crit. Thus by inflating the planet, strong mixing brings the planet's energy balance into equilibrium. An ultimate upper limit to K_zz is that the planet not overinflate and exceed its observed radius. The K_zz values invoked in S09, from 10⁷ to >10¹⁰ cm² s^-1 would imply significant or (on the upper end) excessive inflation. For comparison, AB06 showed (see their Figure 11) that entropy changes of ΔS ≈ k_B/m_p (the scale in our Figure 4) can expand a hot Jupiter's radius by ∼10%–25%. Accurate determination of the maximum K_zz allowed for a given planet requires more detailed modeling (including global structure with realistic opacities and EOS) than we perform. However, our results strongly suggest that K_zz values invoked in the literature have significant, or even excessive, effects on energetics.

Zoom In Zoom Out Reset image size

Figure 4 also shows that K_zz,crit increases with decreasing T_deep. Thus, our constraints on mixing are much more stringent for hot Jupiters than for more distant planets, including Jupiter itself. Recall that thermal inversions lower T_deep and that mixing can sustain thermal inversions by keeping opacity sources aloft in the stratosphere. A planet can accommodate strong mixing with some combination of thermal inversions to lower T_deep and increased internal entropy. It is hard to predict if thermally inverted planets should be more inflated—due to the presumed presence of turbulence—or less inflated—because lower T_deep promotes cooling and inhibits our mechanical greenhouse effect. Observations do not indicate an obvious correlation. Planets with signatures of inversions exhibit varying degrees of inflation. See Miller et al. (2009) for a comparison of observed to model radii of transiting planets.

While strong mixing pushes P_c → ∞, the depth of the isothermal layer, P_deep, is relatively unchanged by mixing. (We explore structure in detail below.) Thus, planets with higher entropy or lower T_deep have shallower isothermal layers. Specifically, $P_{\rm deep} \propto (T_{\rm deep}/T_1)^{1/\nabla _{\rm ad}}$ from Equations (13) and (16).

It is hardly surprising that planets which can accommodate more mixing (larger K_zz,crit) have shallower stratified layers to mix (smaller P_deep). In principle, strong mixing could destroy the deep isotherm altogether by pushing P_deep to optically thin regions. This probably requires unrealistically large core entropies. Alternatively, as interior temperatures rise, blackbody emission at depth may find opacity windows at wavelengths longer than the infrared.

4.1.2. Structure and Energetics of Stirred Atmospheres

We now consider the structure and energetic balance of "stirred atmospheres" with K_zz ≈ K_zz,crit. We compare these solutions to standard RE atmospheres of Section 2.2 with K_zz = 0. Since our model is self-similar, the behavior is independent of the parameters (T₁, T_deep) chosen for illustration.

Figure 5 (top panel) shows that the temperature profile of the stirred atmosphere is very similar to the RE case. In the stirred atmosphere, the RCB lies at much higher pressure below an extended "pseudo-adiabat,"⁷ which lies very close to the original adiabat. The inset to Figure 5 (top panel) focuses on the region near P_deep (which drops slightly to 480 bar from the original 610 bar) where the stirred atmosphere is hotter, by at most 60 K. The stirred atmosphere is slightly colder below 250 bar, though this difference of at most 3 K is not visible. The stirred atmosphere is very modestly thicker, by 0.06 H_deep, where H_deep = RT_deep/g ≈ 0.008 R_J.

Zoom In Zoom Out Reset image size

Figure 5 (bottom panel) plots the lapse rate. Mixing smoothes the transition toward the adiabat. The inset shows the smooth decline of 1 − ∇/∇_ad along the pseudo-adiabat, which gradually reduces the amplitude of the downwelling eddy flux, from Equation (20).

Figure 6 shows that the energetics of the stirred atmosphere differ significantly from the RE case. With K_zz = 0 there is no eddy flux and the radiative flux is constant down to the RCB, here at P_c ≈ 1.1 kbar. The stirred atmosphere has a deeper RCB which reduces the core flux significantly. We could push P_c → ∞ and F_c → 0 with a 0.1% increase of K_zz all the way to K_zz,crit, but choose not to for visualization.

Zoom In Zoom Out Reset image size

The radiative and eddy fluxes change with height. Their sum—the total flux—remains constant because we ignore dissipation. Figure 6 shows that the fluxes behave differently above and below P_deep, i.e., along the isothermal and pseudo-adiabatic regions, respectively. Along the pseudo-adiabat, the radiative flux declines with depth as F_rad ∝ P^−6/7, as we derived for the core flux in Equation (8). The radiative cooling (dF_rad/dP < 0) in this region balances heating by eddy diffusion (dF_eddy/dP > 0). Achieving this energetic balance requires only modest changes to the T–P profile. Since F_eddy scales as 1 − ∇/∇_ad, it is very sensitive to small changes in ∇ along the pseudo-adiabat (see Equation (20) and the bottom inset of Figure 5).

The energy balance along the isotherm, i.e., above P_deep is different. With ∇ ≪ 1, the eddy flux F_eddy ≈ −ρgK_zz ∝ −P grows in magnitude with depth (a different scaling holds if we vary K_zz with height, see Section 4.2). This localized cooling (dF_eddy/dP < 0) balances radiative heating (dF_rad/dP > 0). The decline in radiative flux with height again requires only modest changes to the thermal profile. From Equation (5), F_rad is sensitive to small changes in ∇ ≪ 1.

We can now simply estimate K_zz,crit using our knowledge that F_c → 0 and that mixing only modestly changes the RE profile. The transition region near P_deep is crucial. Here, the eddy flux reaches its peak negative value F_eddy,deep ≈ −ρ_deepgK_zz/2. The thermal profile constrains F_rad to be roughly F_c(K_zz = 0), the core flux of the RE atmosphere. We set F_rad,deep ≈ 2F_c(K_zz = 0) to account for the slightly hotter atmosphere near P_deep. Energetic balance, F_eddy,deep + F_rad,deep = F_c → 0, then gives

$$\begin{equation} K_{zz,{\rm crit}}\approx {4 F_{\rm c}(K_{zz}=0) \over \rho _{\rm deep} g}\propto \left({T_1^{1/ \nabla _{\rm ad}} \over T_{\rm deep}^{1/\nabla _{\rm ad}- 1/\nabla _\infty ^{\prime }}}\right)^{2+\alpha }, \end{equation} \tag{ 32 }$$

where ∇'_∞ = (2 + α)/(5 − β) and our parameters give K_zz,crit ∝ T^21/2₁/T^11/2_deep. These scalings agree with the results of Figure 4.

The above analysis (Section 4.1) shows that upper limits on a constant K_zz are set by the balance of eddy and radiative fluxes near P_deep, which itself scales with internal entropy and T_deep. To test the robustness of this finding, we include a depth dependence to K_zz ∝ P^ζ. For winds driven near the photosphere, i.e., the top of our atmospheres, one might expect stronger diffusion in the upper atmosphere, i.e., ζ < 0. On the other hand, if turbulence is triggered by shear layers with the convective interior, perhaps ζ>0. As discussed in Section 5.1, detailed dynamical simulations can help determine plausible diffusion profiles.

Figure 7 shows the effect of varying ζ with other parameters fixed (at our standard values of T_deep = 1500 K, T₁ = 250 K, α = 1, and β = 0 as in, e.g., Figure 3). These plots show the strongest possible mixing, which (as we found for constant K_zz) drives the RCB to infinite depths and reduces the core flux to zero.

Zoom In Zoom Out Reset image size

The maximum mixing near P_deep is relatively unchanged, except when the mixing at the top of the atmosphere is quite strong. Quantitatively we compare values of K_zz,deep, defined as the maximum value of K_zz at a reference P = 550 bar, which is P_deep of the RE atmosphere. For constant K_zz, we found K_zz,crit = K_zz,deep = 1665 cm² s^-1. For mixing that increases with depth as ζ = 0.5, 1.0, and 1.5, K_zz,deep declines by a modest 5%, 6%, and 5%, respectively. When mixing declines with depth as ζ = −0.5, −1.0, and −1.4, K_zz,deep increases by 14%, 58%, and 300%, respectively. We cannot consider models with ζ ≲ −1.5 because they do not approach an adiabat at depth. The bottom panel of Figure 7 shows that the approach to the adiabat is already quite gradual for ζ = −1.4. We explore this issue further in Appendix A.

The top panel of Figure 7 shows the flux profiles for several ζ values. The plot shows both radiative and eddy fluxes, which obey F_rad = −F_eddy because the net flux, F → 0, when mixing pushes the RCB to infinite depth. We also plot the radiative flux for the reference RE model (horizontal black dotted line) without any mixing. The explanation for these flux profiles mirrors the discussion of Figure 6 in Section 4.1. At high pressures, the flux is controlled by the radiative flux along the adiabat. Increasing or decreasing the mixing with depth has little effect on the deep eddy flux. Changes to K_zz are compensated by (1 − ∇/∇_ad)—see Equation (20)—which is small and sensitive to slight changes in ∇ close to the adiabat.

The flux in the low pressure, isothermal region scales as −F_eddy ∝ ρK_zz ∝ P^1+ζ (see Equation (20)). This explains why the flux components, F_rad = −F_eddy, increase with height if ζ < −1. Driving larger radiative fluxes in the upper atmosphere requires a steeper dT/dP. The temperature profiles in Figure 7 (bottom panel) reflect this. The ζ = −1 and especially the top ζ = −1.4 curves are noticeably hotter at intermediate pressures and have smaller P_deep. We thus find that mixing at the top of the atmosphere is more effective—compared to uniform or bottom-focused mixing—at lifting (i.e., heating) the T–P profile. The additional heat in this case is provided by a downward flux of mechanical energy across the top boundary.

4.2.1. Limits on Mixing Near the Photosphere

We now consider what might constrain K_zz near the top of the atmosphere, since we find that the internal entropy mostly constrains diffusion near P_deep. Our ζ ⩽ −1 solutions in Figure 7 show that strong mixing at the top requires a large flux of mechanical energy at the top of the atmosphere. Weather-layer winds are a plausible source of mechanical energy, and they are generated by the atmospheric heat engine driven by insolation. The thermodynamic efficiency of all planetary atmospheres in the solar system is of order 1%. If we restrict the magnitude of F_eddy to a fraction f_* ∼ 1% of the insolation F_* ∼ σT⁴_deep, we get

$$\begin{eqnarray} K_{zz,{\rm top}} &<& {f_\ast F_\ast \over \rho _{\rm top} g} \nonumber\\ &\approx & 10^9 {\rm cm^2 \over s}\left(P_{\rm top} \over 0.1 {\rm bar}\right)^{-1} \left(T_{\rm deep} \over 1500\, {\rm K}\right)^5{f_\ast \over 1\%}, \end{eqnarray} \tag{ 33 }$$

scaled for a downwelling flux that originates at a 0.1 bar photosphere. While the efficiency f_* and mechanisms of generating a downward mechanical flux are uncertain, the energetic difficulties of mixing at K_zz ≫ 10⁹ cm² s^-1 are evident.

Alternatively, we could attempt to constrain mixing in the upper atmosphere by appealing to our ζ = −1.4 solution. This gives the largest downward eddy flux (Figure 7, top). Smaller ζ values would be needed for a larger flux, but these do not give consistent solutions. In Appendix A, we analyze the ζ ≈ −1.4 limit and argue that it may not be physically significant. Ignoring this concern, the diffusion near the top of our ζ = −1.4 solutions is K_zz ≈ 5 × 10⁷(P/bar)^−1.4 cm² s^-1. A larger internal entropy could support a larger K_zz (as we showed for constant K_zz models in Figure 4). Therefore, this constraint is not inconsistent with Equation (33), which is a more robust constraint.

We now consider the effect of adding dissipation to our models with eddy diffusion. The total flux F will now increase with height due to dissipation. The coupled Equations (23) and (24) govern the steady state structure. To understand the effect of dissipation on the turbulent heat flux, we return to the simpler case of spatially uniform K_zz.

With dissipation we still find an upper limit to diffusion, K_zz,crit, for a given T₁ and T_deep. However, K_zz,crit declines with increasing dissipation. Figure 8 shows this for both a dissipation rate _o that is constant with height (solid blue curve) and the full dissipation prescription (dashed red curve; see Equation (31)). The full prescription includes _buoy, our estimate of the minimum dissipation due to stratified turbulence. This additional dissipation further reduces K_zz,crit.

Zoom In Zoom Out Reset image size

We confirm that the dissipation-free estimates of K_zz,crit in previous sections represent a conservative upper bound. Lowering K_zz,crit means that weaker turbulent diffusion will inflate the planet. Admittedly, the cases shown in Figure 8 do not prove that all dissipation profiles will lower K_zz,crit. However, we investigated the effects of both spatially varying K_zz (as in Section 4.2) and also different profiles of . In all cases, adding dissipation reduced K_zz,crit from the dissipation-free value.

4.3.1. How Dissipation Lowers K_zz,crit

We explore the energetics of how dissipation lowers K_zz,crit. Figure 9 shows the flux balance with dissipation and can be compared to Figure 6. Note that the peak value of −F_eddy now falls well short of F_rad at the relevant pressure, P_deep. This is because the total flux

$$\begin{equation} F = F_{\rm c} + \int ^{P_{\rm c}}_P {\epsilon \over g} dP \end{equation} \tag{ 34 }$$

now includes the integrated dissipation, causing F to greatly exceed F_c, the small loss of heat from the core. From the RCB to P_deep, F_rad also increases with height. This rise is not significantly affected by dissipation, with F_rad ∝ P^−6/7_c along the pseudo-adiabat as before. The magnitude of −F_eddy = F_rad − F is smaller at P_deep because dissipation increases F without comparably increasing F_rad.

Zoom In Zoom Out Reset image size

We find that P_deep (and also ρ_deep) does not significantly change when we add dissipation. Indeed, the T–P or ∇ profiles for the solution in Figure 6 are indistinguishable from the stirred atmospheres in Figure 5, except the RCB is not pushed as deep, "only" to P_c ≈ 11 kbar.

The limiting value of K_zz = 2|F_eddy|/(ρ_deepg) drops because the smaller eddy flux is not compensated by a lower ρ_deep. It seems possible that some dissipation profile could heat the atmosphere and lower P_deep and ρ_deep enough to increase K_zz,crit. We did not find this to be the case. One reason is that too much dissipation can affect the location of the RCB, a subject we address below.

For completeness, we explain the flux balance at small pressure illustrated in Figure 9. The decline in −F_eddy ≈ ρgK_zz with height in isothermal regions again results from declining density. However, F_rad does not decline toward low pressure (as shown in Figure 6), because the total flux F is larger with dissipation. The fact that the escaping flux matches the flux from the RE solution (dotted black line) is a coincidence (with some significance, see below). This coincidence occurs because the dominant dissipation is ≈ _buoy. Larger or smaller choices of dissipation would give a larger or smaller (respectively) net F and escaping F_rad.

This coincidence has some significance. Downward eddy fluxes reduce the loss of heat from the core, as we have discussed extensively. However, we now see that this loss is matched by the flux due to the dissipation of that turbulence—provided our prescription for the minimum _buoy is correct. This replacement is intriguing, but does not alter our discussions of evolutionary consequences: turbulent dissipation in radiative regions is powered not by interior heat, but by external means (such as forced atmospheric circulation).

4.3.2. The Effect of Dissipation on the RCB

All of our plots and numerical estimates have used an opacity law κ ∝ P and an ideal gas EOS with ∇_ad = 2/7. We provide general scalings for many results to show the effect of varying these parameters.

Our qualitative results hold for all "reasonable" choices of power-law opacities and EOS. As discussed in Section 2.1, reasonable means that ∇_∞ > ∇_ad and α> − 1 so that RE solutions become convectively unstable at depth. We tested our results with an alternate opacity law κ ∝ T² (α = 0 and β = 2) that approximates H⁻ opacities for T > 2000 K or dust grain opacities at colder temperatures. Like κ ∝ P, this law also has ∇_∞ = (1 + α)/(4 − β) = 1/2>∇_ad.

The behavior of K_zz,crit—shown in Equation (32)—is particularly important for interpreting our results. The scaling with internal entropy is quite steep with K_zz,crit ∝ T^10.5₁ and K_zz,crit ∝ T⁷₁ for our standard and alternate opacities, respectively. Generally, the entropy dependence becomes steeper for larger α (as in the above example) and also for a smaller ∇_ad, but does not depend on β. Recall that the burial of the turbulent heat flux into the convective interior can bring an atmosphere with K_zz > K_zz,crit toward energetic equilibrium. A steeper dependence of K_zz,crit on entropy means that less inflation is needed to enforce this equilibrium.

The scaling of K_zz,crit with T_deep controls how our inflation mechanism depends on the level of irradiation. Also, the development of thermal inversions will lower T_deep for a fixed level of irradiation. We find K_zz,crit ∝ T^−5.5_deep or K_zz,crit ∝ T⁻⁴_deep again for the standard and alternate opacities, respectively. The T_deep dependence becomes more steeply negative for larger α (again the dominant effect in our example) and also for larger β (less important in our example) and smaller ∇_ad.

These scalings emphasize that a more detailed treatment of turbulent heat fluxes in hot Jupiters should include realistic opacities and EOSs. Non-power-law behavior could have significant consequences. For instance, Guillot et al. (1994) demonstrated the importance of an opacity window near ∼2000 K. This window can give rise to an isolated convective layer sandwiched between two radiative zones. It would be interesting to consider how a downward turbulent heat flux would interact with such a region. Future work that generalizes our self-similar approach can address these more detailed issues.

Hydrodynamic simulations of hot Jupiter atmospheres have been studied in local (Burkert et al. 2005; Li & Goodman 2010, hereafter LG10) and global (e.g., Cooper & Showman 2005; Rauscher & Menou 2010; Showman et al. 2009b; Dobbs-Dixon & Lin 2008) models. A major goal of these studies is to determine the circulation induced by stellar irradiation. We first discuss estimates of K_zz from these simulations and then address the question of whether radiatively forced turbulence extends throughout the radiative zone, as we have assumed. Other sources of turbulence—perhaps involving magnetic fields—may also exist, but we do not analyze them here.

The global circulation models (GCMs) of Showman et al. (2009b) estimate K_zz ∼ 10¹¹ cm² s^-1 at mbar pressures. This does not contradict our constraint on upper atmospheric mixing from the efficiency of radiative forcing (Equation (33)) when extrapolated to such low pressures. Moreover, at low pressures radiative losses can lower the eddy flux (see Section 3.1) and further weaken our constraint on K_zz. We caution that K_zz estimates from GCMs are rough, since they only resolve relatively large scale flow patterns. The Showman et al. (2009b) K_zz estimate arises from multiplying measured vertical speeds by the scale height H. While a useful guide to what is possible, this does not constitute a direct measurement of diffusion.

Local hydrodynamic simulations can better resolve turbulent flows, though as usual with Reynolds numbers far lower than reality. The calculations of LG10 find an effective turbulent viscosity of ν_t ∼ 0.001–0.01c_sH ∼ 10¹⁰–10¹¹ cm² s^-1. It is unclear if this viscosity should be interpreted as a mixing coefficient. Their simulation with ν_t = 0.015c_sH has an rms vertical speed w ∼ 0.3 c_s. Thus, assuming ν_t ∼ K_zz is not consistent with the simple estimate K_zz ∼ wℓ because it gives a length scale ℓ ∼ H/20, smaller than the grid spacing of H/10. It should not be surprising that simple estimates based on three-dimensional isotropic turbulence fail, given the organized structure in their two-dimensional "turbulent" state.

Moreover, the forcing in LG10 was not by irradiation, but chosen to be large enough to drive super-sonic flows despite artificial viscosities that are large for numerical reasons. Thus, attempting to interpret the Carnot efficiency of stellar irradiation may overextend their results. Despite these caveats, if the LG10 simulations apply near the millibar weather layer, there is again no contradiction with Equation (33).

Global simulations indicate that the shear layer may not extend throughout the radiative zone. For instance, Showman et al. (2009b) find that strong (∼ km s⁻¹) zonal winds terminate at ∼10 bar. They argue that circulation stops because the planet is horizontally isothermal at these depths, removing the local forcing. However, hot Jupiter atmospheres have a large separation between the radiative timescales in the weather layer and the deep radiative layer. It is possible (and arguably likely) that unresolved or long timescale dynamics could push the shear layer even deeper. This possibility should be explored in future modeling studies.

How far turbulence can extend below the shear layer is also uncertain. LG10 force a shear layer of with ∼2H, but find that turbulence (or at least some disordered motion) extends throughout their box of size ∼5H (see their Figure 10). Turbulence that extends a full 5H below a shear layer could thus extend to P ∼ e⁵ · 10 bar ∼ 1.5 kbar. These issues deserve further study, but it cannot be ruled out that the radiative zone is turbulent down to the RCB.

We now provide a more detailed interpretation of our results in terms of the S09 constraints on the mixing required to loft TiO and create thermal inversions. We emphasize that the constraints on K_zz and especially on the depth dependence of K_zz depend sensitively on whether or not there is a cold trap. Wherever TiO condenses, one must consider the mixing of dust grains, not just molecules.

First, consider the case where there is no cold trap, and TiO is always in the vapor phase. In the S09 analysis, this is only possible for the most intensely irradiated planet, WASP 12b, in part because thermal inversions lower T_deep and favor condensation at depth. The mixing of TiO vapor requires K_zz ∼ 10⁷ cm² s^-1 at P ∼ 1 mbar, i.e., the height of the inversion where TiO is needed. The specification of pressure level is important because the constraint on molecular diffusion scales as K_zz ∝ P⁻¹ in an isothermal atmosphere. Thus, the constraint is weaker at depth. In a model where the actual K_zz ∝ P⁻¹ (as in Section 4.2), we only require K_zz ∼ 10 cm² s^-1 at the kbar RCB. In this case, our model predicts a minimal effect of the eddy flux on the planet's evolution, though the dissipation of turbulence above the RCB could still be significant.

We briefly summarize how this K_zz ∼ 10⁷ cm² s^-1 limit and the depth dependence arise, and refer the reader to S09 for details. With only molecular viscosity, i.e., collisions, the TiO scale height would be hydrostatic, and thus ∼30 times smaller than the dominant H₂ species. S09 conclude that the turbulent diffusion must exceed the molecular diffusion D_TiO by a factor ∼100, due to the many (∼14) scale heights between P∼ mbar and the kbar RCB. Estimating D_TiO as usual—the product of the mean free path and thermal speed—gives an inverse scaling with density, and thus also with pressure in isothermal regions. The numerical value of D_TiO ∼ 10⁵ cm² s^-1 at P∼ mbar, combined with the factor of 100 excess needed for efficient turbulent mixing, gives the K_zz ∼ 10⁷ cm² s^-1 constraint.

Now consider the case where grains do condense somewhere in the radiative zone, which is true for most hot Jupiters (especially those with inversions). The K_zz required to loft TiO increases to 10⁸–10¹¹ cm² s^-1 depending on the size of grains and the depth of the cold trap. A rough estimate of K_zz ∼ v_termL follows from equating the diffusion timescale, L²/K_zz, to the dust settling timescale L/v_term, where L is the depth of the cold trap and v_terms the grain's terminal speed. While the terminal speed will increase with particle size, it does not depend on atmospheric density for the relevant viscous (i.e., Stokes') drag. Thus, unlike the case of molecular viscosity, the K_zz required for mixing condensed grains does not decline with depth.

We can thus conclude that the K_zz required to mix TiO and create thermal inversions in a hot Jupiter are likely excessive, provided TiO condenses somewhere in the radiative zone. This follows from Figure 4 which shows that K_zz ∼ 10⁹ cm² s^-1 appear off-scale. The internal (specific) entropy of a planet would have to increase by ≳2k_B/m_p for a planet with T_deep ≳ 1500 K, as is the case for all planets in S09. It seems likely that such large entropies would overinflate a planet, based on studies such as in AB06. This is especially true because Figure 4 neglects dissipation which can further inflate a planet, which can occur by lowering K_zz,crit as shown in Figure 8.

However, we cannot firmly declare that the TiO hypothesis fails. This is largely due to the approximate nature of our treatment of opacities and the EOS. A more conclusive analysis of the TiO hypothesis would require a detailed atmospheric model that includes the eddy fluxes and turbulent dissipation described in this work. Such a study would also have to abandon the fixed flux bottom boundary condition used in most atmospheric models (including S09). Instead, the bottom boundary should be a fixed adiabat, chosen to match the observed radius. This is subject as usual to assumptions about composition and presence of a core. However, allowing the flux to float is needed for a consistent determination of the RCB. This is crucial for understanding the coupled relationship between compositional mixing and cooling history.

We have investigated how forced turbulent mixing affects the energetic balance and structure of the stably stratified, radiative layers of hot Jupiters. It is crucial to understand how the radiative layer matches onto the convective interior at the RCB. This regulates the rate at which the planet cools and therefore controls the evolution of the planet's radius.

Previous work has invoked turbulent eddy diffusion of molecular species (and dust). This mixing changes the opacity to provide better-fitting model spectra of transiting planets, especially those that appear to have thermal inversions (S09). We find that turbulent mixing of this kind does not just redistribute chemical species but also significantly affects energetics.

Forced turbulent mixing in stable, radiative regions drives a downward flux of energy that pushes the RCB deeper in the atmosphere, lowering the planet's cooling rate. We found an upper limit to the strength of turbulent diffusion, K_zz,crit, that can be achieved in steady state. Beyond this limit, the downward flux of energy will heat the convective interior and inflate the planet. We did not directly model this entropy growth because our model was steady state and did not include overshoot across the RCB. The deep, marginally stable layer in our solutions would not strongly inhibit overshoot, so heat burial by this mechanism is a likely outcome. Our solutions indicate that interior heating brings the planet back toward steady state, because higher entropy planets have larger K_zz,crit (see Figure 4). Our mechanism is thus a "mechanical greenhouse effect" with the role of sunlight in the traditional greenhouse being played by forced turbulent mixing.

For non-uniform turbulent mixing, we find that our constraint on K_zz,crit applies near the RCB. More specifically, this constraint applies at a pressure P_deep, which lies below the deep isothermal region where the radiative layer is transitioning toward convective instability. Our constraints on turbulence in the upper atmosphere are less stringent. However, the downward flux of mechanical energy likely cannot exceed a small fraction (probably at the percent level) of the stellar irradiation if it is supplied by weather-layer winds. If turbulence is too weak at the bottom of the radiative layer, it will not dredge up heavy molecular species—particularly those that condense onto dust grains—to serve as opacity sources near the observable photosphere.

Turbulence also deposits heat in the radiative atmosphere when it decays. Non-turbulent sources of energy dissipation—including nonlinear wave breaking and ohmic dissipation—also affect energetic balance. We find that including energy dissipation reduces K_zz,crit, and thus makes it easier to inflate a planet for a given level of forced turbulence.

We find a characteristic scale of K_zz,crit ∼ 10³–10⁵ cm² s^-1 for typical hot Jupiter parameters, even ignoring dissipation. This is orders of magnitude below values quoted in the literature of 10⁷–10¹¹ cm² s^-1 for the mixing of chemical species (Spiegel et al. 2009; Zahnle et al. 2009). We caution that our quantitative results should be taken as illustrative, due to the approximations detailed in Section 2.

Thus, when turbulence is important for the redistribution of opacity sources, it also has significant energetic consequences. Our model represents a step toward more self-consistent modeling of exoplanet atmospheres.

Guillot & Showman (2002) first noted that dissipation in radiative layers can slow planetary contraction by pushing the RCB to higher pressures. They imagined dissipation concentrated in the upper regions of the atmosphere, sourced by insolation driven winds. Other authors (e.g., Bodenheimer et al. 2003) have included the Guillot & Showman (2002) dissipation prescription to model exoplanet radii.

However, we are unaware of previous works that include the mechanical flux of energy due to turbulence or waves—our F_eddy—in radiative layers. Incorporating this effect in detailed planetary evolution models would require a more precise treatment than this exploratory study. Notably it would require a global model with a realistic EOS and opacity, not the power laws considered here.

In this paper, we treat diffusion and dissipation by forced turbulence as free parameters to allow a general analysis. This approach is justified since we currently lack a detailed understanding of these processes. Thus, our model can be used to test how any specific model for turbulence and/or energy dissipation affects the structure and evolution of the planet. For instance, the work of Batygin & Stevenson (2010) considers ohmic dissipation in radiative and convective regions, but does not recompute the effect of the dissipation on the structure of the radiative layer. Our point is not to critique, but to emphasize that a more consistent treatment of energetics may make it easier to inflate a planet—by any number of mechanisms.

We also hope to include the effects of our study in detailed three-dimensional radiation-hydrodynamical simulations of exoplanet atmospheres. This would involve adding sub-grid physical prescriptions for eddy fluxes to GCMs such as those by Showman et al. (2009b) and Rauscher & Menou (2010) and/or non-hydrostatic simulations (Dobbs-Dixon et al. 2010). Goodman (2009) discusses the need to include explicit sources of dissipation. Li & Goodman (2010) present a first step toward sub-grid modeling of turbulence due to Kelvin–Helmholtz instabilities. Breaking of vertically propagating gravity waves represents another possible source of turbulence (Showman et al. 2009a). When all these effects are taken into account, the inflated radii of transiting planets may not be surprising after all.

We have focused on hot Jupiters, but the physics we describe in principle applies to other atmospheres. Our limits on eddy diffusion become less severe as objects are more weakly irradiated (see Figure 4). Thus, our mechanism would not inflate longer period gas giants, including Jupiter and Saturn. More intrinsically luminous objects like brown dwarfs will similarly be less affected. The expanding inventory of transiting exoplanets, at wider separations, is an excellent test of radius evolution models.

Finally, we comment on a possible connection to the atmosphere of Venus. Venus has a marginally stable pseudo-adiabat beneath a thick cloud deck (Schubert et al. 1980). The clouds shield sunlight from warming the surface of Venus. Because of this, Venus' atmosphere would be nearly isothermal were it not for some poorly understood mechanical stirring process. Our well-stirred hot Jupiter atmospheres also have marginally stable pseudo-adiabats at depth. Perhaps Venus is a very well-stirred analog of a hot Jupiter. Obvious differences exist, for instance, the non-negligible fraction of sunlight that reaches the Venusian surface versus the radiant flux escaping the core of hot Jupiters. Pursuing analogies such as these should improve our understanding of worlds near and far.