A Scaling for Atmospheric Heat Redistribution on Tidally Locked Rocky Planets

Telescope observations have recently begun to probe the atmospheres of small, rocky, exoplanets around nearby stars (Demory et al. 2016; de Wit et al. 2018; Diamond-Lowe et al. 2018; Kreidberg et al. 2019). These observations will only get better in the near future with the launch of the James Webb Space Telescope (JWST) and the construction of extremely large ground-based telescopes. The data that we will be able to obtain with these instruments promise insight into fundamental questions such as how rocky exoplanets form, how many of them host atmospheres, and whether any of them might be habitable.

To interpret these observations correctly, however, we need to understand the chemical and physical processes that shape the observable features of rocky exoplanets. These processes include radiative transfer and molecular absorption (Seager & Sasselov 2000), gas-phase chemistry (Moses et al. 2011; Hu & Seager 2014), clouds and hazes (Hörst et al. 2018; Moran et al. 2018), as well as large-scale atmospheric dynamics (Joshi et al. 1997; Merlis & Schneider 2010). Of these processes, atmospheric dynamics has an outsized impact on thermal observations, because it determines a planet’s global heat redistribution—which sets the depth of the planet’s secondary eclipse as well as the shape of the planet’s thermal phase curve.

Unfortunately, heat redistribution remains poorly represented in most models that are being used to match and interpret exoplanet observations. The underlying reason is that these retrieval models have to be fast enough to be run ∼10⁵–10⁶ times (Madhusudhan & Seager 2009; Line et al. 2013), which is required to comprehensively map out which atmospheric scenarios can, or cannot, fit an observed data set. As a consequence, most retrieval models are idealized one-dimensional representations of a planet’s atmosphere (e.g., Tinetti et al. 2007; Madhusudhan & Seager 2009; Benneke & Seager 2012; Line et al. 2013; Kempton et al. 2017). By virtue of being 1D, these models cannot resolve 3D processes such as clouds and atmospheric heat redistribution, which has been shown to bias the results of atmospheric retrievals as well as model mean states (Feng et al. 2016; Line & Parmentier 2016; Fauchez et al. 2018).

In principle, 3D processes including atmospheric heat redistribution can be resolved by more complex models, such as general circulation models (GCMs). In practice, 3D models are computationally far too costly to be used in place of 1D retrieval models. For example, the GCM used later in this paper requires about 24 hr of computation time on eight processors. Even a dedicated supercomputer, running 32 GCM simulations in parallel for an entire month, would thus only be able to explore ∼10³ parameter combinations, falling far short of being useful for standard retrieval techniques.

An important question is therefore how exoplanet retrieval models should represent the basic physics that determines observable quantities such as a planet’s dayside emission or its thermal phase curve, while remaining computationally cheap. Most 1D models represent the effects of heat redistribution by adjusting the planet’s dayside energy budget, which can be written as (Burrows 2014)

$$\begin{eqnarray}\begin{array}{rcl}{T}_{\mathrm{day}} & = & {T}_{* }\sqrt{\displaystyle \frac{{R}_{* }}{d}}{\left(1-{\alpha }_{B}\right)}^{1/4}{f}^{1/4},\\ & = & {\left(4f\right)}^{1/4}{T}_{\mathrm{eq}}.\end{array}\end{eqnarray} \tag{ 1 }$$

Here, T_day is the observed dayside brightness temperature, ¹ T_* is the stellar temperature, R_* is the stellar radius, d is the planet’s semimajor axis, α_B is the planet’s bond albedo, and f is the so-called heat redistribution factor. It is useful to re-express T_day in terms of the equilibrium temperature ${T}_{\mathrm{eq}}={T}_{* }\sqrt{{R}_{* }/\left(2d\right)}{\left(1-{\alpha }_{B}\right)}^{1/4}$, which is the temperature of an isothermal sphere with the same bond albedo and semimajor axis as the planet. The temperature T_eq is closely related to the planet’s characteristic fluxes, in particular $\sigma {T}_{\mathrm{eq}}^{4}$, the flux received from the host star averaged across the entire planet, $2\sigma {T}_{\mathrm{eq}}^{4}$, the stellar flux averaged across the planet’s dayside only, and $4\sigma {T}_{\mathrm{eq}}^{4}$, the incoming flux at the planet’s substellar point. The redistribution factor f has to lie between 2/3 for a planet without an atmosphere and 1/4 for a planet that is extremely efficient at redistributing heat (Hansen 2008), but currently there no easy way of expressing f in terms of a planet’s physical properties.

The goal of this paper is to develop a simple scaling theory for f. The derivation is shown in Section 2. The results demonstrate that f is sensitive to a planet’s equilibrium temperature, surface pressure, and its broadband optical thickness, although the scaling also describes the dependency on other physical parameters. In Section 3, the scaling is compared against GCM simulations of three nearby rocky planets: TRAPPIST-1b, GJ1132b, and LHS3844b (Berta-Thompson et al. 2015; Gillon et al. 2016; Delrez et al. 2018; Vanderspek et al. 2019), which are among the highest-priority targets for upcoming thermal observations with JWST. The scaling successfully captures the dominant processes that determine atmospheric heat redistribution on these planets. The scaling therefore lends insight into the atmospheric dynamics of tidally locked planets. Moreover, it can be used as a computationally efficient parameterization of large-scale dynamics in 1D models. Because the derivation and GCM simulations assume idealized semi-gray radiative transfer, Section 4 considers how this work can be applied to real gases with non-gray absorption. The results are discussed in Section 5, and the conclusions are in Section 6.

This section derives a day–night heat redistribution scaling by combining two arguments. First, tidally locked planets are relatively slow rotators, so their atmospheric dynamics are to first order in weak-temperature-gradient (WTG) balance (Pierrehumbert 2011; Mills & Abbot 2013). Second, planetary atmospheres resemble heat engines, which means their day–night circulation strength and wind speeds can be approximately predicted using Carnot’s theorem (Koll & Abbot 2016; Koll & Komacek 2018). The combination of these two arguments allows us to analytically predict a planet’s day–night heat redistribution, at least in the limit of thin atmospheres with inefficient heat redistribution. The prediction is then extended to arbitrary atmospheric thicknesses using an Ansatz. The choice of an Ansatz solution can only be justified retrospectively. It does appear useful, however, as the scaling’s predictions compare favorably against a wide range of GCM simulations.

The core idea of WTG balance is to reduce the atmosphere’s full equation of energy conservation into a balance between radiative cooling and warming by adiabatic compression. We can expect this balance to hold for parcels of air that are undergoing large-scale descent, in particular the atmosphere on the nightside of a slowly rotating tidally locked planet for which horizontal temperature gradients are small. Using the gray optical depth τ as the vertical coordinate, the balance can be written as (see Koll & Abbot 2016)

$$\begin{eqnarray}&&\displaystyle \frac{{c}_{p}\omega }{g}\left(\displaystyle \frac{{dT}}{d\tau }-\beta \displaystyle \frac{T}{\tau }\right)=\displaystyle \frac{{dF}}{d\tau }.\end{eqnarray} \tag{ 2 }$$

Here, T(τ) is the atmosphere’s nightside vertical temperature profile, F(τ) is the net infrared flux, ω is the vertical velocity, c_p is the atmosphere’s specific heat capacity, g is the acceleration of gravity, β ≡ R/(c_p n_LW) is the dry adiabatic lapse rate in optical depth coordinates, R is the atmosphere’s specific heat constant, and n_LW = 1 if opacity is independent of pressure (e.g., if molecular line widths are set by thermal broadening) while n_LW = 2 for pressure broadening. The vertical velocity is defined in terms of a parcel’s rate of pressure change, ω ≡ Dp/Dt, and for sinking air ω > 0. The optical depth τ varies from zero at the top of the atmosphere to τ_LW at the surface. If one evaluates F at the top of the atmosphere, one obtains the nightside’s flux of outgoing longwave radiation, F_n ≡ F(τ = 0). As long as internal heat fluxes (e.g., tidal heating) are negligible, this has to equal the day–night heat transport by the atmosphere. Therefore, if we can use Equation (2) to constrain F_n , we will also know how much heat the atmosphere is redistributing.

It is difficult to solve Equation (2) in general, because the infrared flux F(τ) and temperature profile T(τ) are closely coupled. However, previous work showed from simple energy-balance considerations that, in the limit of a thin atmosphere, the atmosphere becomes roughly isothermal and approaches the gray skin temperature, so T(τ) ≈ T_skin = 2^−1/4 T_eq (Wordsworth 2015). The underlying reason is that an optically thin atmosphere still absorbs a trickle of radiation from the warm dayside surface—which, as long as horizontal temperature gradients are small, allows the entire atmosphere to equilibrate at the temperature T_skin. A roughly isothermal atmosphere also significantly simplifies Equation (2),

$$\begin{eqnarray}&&-\displaystyle \frac{{c}_{p}\omega }{g}\beta \displaystyle \frac{{T}_{\mathrm{skin}}}{\tau }\approx \displaystyle \frac{{dF}}{d\tau }.\end{eqnarray} \tag{ 3 }$$

If we specify ω, we can integrate Equation (3) to find the top-of-atmosphere flux F_n . In general, one can expect ω = f(τ) × ω₀(p_s , …), where f(τ) needs to go to zero at the top of the atmosphere and at the surface to ensure that the atmospheric flow vanishes at those boundaries, while ω₀ captures the still-to-be-determined dependence of the circulation on all other atmospheric parameters (such as surface pressure p_s ). The detailed structure of f(τ) does not matter much, because it only ends up yielding a numerical constant. For example, assuming $f(\tau )=\sin (\pi \tau /{\tau }_{\mathrm{LW}})$ and using the boundary condition that the nightside surface has to be in radiative equilibrium, F(τ_LW) = 0, integration leads to:

$$\begin{eqnarray}\begin{array}{rcl}F({\tau }_{\mathrm{LW}})-{F}_{n} & \approx & -\displaystyle \frac{{c}_{p}{\omega }_{0}}{g}\beta {T}_{\mathrm{skin}}{\displaystyle \int }_{0}^{{\tau }_{\mathrm{LW}}}\displaystyle \frac{\sin (\pi \tau /{\tau }_{\mathrm{LW}})}{\tau }d\tau ,\\ {F}_{n} & \approx & \displaystyle \frac{{c}_{p}{\omega }_{0}}{g}\beta {T}_{\mathrm{skin}}\times 1.85,\\ & \approx & \displaystyle \frac{{c}_{p}{\omega }_{0}}{g}\beta {T}_{\mathrm{skin}},\end{array}\end{eqnarray} \tag{ 4 }$$

where we have dropped the numerical factor of order unity.

Next, we introduce Carnot’s theorem to constrain the circulation strength ω₀. In steady state, an atmospheric heat engine needs to balance the generation of work against frictional dissipation, so Carnot’s theorem can be written as

$$\begin{eqnarray}&&{\rho }_{s}{C}_{d}{U}^{3}=\eta Q.\end{eqnarray} \tag{ 5 }$$

Here, the left-hand side expresses the rate of frictional dissipation, ρ_s = p_s /(RT_s ) is the air density near the surface, C_d is a drag coefficient, U is the near-surface wind speed, η = (T_hot − T_cold)/T_hot is the heat engine’s efficiency, and Q is the heat flux absorbed by the atmosphere. In the limit of a thin atmosphere, the dayside surface is heated by stellar radiation, but can re-emit most of this energy directly back to space as longwave radiation. The remaining flux that is available to heat the atmosphere is therefore only $Q=(1-{e}^{-{\tau }_{\mathrm{LW}}})\times 2\sigma {T}_{\mathrm{eq}}^{4}\approx {\tau }_{\mathrm{LW}}\times 2\sigma {T}_{\mathrm{eq}}^{4}$, where $(1-{e}^{-{\tau }_{\mathrm{LW}}})$ is the fraction of radiation not directly re-emitted to space according to Beer’s law and $2\sigma {T}_{\mathrm{eq}}^{4}$ is the incoming stellar flux averaged across the planet’s dayside. If the areas of upwelling air on the dayside and downwelling air on the nightside are roughly equal, the vertical velocity on the nightside is then simply determined by mass balance, U/L ≈ ω₀/p_s , where L is the circulation’s horizontal scale. Combining the expressions in this paragraph with Carnot’s theorem (Equation (5)) yields the following vertical velocity scale:

$$\begin{eqnarray}&&{\omega }_{0}=\displaystyle \frac{{p}_{s}}{L}{\left(\displaystyle \frac{\eta {\tau }_{\mathrm{LW}}2\sigma {T}_{\mathrm{eq}}^{4}}{{\rho }_{s}{C}_{d}}\right)}^{1/3}.\end{eqnarray} \tag{ 6 }$$

Plugging back into Equation (4), the day–night heat redistribution of a thin atmosphere is thus

$$\begin{eqnarray}\begin{array}{rcl}{F}_{n} & \approx & \chi \displaystyle \frac{{c}_{p}\beta }{{gL}}{\left(\displaystyle \frac{{{Rp}}_{s}^{2}{\tau }_{\mathrm{LW}}\sigma {T}_{\mathrm{eq}}^{8}}{{C}_{d}}\right)}^{1/3},\\ {F}_{n} & \approx & \chi \displaystyle \frac{{c}_{p}\beta }{{gL}}{\left(\displaystyle \frac{{{Rp}}_{s}^{2}{\tau }_{\mathrm{LW}}}{{C}_{d}{\sigma }^{2}{T}_{\mathrm{eq}}^{4}}\right)}^{1/3}\sigma {T}_{\mathrm{eq}}^{4},\end{array}\end{eqnarray} \tag{ 7 }$$

where we have absorbed all numerical constants as well as the heat engine efficiency η into the factor χ. Because the values that enter χ are all roughly order unity, and because most of these values are further raised to the one-third power, one might expect that χ should also be of order one. In practice, however, we do not try to predict the value of χ from first principles. Instead, we treat χ as the major source of uncertainty when comparing the scaling to numerical simulations (see below).

Finally, we seek to generalize the asymptotic solution for thin atmospheres in Equation (7) using an Ansatz. To do so, note that Equation (7) cannot be generally valid because it predicts that F_n → ∞ as p_s or τ_LW become large, whereas energy conservation requires that ${F}_{n}\to \sigma {T}_{\mathrm{eq}}^{4}$ for a thick atmosphere. A general solution needs to reduce to the asymptotic solution as p_s and τ_LW become small, but tend toward $\sigma {T}_{\mathrm{eq}}^{4}$ as one or both parameters become large. The following expression satisfies both requirements:

$$\begin{eqnarray}&&{F}_{n}=\displaystyle \frac{\chi \tfrac{{c}_{p}\beta }{{gL}}{\left(\tfrac{{{Rp}}_{s}^{2}{\tau }_{\mathrm{LW}}}{{C}_{d}{\sigma }^{2}{T}_{\mathrm{eq}}^{4}}\right)}^{1/3}}{1+\chi \tfrac{{c}_{p}\beta }{{gL}}{\left(\tfrac{{{Rp}}_{s}^{2}{\tau }_{\mathrm{LW}}}{{C}_{d}{\sigma }^{2}{T}_{\mathrm{eq}}^{4}}\right)}^{1/3}}\sigma {T}_{\mathrm{eq}}^{4}.\end{eqnarray} \tag{ 8 }$$

There is no guarantee that Equation (8) will be exact when the factor preceding $\sigma {T}_{\mathrm{eq}}^{4}$ is of order unity. Nevertheless, because the expression reduces to the correct limits for both a thin and a thick atmosphere, it is physically motivated.

To transform F_n into a prediction for the dayside flux F_d , one can use the planet’s energy budget, $1/2({F}_{d}+{F}_{n})=\sigma {T}_{\mathrm{eq}}^{4}$, so ${F}_{d}=2\sigma {T}_{\mathrm{eq}}^{4}-{F}_{n}$. However, this argument does not yet account for the observer’s skewed viewing geometry. At secondary eclipse, the observer’s view is weighted toward the hot substellar point, which directly faces the observer, and is less sensitive to colder regions that lie close to the terminator, which tilt away from the observer (Cowan & Agol 2008). The limiting expressions for the observed dayside flux are given by Equation (1), namely ${F}_{d,\mathrm{obs}}=\sigma {T}_{\mathrm{day}}^{4}=8/3\times \sigma {T}_{\mathrm{eq}}^{4}$ for a bare rock, and ${F}_{d,\mathrm{obs}}=\sigma {T}_{\mathrm{eq}}^{4}$ for a planet with uniform heat redistribution. Since the scaling that interpolates between the bare-rock and uniform limits on the nightside should also apply on the dayside, the observed dayside flux is thus

$$\begin{eqnarray}&&{F}_{d,\mathrm{obs}}=\left(\displaystyle \frac{8}{3}-\displaystyle \frac{5}{3}\displaystyle \frac{{\tau }_{\mathrm{LW}}^{1/3}{\left(\tfrac{{p}_{s}}{1\,\mathrm{bar}}\right)}^{2/3}{\left(\tfrac{{T}_{\mathrm{eq}}}{600\,{\rm{K}}}\right)}^{-4/3}}{k+{\tau }_{\mathrm{LW}}^{1/3}{\left(\tfrac{{p}_{s}}{1\,\mathrm{bar}}\right)}^{2/3}{\left(\tfrac{{T}_{\mathrm{eq}}}{600\,{\rm{K}}}\right)}^{-4/3}}\right)\sigma {T}_{\mathrm{eq}}^{4}.\end{eqnarray} \tag{ 9 }$$

Here, $k={Lg}/{(\chi \beta {c}_{p})\times ({C}_{d}{\sigma }^{2}/R)}^{1/3}{(1\,\mathrm{bar})}^{-2/3}{(600\,{\rm{K}})}^{4/3}$ captures all planetary parameters other than the optical thickness, surface pressure, and equilibrium temperature. Detailed evaluation below shows that, as long as R and c_p correspond to high mean-molecular-weight (MMW) atmospheres, k is roughly of order unity and tends to vary little between common planetary scenarios (e.g., radius and surface gravity tend to vary by less than a factor of two between different rocky planets). We therefore separate out this parameter to underline the dominant dependency of heat redistribution on p_s , τ_LW, and T_eq.

To relate the main result back to the heat redistribution factor f used in 1D models, one simply has to divide the factor preceding $\sigma {T}_{\mathrm{eq}}^{4}$ in Equation (9) by four (see Equation (1)):

$$\begin{eqnarray}&&f=\displaystyle \frac{2}{3}-\displaystyle \frac{5}{12}\times \displaystyle \frac{{\tau }_{\mathrm{LW}}^{1/3}{\left(\tfrac{{p}_{s}}{1\,\mathrm{bar}}\right)}^{2/3}{\left(\tfrac{{T}_{\mathrm{eq}}}{600\,{\rm{K}}}\right)}^{-4/3}}{k+{\tau }_{\mathrm{LW}}^{1/3}{\left(\tfrac{{p}_{s}}{1\,\mathrm{bar}}\right)}^{2/3}{\left(\tfrac{{T}_{\mathrm{eq}}}{600\,{\rm{K}}}\right)}^{-4/3}}.\end{eqnarray} \tag{ 10 }$$

As expected, Equation (10) recovers the no-redistribution limit f → 2/3 as τ_LW, p_s → 0, and the uniform-redistribution limit f → 1/4 as τ_LW, p_s → ∞ .

3.1. Numerical Setup

3.2. GCM Results

Figure 1 shows that thermal phase curves start to significantly deviate from a bare rock once surface pressure exceeds about 1 bar, in agreement with the analytical scaling. Shown are results for simulations with constant optical thickness, τ_LW = 1. The thermal phase curves are computed across the JWST MIRI bandpass (5–12 μm). To convert a planet’s simulated thermal emission to a planet–star flux ratio, the stellar spectrum is additionally assumed to be a blackbody. At surface pressures of 0.1 bar or less, phase curves are effectively indistinguishable from a bare rock. This is particularly the case once JWST's instrumental precision is included, which cannot be any better than the photon noise limit but could optimistically be comparable to it (Figure 1). Conversely, at surface pressures above 10 bar, phase curves become increasingly uniform. This matches the theoretical prediction expressed in Equation (9), namely that as long as k is of order unity, the transition between a bare rock and a uniform planet also happens at ${ \mathcal O }(1)$ bar. Moreover, for a fixed surface pressure, LHS3844b has the least efficient redistribution while TRAPPIST-1b has the most efficient redistribution. This matches the expectation that, at fixed surface pressure and optical thickness, redistribution becomes less efficient on hotter planets.

Zoom In Zoom Out Reset image size

Figure 2 compares the theoretical prediction against the dayside brightness temperatures that an observer would see at secondary eclipse. To evaluate k in Equation (10), we assume a high-MMW atmosphere with $(R,{c}_{p})=({R}_{{{\rm{N}}}_{2}},{c}_{p,{{\rm{N}}}_{2}})$, C_d = 1.9 × 10⁻³, and L = a. To constrain χ, note that previous work showed some GCM simulations of tidally locked planets can indeed approach the ideal limit described by Carnot’s theorem; more typical simulations, however, find wind speeds that fall a factor of two to four below those predicted for an ideal heat engine (Koll & Abbot 2016). The underlying reason is that other atmospheric processes, such as small-scale convection, also redistribute heat but do so without generating large-scale kinetic energy. Based on comparison with the GCM experiments, this paper adopts χ = 0.2, but the value should be treated as uncertain by at least a factor of two since χ accounts not only for the atmosphere’s thermodynamic efficiency but also the circulation geometry.

Zoom In Zoom Out Reset image size

Adopting the above values, k only depends weakly on a planet’s exact properties, with k = 1.2 for TRAPPIST-1b, k = 1.9 for GJ1132b, and k = 2.3 for LHS3844b. At least across the range of planets considered here, the factors determining k (e.g., planetary radius) are thus of secondary importance for a planet’s heat redistribution compared to variations in p_s , τ_LW, and T_eq. The notable exception to this rule would be low-MMW atmospheres. Because k ∝ 1/(c_p R^1/3) ∼ MMW^4/3, the scaling predicts that planets with H₂-dominated atmospheres will have significantly lower k, and thus more efficient heat redistribution, than planets with CO₂-, N₂-, or H₂O-dominated atmospheres. This result matches previous theoretical work on exoplanet heat transport which also found that heat redistribution is significantly more efficient in low-MMW atmospheres (Heng & Kopparla 2012; Menou 2012; Zhang & Showman 2017).

Figure 2 shows that the scaling correctly captures the main variation in dayside thermal emission across variations of several orders of magnitude in surface pressure, and more than one order of magnitude variation in optical thickness. The scaling is not perfect, however. For example, the scaling tends to slightly underpredict thermal emission at low optical thickness, while it overpredicts thermal emission at high optical thickness. This suggests that day–night redistribution is somewhat more sensitive to changes in τ_LW than is implied by the scaling, which in the optically thin limit predicts ${F}_{n}\propto {\tau }_{\mathrm{LW}}^{1/3}$ (Equation (7)). There are also significant differences between scaling and simulations for atmospheres with high surface pressure and low optical thickness, in particular TRAPPIST-1b with τ_LW = 0.1 and p_s > 10 bar. Although the combination of such high surface pressure and low optical thickness is probably unrealistic (see Section 4), the scaling becomes inaccurate in this regime. One potential reason is that, at high surface pressure, convective heat fluxes should become important. Convection implies that the atmosphere ceases being isothermal and instead tends toward a vertically adiabatic profile, which would reduce the efficiency with which adiabatic compression can heat the nightside (see Equation (2)).

Nevertheless, the majority of simulations in Figure 2 fall inside the red envelope around the theoretical scaling. This indicates that the differences between scaling and simulations can be largely explained by a factor-of-two inaccuracy in the derivation. One potential source of inaccuracy is the inefficiency factor χ, which the scaling assumes to be constant but which could vary if the atmosphere’s heat engine efficiency shifted in response to changes in surface pressure or optical thickness. Another potential source of inaccuracy is the effect of planetary rotation, which is not captured by the scaling. For a given planet, the rotation rate is constant, so one might naively expect a roughly constant mismatch between scaling and simulations. A planet’s rotation interacts with the large-scale circulation, however, to induce spatially varying flow and temperature patterns that could shift with changes in pressure or optical thickness and thus create a varying mismatch between scaling and simulations.

Figure 3 shows that the scaling also captures the main trends in nightside brightness temperatures. This is not surprising—given that the scaling largely captures the dayside’s thermal emission, energy conservation implies that the scaling should match the nightside thermal emission with similar accuracy. As on the dayside, the majority of nightside brightness temperatures agree with the scaling to within a factor-of-two change in k. This means the scaling is able to roughly capture how changes in atmospheric and planetary properties affect a planet’s dayside and nightside thermal emission, which is important because one or both of these quantities can be directly measured via secondary eclipse and thermal phase curve observations.

Zoom In Zoom Out Reset image size

Up to now, the derivation has assumed gray radiative transfer, but real gases have an optical thickness that varies strongly as a function of wavelength. How can one assign an equivalent optical thickness τ_LW to atmospheres with realistic compositions? We tested Rosseland and Planck mean opacities but found that they provide a poor fit. Instead, this section will focus on the degree to which an atmosphere diminishes the surface’s thermal emission and thus reduces a planet’s observable surface emission below that of a bare rock.

For an atmosphere with spectrally varying absorption, the planet’s top-of-atmosphere thermal flux is equal to

$$\begin{eqnarray}&&F=\pi \int {B}_{\lambda }({T}_{s}){e}^{-{\tau }_{\lambda }}d\lambda +\left(\mathrm{atm}.\,\mathrm{emission}\right).\end{eqnarray} \tag{ 11 }$$

Here, λ is wavelength, T_s is the surface temperature, the first term is the surface’s blackbody emission attenuated by the overlying atmosphere, and the second term is the atmosphere’s emission. For comparison, with a gray absorber, the top-of-atmosphere thermal flux is equal to

$$\begin{eqnarray}&&F=\sigma {T}_{s}^{4}{e}^{-{\tau }_{\mathrm{LW}}}+\left(\mathrm{atm}.\,\mathrm{emission}\right),\end{eqnarray} \tag{ 12 }$$

where τ_LW is now the atmosphere’s gray optical thickness.

By setting the two surface terms equal to each other, one can define the equivalent gray optical thickness for any atmospheric composition as

$$\begin{eqnarray}&&{\tau }_{\mathrm{LW}}\equiv -\mathrm{ln}\left[\displaystyle \frac{\int {e}^{-{\tau }_{\lambda }}{B}_{\lambda }({T}_{s})d\lambda }{\int {B}_{\lambda }({T}_{s})d\lambda }\right].\end{eqnarray} \tag{ 13 }$$

This equivalent gray optical thickness has the desirable property that it exactly matches the extent to which an atmospheric column attenuates the surface’s thermal flux. It is therefore most appropriate for the thin-atmosphere limit, in which most of a planet’s thermal emission originates at the surface, and it correctly captures how the gradual increase of atmospheric mass then reduces the planet’s surface thermal emission below that of a bare rock.

Figure 4 shows τ_LW for four representative atmospheric compositions, namely pure CO₂, H₂O, O₃, and CO. To calculate τ_LW, a 1D radiative transfer model with line-by-line spectral resolution was used (Koll & Cronin 2018). ² The surface temperature is set to the equilibrium temperature of LHS3844b, and the overlying atmosphere follows a dry adiabat. All molecular opacities are taken from the HITRAN2016 database (Gordon et al. 2017). Collision-induced absorption (CIA) is included for H₂O and CO₂ using the fits from Pierrehumbert (2010). The HITRAN2016 database does not include weak lines, which become important opacity sources at high temperatures. Similarly, we are not aware of CO–CO or O₃–O₃ CIA data, even though CIA should become important at high surface pressures. The optical thicknesses shown in Figure 4 should therefore be treated as a lower bound.

Zoom In Zoom Out Reset image size

H₂O is a very strong absorber, with 10⁻² bar of H₂O enough to make the atmosphere optically thick. Next follow CO₂ and O₃, which become optically thick once surface pressure exceeds about 1 bar. Finally, CO is the poorest absorber out of these gases, and is optically thin even at a surface pressure of 1 bar. The rate at which τ_LW increases with pressure can be roughly understood as follows. At low surface pressure, τ_LW is dominated by the behavior of individual molecular absorption lines, so absorption is proportional to pressure (see strong and weak line limits in Pierrehumbert 2010). H₂O and CO both follow this expectation below about 10⁻² bar. At high surface pressures, absorption becomes dominated by CIA, which scales with pressure squared. H₂O follows this expectation above 1 bar. Finally, for CO₂ and O₃, as well as CO above 0.1 bar, τ_LW is set by the closing of a molecule’s window regions and thus a molecule’s overall band shape. For a molecular band for which the absorption cross-section decays roughly exponentially away from the band center, this implies a logarithmic scaling with pressure (Pierrehumbert 2010; Koll & Cronin 2018).

Figure 4 shows that, for an atmosphere with a moderately strong greenhouse effect, such as CO₂, the equivalent optical thickness of a 1 bar atmosphere is of order unity. Combined with the heat transport scaling (Equation (10)), and assuming a high MMW atmosphere for which k is of order unity, this implies that a typical surface pressure above which the thermal phase curve of a rocky exoplanet with T_eq ∼ 600 K starts to deviate from that of a bare rock is also about 1 bar. However, the scaling underlines that there are clear exceptions to this rule of thumb: in addition to low-MMW atmospheres, for which k becomes much smaller than unity, atmospheres with a particularly strong or weak greenhouse effect will also exhibit a more or less efficient heat transport (e.g., steam or pure CO atmospheres).

The scaling proposed here provides first-order insight into the processes that determine the day–night thermal contrasts of tidally locked, rocky exoplanets. The scaling builds on the numerical results of Koll & Abbot (2015) and the theory of Koll & Abbot (2016), and translates them into observable quantities that can be measured in the near future with JWST via secondary eclipse and phase curve measurements.

Although the scaling relies on gray radiation for its derivation, it appears to agree with previous GCM studies that used comprehensive radiative transfer. Selsis et al. (2011) found that the transition between large day–night contrasts and uniform emission on short-period rocky planets with CO₂ atmospheres occurs somewhere between 0.1 and 1 bars. The scaling here agrees with their results and predicts a transition at 0.4 bar, assuming their planetary parameters and τ_LW ∼ 1. Similarly, Yang et al. (2013) and Wolf et al. (2019) simulated habitable planets with 1 bar N₂–H₂O atmospheres. These studies found thermal phase curves that are fairly uniform relative to a bare rock, and thus also support the scaling here, which predicts a transition at about 0.2 bar. Note that the last two studies simulated habitable planets with a full hydrologic cycle, which allows the atmosphere to transport more heat than predicted by the scaling (see below). Habitable planets might thus exhibit efficient heat redistribution at surface pressures even below 0.2 bar.

The scaling still contains some fundamental uncertainties. In addition to the inefficiency factor χ, the derivation relies on an Ansatz to generalize from thin atmospheres to arbitrary surface pressures. Furthermore, the derivation does not account for the effects of planetary rotation. Figures 2 and 3 show that the basic approach is broadly justified, but also suggest that the match with the GCM simulations tends to be better at low surface pressures than at high pressures. Future analytical work could therefore focus on the limit of extremely thick atmospheres, and consider a wider range of planetary rotation rates, to improve the scaling in those regimes.

Another outstanding issue is that the thermal emission of rocky exoplanets can be affected by additional atmospheric and planetary physics not included here. Broadly speaking, we expect that the scaling should overpredict a planet’s day–night thermal contrast. First, the scaling assumes dry thermodynamics, but an atmosphere with condensation can additionally transport heat via latent heat transport (Ding & Pierrehumbert 2018). The derivation here does not capture this process, so it should tend to overpredict day–night thermal contrasts in atmospheres in which condensation starts to become important. Second, some exoplanets might be able to sustain a surface ocean, in particular planets inside their host star’s habitable zone. In addition to supporting a hydrologic cycle and latent heat transport, an ocean can additionally transport energy via its own fluid motions, which will again reduce the planet’s day–night thermal contrast (Yang et al. 2019).

Finally, the scaling proposed here does not account for the radiative effect of clouds, but it is unclear whether clouds should systematically increase or decrease a planet’s day–night emission contrast. For clouds that are linked to atmospheric convection, such as H₂O clouds on habitable planets, clouds will occur preferentially on a planet’s dayside and thus reduce the planet’s day–night emission contrast (Yang et al. 2013). However, on extremely hot planets with large day–night temperature contrasts, such as 55 Cnc e, clouds can also form by advection of Na or SiO vapor to the planet’s cold nightside (Hammond & Pierrehumbert 2017). In this case, clouds would preferentially amplify the planet’s day–night emission contrast.

This paper proposes an analytical scaling for atmospheric heat redistribution on tidally locked rocky exoplanets. The surface pressure at which a planet’s broadband secondary eclipse and thermal phase curve transitions from inefficient to efficient heat redistribution depends on the atmosphere’s bulk composition, its radiative properties, and the planet’s equilibrium temperature.

For “warm” planets like TRAPPIST-1b, GJ1214b, and LHS3844b, a typical threshold pressure is ${ \mathcal O }(1)$ bar, but there are notable exceptions. H₂-dominated atmospheres redistribute heat much more efficiently than high-MMW atmospheres. Heat redistribution also depends on optical thickness, with optically thick H₂O atmospheres better at redistributing heat than optically thin CO atmospheres. Finally, JWST will also be able to observe thermal emission from hot “lava ocean” worlds, such as 55 Cancri e and K2-141b. Because heat redistribution is less efficient on hotter planets, the scaling predicts that these planets require atmospheres thicker than ${ \mathcal O }(1)$ bar surface pressure to exhibit significant heat redistribution.

In addition to improving our understanding of atmospheric heat transport, the scaling is useful in the context of interpreting thermal observations of rocky exoplanets. For favorable targets, JWST will be able to measure broadband thermal fluxes as well as emission spectra (e.g., Greene et al. 2016; Morley et al. 2017; Batalha et al. 2018; Kempton et al. 2018). To make sense of these measurements, models are needed that properly account for the physical processes that shape the observable features of rocky exoplanets. The scaling derived here offers a way of parameterizing the large-scale atmospheric heat transport in 1D models, and was used as such in a number of related studies (Koll et al. 2019; Kreidberg et al. 2019; Malik et al. 2019; Mansfield et al. 2019). Doing so is attractive because parameterization is many orders of magnitude faster than explicitly simulating an atmosphere’s large-scale fluid dynamics.

This work was supported by a James McDonnell Foundation postdoctoral fellowship. The paper was improved by constructive comments from an anonymous reviewer.