Why Is it So Cold in Here? Explaining the Cold Temperatures Retrieved from Transmission Spectra of Exoplanet Atmospheres

A transmission spectrum is given by the wavelength-dependent area ratio between a transiting planet and its star

$$\begin{eqnarray}&&{{\rm{\Delta }}}_{\lambda }=\displaystyle \frac{\pi {R}_{{\rm{p}}}^{2}+{\int }_{0}^{2\pi }{\int }_{{R}_{{\rm{p}}}}^{\infty }\left(1-{e}^{-{\tau }_{\lambda }(b,\theta )}\right)\,b\,{db}\,d\theta }{\pi {R}_{* }^{2}}\end{eqnarray} \tag{ 1 }$$

where ${R}_{* }$ and ${R}_{{\rm{p}}}$ are the stellar and planetary radii, respectively, b is the impact parameter, θ is the azimuthal angle, and ${\tau }_{\lambda }(b,\theta )$ is the slant optical depth—the extinction coefficient integrated along the line of sight. The first term represents the disk area of the opaque deep atmosphere at a reference pressure $P(r={R}_{{\rm{p}}})={P}_{0}$. The second term gives the effective area of successive atmospheric elements in a polar coordinate system, weighted by the absorptivity of each element.

For a 1D atmosphere, this expression reduces to

$$\begin{eqnarray}&&{{\rm{\Delta }}}_{\lambda ,1{\rm{D}}}=\displaystyle \frac{{R}_{{\rm{p}}}^{2}+2{\int }_{{R}_{{\rm{p}}}}^{\infty }\left(1-{e}^{-{\tau }_{\lambda }(b)}\right)\,b\,{db}}{{R}_{* }^{2}}.\end{eqnarray} \tag{ 2 }$$

Analytical tractability arises via the following assumptions: (i) constant pressure scale height with altitude (i.e., an isothermal, isocompositional, isogravitational atmosphere); (ii) hydrostatic equilibrium and the ideal gas law hold; (iii) cross sections vary weakly with pressure (i.e., ${\sigma }_{\lambda }(P,T)\approx {\sigma }_{\lambda }(T)$); (iv) atmospheric altitudes satisfy $z/{R}_{p}\ll 1$. With these assumptions, it is well established that Equation (2) can be simply written as (e.g., Fortney 2005; Lecavelier des Etangs et al. 2008; de Wit & Seager 2013; Bétrémieux & Swain 2017; Heng & Kitzmann 2017)

$$\begin{eqnarray}&&{{\rm{\Delta }}}_{\lambda ,1{\rm{D}}}=\displaystyle \frac{{R}_{{\rm{p}}}^{2}+2{R}_{{\rm{p}}}{H}_{1{\rm{D}}}(\gamma +\mathrm{ln}{\tau }_{0,\lambda ,1{\rm{D}}})}{{R}_{* }^{2}}\end{eqnarray} \tag{ 3 }$$

where ${H}_{1{\rm{D}}}={{kT}}_{1{\rm{D}}}/\mu g$ is the scale height of the 1D atmosphere (T_1D, μ, g, and k being respectively the 1D temperature, mean molecular mass, surface gravity, and Boltzmann constant), $\gamma \approx 0.57722$ is the Euler–Mascheroni constant, and ${\tau }_{0,\lambda ,1{\rm{D}}}$ is given by¹

$$\begin{eqnarray}&&{\tau }_{0,\lambda ,1{\rm{D}}}=\displaystyle \frac{{P}_{0}}{{{kT}}_{1{\rm{D}}}}\sqrt{2\pi {R}_{{\rm{p}}}{H}_{1{\rm{D}}}}\sum _{i}{X}_{1{\rm{D}},i}\,{\sigma }_{\lambda ,i}({T}_{1{\rm{D}}})\end{eqnarray} \tag{ 4 }$$

where ${X}_{1{\rm{D}},i}$ and ${\sigma }_{\lambda ,i}$ are the (1D) volume-mixing ratio and absorption cross section of chemical species i, respectively.

For the 2D atmosphere depicted in Figure 1, we can write Equation (1) as

$$\begin{eqnarray}\begin{array}{rcl}{{\rm{\Delta }}}_{\lambda ,2{\rm{D}}} & = & \displaystyle \frac{{R}_{{\rm{p}}}^{2}+{\int }_{{R}_{{\rm{p}}}}^{\infty }\left(1-{e}^{-{\tau }_{\lambda ,{\rm{M}}}(b)}\right)\,b\,{db}+{\int }_{{R}_{{\rm{p}}}}^{\infty }\left(1-{e}^{-{\tau }_{\lambda ,{\rm{E}}}(b)}\right)\,b\,{db}}{{R}_{* }^{2}}\end{array}\end{eqnarray} \tag{ 5 }$$

where "E" and "M" denote the evening and morning terminators. To clearly isolate the effect of interest for the present study, we assume the atmosphere to be uniform within each terminator sector, along the day–night slant path (see Caldas et al. 2019), and a morning–evening transition with negligible width (i.e., "M" and "E" both span ${\rm{\Delta }}\theta =\pi$). Making the same simplifying assumptions as the 1D case (but now with two isotherms, T_E and T_M), one obtains²

$$\begin{eqnarray}&&{{\rm{\Delta }}}_{\lambda ,2{\rm{D}}}=\displaystyle \frac{{R}_{{\rm{p}}}^{2}+{R}_{{\rm{p}}}{H}_{{\rm{M}}}(\gamma +\mathrm{ln}{\tau }_{0,\lambda ,{\rm{M}}})+{R}_{{\rm{p}}}{H}_{{\rm{E}}}(\gamma +\mathrm{ln}{\tau }_{0,\lambda ,{\rm{E}}})}{{R}_{* }^{2}}.\end{eqnarray} \tag{ 6 }$$

Fitting a 1D model to a 2D spectrum is equivalent to setting ${{\rm{\Delta }}}_{\lambda ,2{\rm{D}}}={{\rm{\Delta }}}_{\lambda ,1{\rm{D}}}$. One may expect this condition to result in equivalent 1D properties given by ${T}_{1{\rm{D}}}=\bar{T}\equiv \tfrac{1}{2}({T}_{{\rm{E}}}+{T}_{{\rm{M}}})$ and ${X}_{1{\rm{D}},i}=\bar{{X}_{i}}\equiv \tfrac{1}{2}({X}_{{\rm{E}},i}+{X}_{{\rm{M}},i})$—i.e., terminator average temperature and mixing ratios. However, this is not the case. In Appendix A, we derive that the actual equivalent 1D temperature is³

$$\begin{eqnarray}&&{T}_{1{\rm{D}}}=\bar{T}\left[\displaystyle \frac{{{\rm{\Psi }}}_{\lambda }}{{W}_{-1}\left({{\rm{\Psi }}}_{\lambda }{\left(\bar{X}/{X}_{1{\rm{D}}}\right)}^{2}{e}^{-2(\gamma +\mathrm{ln}{\bar{\tau }}_{0,\lambda })}\right)}\right]\end{eqnarray} \tag{ 7 }$$

where

$$\begin{eqnarray}&&{\bar{\tau }}_{0,\lambda }\,=\displaystyle \frac{{P}_{0}}{k\bar{T}}\sqrt{2\pi {R}_{{\rm{p}}}\bar{H}}\bar{X}{\sigma }_{\lambda }\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&{{\rm{\Psi }}}_{\lambda }=-2\left[\gamma +\mathrm{ln}{\bar{\tau }}_{0,\lambda }+\tilde{f}\left(\displaystyle \frac{{\rm{\Delta }}T}{\bar{T}}\right)+\tilde{g}\left(\displaystyle \frac{{\rm{\Delta }}T}{\bar{T}},{\rm{\Delta }}\mathrm{ln}X\right)\right]\end{eqnarray} \tag{ 9 }$$

and ${W}_{-1}(x)$ is the lower real branch of the Lambert W function.⁴ $\tilde{f}$ and $\tilde{g}$ are dimensionless functions of the temperature and compositional differences between the evening and morning terminators: ${\rm{\Delta }}T\equiv \tfrac{1}{2}({T}_{{\rm{E}}}-{T}_{{\rm{M}}})$, ${\rm{\Delta }}\mathrm{ln}X\equiv \tfrac{1}{2}(\mathrm{ln}{X}_{{\rm{E}}}-\mathrm{ln}{X}_{{\rm{M}}})$. They are given by

$$\begin{eqnarray}\begin{array}{rcl}\tilde{f} & = & -\displaystyle \frac{1}{4}\left\{\left(1-\displaystyle \frac{{\rm{\Delta }}T}{\bar{T}}\right)\mathrm{ln}\left(1-\displaystyle \frac{{\rm{\Delta }}T}{\bar{T}}\right)\right.\\ & & \left.+\left(1+\displaystyle \frac{{\rm{\Delta }}T}{\bar{T}}\right)\mathrm{ln}\left(1+\displaystyle \frac{{\rm{\Delta }}T}{\bar{T}}\right)\right\}\end{array}\end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray}&&\tilde{g}=\left(\displaystyle \frac{{\rm{\Delta }}T}{\bar{T}}\right){\rm{\Delta }}\mathrm{ln}X-\mathrm{ln}\left[\cosh ({\rm{\Delta }}\mathrm{ln}X)\right].\end{eqnarray} \tag{ 11 }$$

We now demonstrate that Equation (7) predicts ${T}_{1{\rm{D}}}\lt \bar{T}$ readily occurs for ${\rm{\Delta }}\mathrm{ln}X\ne 0$ when ${X}_{1{\rm{D}}}=\bar{X}$ (i.e., assuming mixing ratios are correctly retrieved—we revisit this in Section 2.3).

Equation (7) is graphically rendered in Figure 2 for a typical hot Jupiter ($\bar{T}=1200$ K) over a range of ${\bar{\tau }}_{0,\lambda }$. Specifically, we show three ${\bar{\tau }}_{0,\lambda }$ surfaces considering different potential opacity sources for a hot Jupiter with ${R}_{{\rm{p}}}=1.4{R}_{{\rm{J}}}$, g = 10 m s⁻², and $\mu =2.3\,{m}_{{\rm{u}}}$. The first surface (${\bar{\tau }}_{0,\lambda }={10}^{4}$) corresponds to H₂O absorption around $1.4\,\mu {\rm{m}}$ (${\sigma }_{\lambda }\sim {10}^{-21}$ cm²; Sharp & Burrows 2007) with ${X}_{{{\rm{H}}}_{2}{\rm{O}}}={10}^{-4}$. The second surface (${\bar{\tau }}_{0,\lambda }={10}^{6}$) corresponds to TiO absorption around $0.7\,\mu {\rm{m}}$ (${\sigma }_{\lambda }\sim {10}^{-16}$ cm²; Sharp & Burrows 2007) with ${X}_{\mathrm{TiO}}={10}^{-7}$. Finally, the third surface (${\bar{\tau }}_{0,\lambda }={10}^{8}$) corresponds to Na doublet line core absorption around $0.6\,\mu {\rm{m}}$ (${\sigma }_{\lambda }\sim {10}^{-15}$ cm²; Allard et al. 2019) with ${X}_{\mathrm{Na}}={10}^{-6}$. From examining Figure 2, one notices three important takeaways:

• 1.  
  pure temperature differences (${\rm{\Delta }}\mathrm{ln}X=0$) result in negligible biases to retrieved temperatures;
• 2.  
  compositional differences exceeding a factor of 2 (${\rm{\Delta }}{\mathrm{log}}_{10}(X)\gt 0.3$) result in T_1D biases to many hundreds of degrees colder than $\bar{T}$;
• 3.  
  the wavelength dependency of T_1D (from ${\bar{\tau }}_{0,\lambda }$) implies that no one equivalent temperature can perfectly reproduce a 2D spectrum using a 1D model.

Observations 1 and 2 arise from properties of $\tilde{f}$ and $\tilde{g}$, elucidated in Appendix B. The final observation yields an important consequence for fitting 2D transmission spectra with 1D models: as a retrieval can only chose one value of ${T}_{1{\rm{D}}}\,\forall \lambda$, the chosen value will balance the different ${\bar{\tau }}_{0,\lambda }$ surfaces to minimize residuals between the 2D data and 1D model. The chosen T_1D, considered as a wavelength-average, will however still be biased to colder temperatures than $\bar{T}$ (Figure 2). It may be possible to use such residuals, with sufficiently precise observations, to identify a preference for 2D atmospheric models.

The conceptual picture provided by Equation (7) and Figure 2 will be altered in regimes where the underlying assumptions break down. We highlight two important limitations: (i) if T_1D becomes sufficiently cooler than $\bar{T}$, neglected temperature dependencies in absorption cross sections, ${\sigma }_{\lambda }(T)$, will alter the shapes of spectral features; (ii) as ${\rm{\Delta }}\mathrm{ln}X$ grows, the assumption that only one species dominates the opacity will be violated on the terminator side deficient in the given species. Taken together, the breakdown of these assumptions will place a lower limit on how cold T_1D may become.

1D models have two additional degrees of freedom to compensate for such higher-order effects: X_1D and ${R}_{{\rm{p}}}$. Varying the former is already encapsulated by Equation (7). Varying the latter corresponds to the retrieved base planet radius, ${R}_{{\rm{p}},1{\rm{D}}}$, differing from the actual base radius, ${R}_{{\rm{p}},2{\rm{D}}}$ (equated until now). A generalization of Equation (7) for ${R}_{{\rm{p}},1{\rm{D}}}\ne {R}_{{\rm{p}},2{\rm{D}}}$ is presented in Appendix A. Perturbing either X_1D or ${R}_{p,1{\rm{D}}}$ (therefore biasing these quantities) effectively translates the surfaces in Figure 2 along the T_1D axis (not shown). We thus expect that some of the T_1D bias will be "shifted" into X_1D and ${R}_{{\rm{p}},1{\rm{D}}}$, each attaining their own bias.

To relax many of the aforementioned assumptions, and establish the extent of T_1D, X_1D, and ${R}_{{\rm{p}},1{\rm{D}}}$ biases, we turn to more physically realistic numerical models.

As the cold retrieved temperatures of exoplanets span a wide range of equilibrium temperatures (Table 1), so too must our proposed explanation. We therefore consider three diverse case studies of atmospheres expected to posses morning–evening compositional differences.

• 1.  
  Warm Jupiter: ${\bar{T}}_{1\mathrm{mbar}}\sim 1000$ K, ${\rm{\Delta }}T=100$ K (i.e., a 200 K morning–evening temperature difference). The warmer (evening) terminator has Na and K abundances representative of solar elemental abundances: log(${X}_{{\rm{E}},\mathrm{Na}}$) = −6, log(${X}_{{\rm{E}},{\rm{K}}}$) = −7 (Asplund et al. 2009). The cooler (morning) terminator is assumed depleted in Na and K by two orders of magnitude (a proxy for condensation). CH₄ roughly follows equilibrium abundances for a solar-composition atmosphere: log(${X}_{{\rm{E}},{\mathrm{CH}}_{4}}$) = −6, log(${X}_{{\rm{M}},{\mathrm{CH}}_{4}}$) = −4 (Heng & Tsai 2016). H₂O takes a solar abundance, assumed constant around the terminator: log(${X}_{{\rm{E}},{{\rm{H}}}_{2}{\rm{O}}}$) = log(${X}_{{\rm{M}},{{\rm{H}}}_{2}{\rm{O}}}$) = −3.3.
• 2.  
  Hot Jupiter: ${\bar{T}}_{1\mathrm{mbar}}\sim 1600$ K, ${\rm{\Delta }}T=100$ K. Both terminators possess constant Na, K, and H₂O abundances: log(${X}_{\mathrm{Na}}$) = −6, log(${X}_{{\rm{K}}}$) = −7, log(${X}_{{{\rm{H}}}_{2}{\rm{O}}}$) = −3.3. The evening terminator additionally contains TiO and VO with roughly solar abundances: log(${X}_{{\rm{E}},\mathrm{TiO}}$) = −7, log(${X}_{{\rm{E}},\mathrm{VO}}$) = −8. The morning terminator is assumed sufficiently cool for all TiO and VO to have condensed out of the gas-phase in the observable atmosphere.
• 3.  
  Ultra-hot Jupiter: ${\bar{T}}_{1\mathrm{mbar}}\sim 2200$ K, ${\rm{\Delta }}T=250$ K. Both terminators possess Na, K, and H₂O abundances as in the hot Jupiter case. However, here the evening terminator is warm enough for H₂ to partially dissociate and form an inventory of the hydrogen anion (Parmentier et al. 2018), for which we take: log(${X}_{{\rm{E}},{{\rm{H}}}^{-}}$) = −8. The morning terminator is assumed too cold to support ${{\rm{H}}}^{-}$.

In all three cases, physical properties are representative of HD 209458b: ${R}_{{\rm{p}}}=1.359{R}_{{\rm{J}}}$,⁵ ${M}_{{\rm{p}}}=0.6845{M}_{{\rm{J}}}$. Each planet is assumed H₂+He dominated, with a solar-proportion He/H₂ ratio of 0.17. The pressure–temperature (P–T) profiles and morning–evening temperature differences are inspired by literature GCM profiles (e.g., Kataria et al. 2016; Helling et al. 2019), constructed parametrically (Madhusudhan & Seager 2009).⁶ Each profile has an "anchor" temperature at 10 bar, ${T}_{\mathrm{deep}}$, below which the atmosphere is homogeneous: ${T}_{\mathrm{deep},\mathrm{warm}}=1600$ K, ${T}_{\mathrm{deep},\mathrm{hot}}\,=2200$ K, ${T}_{\mathrm{deep},\mathrm{ultra} \mbox{-} \mathrm{hot}}=3000$ K. To isolate biases arising solely from the compositional and temperature differences, the terminators are assumed cloud-free. We note that the morning–evening differences assumed here are intended as illustrative of their corresponding biases, with strictly self-consistent profiles considered in Section 4.

Our model atmospheres, radiative transfer, and retrievals are computed using the POSEIDON atmospheric retrieval code (MacDonald & Madhusudhan 2017). The atmospheric column in each terminator is discretized uniformly in log-pressure with 10 layers per decade, for 81 levels from 10⁻⁶ to 10² bar. The deep atmosphere ($P\geqslant 10$ bar) has a homogeneous temperature, by construction, hence the terminators share a spherical radial grid below the reference radius ${R}_{{\rm{p}},10\mathrm{bar}}$. Above this, separate radial grids are constructed for each terminator using P–T profiles, abundances, and planetary properties under hydrostatic equilibrium. The warmer evening terminator hence extends to greater radii (see Figure 1). The opacities of Na, K, H₂O, CH₄, TiO, and VO, cross section computations, broadening parameters, and line list references are described in MacDonald (2019). ${{\rm{H}}}^{-}$ continuum opacity is included (John 1988). High-resolution ($R\approx {10}^{6}$) opacities are sampled onto an R = 2000 wavelength grid from 0.3 to 2.0 μm. Radiative transfer is solved separately for each terminator, with 2D transmission spectra constructed by a linear superposition of each terminator spectrum.

Synthetic Hubble observations are generated by convolving each model spectrum to the resolving power of the STIS G430/G750 and WFC3 G141 grisms, before integrating over their respective sensitivity functions. We choose spectral resolutions and precisions typical of current HST observations (e.g., Sing et al. 2016): R = 20 and 100 ppm for STIS and R = 60 and 50 ppm for WFC3. The synthetic data are placed on the true transit depths (i.e., without Gaussian scatter), such that any posterior deviations from the true parameter values are attributable to a retrieval bias rather than a specific noise instance (see Feng et al. 2018). The resulting data are shown in Figure 3 (top row).

Zoom In Zoom Out Reset image size

We subject each synthetic data set to a Bayesian atmospheric retrieval. The retrievals assume a 1D forward model with a single six-parameter P–T profile (Madhusudhan & Seager 2009), a single abundance for each species, and a 10 bar planetary radius. The warm and ultra-hot Jupiter retrievals have 11 free parameters, while the hot Jupiter retrieval has 12. The P–T parameter priors are as described in MacDonald & Madhusudhan (2019), with ${T}_{\mathrm{deep}}$ ascribed a uniform prior from 400–3000 K. The logarithm of each mixing ratio has a uniform prior from −12 to −0.3. ${R}_{{\rm{p}},10\mathrm{bar}}$ has a uniform prior from $0.85\mbox{--}1.15\,{R}_{{\rm{p}}}$. The parameter space is explored via the nested sampling algorithm MultiNest (Feroz & Hobson 2008; Feroz et al. 2009, 2019), implemented by the python wrapper PyMultiNest (Buchner et al. 2014).

Our retrieved spectra, P–T profiles, and mixing ratios are shown⁷ in Figure 3. The 1D models achieve an excellent fit for the warm and hot Jupiters, with the median models lying within 1σ of all data points (i.e., discrepancies <50 ppm). The ultra-hot Jupiter spectral fit is the least accurate, with around 15% of the data incorrectly fit to 1σ. The last observation, most prominent at visible wavelengths, arises from a 1D model dominated by ${{\rm{H}}}^{-}$ being unable to reproduce the superposition of ${{\rm{H}}}^{-}$ and H₂ Rayleigh scattering encoded in the 2D model. However, for real observations with Gaussian scatter, it would be difficult to recognize such a model mismatch.

The retrieved P–T profiles are colder than the terminator-averaged profiles for $P\gt \,1$ mbar. At a 10 mbar reference level (∼the median photosphere), biases of −100, −200, and −1000 K result for the warm, hot, and ultra-hot cases, respectively. The retrieved profiles exhibit shallower temperature gradients than the true profiles, despite the ability of 1D retrievals to retrieve temperature gradients (Rocchetto et al. 2016; MacDonald & Madhusudhan 2017), possibly explaining why many retrieved P–T profiles appear near-isothermal (e.g., MacDonald & Madhusudhan 2019; Pinhas et al. 2019). Our findings match the general trend seen in Table 1: retrieved temperatures become far colder than expectations (compared with skin temperatures) as ${T}_{\mathrm{eq}}$ rises.

The mixing ratios for all significant opacity sources⁸ are incorrectly retrieved to 1σ. Many chemical species are only retrieved accurately to the 3σ level (e.g., CH₄, TiO, and ${{\rm{H}}}^{-}$). Those species exhibiting compositional differences have retrieved 1D abundances biased lower than the true terminator-averaged values. Even species uniform around the terminator (here, H₂O) are biased, though to higher abundances. Compositional biases become more severe as the retrieved P–T profile deviates further from the true average terminator temperature. In the most extreme case, the retrieved H₂O abundance is biased by over an order of magnitude, such that one would incorrectly believe a solar-metallicity atmosphere was $15\,\times$ super-solar at $\gt 3\sigma$ confidence.

We compute self-consistent atmospheric P–T and compositional profiles under the assumptions of radiative–convective and chemical equilibrium with rainout condensation and ionization. An initial P–T profile is iteratively perturbed until a solution satisfying hydrostatic equilibrium, energy conservation, and equilibrium chemistry is obtained. Our models include all the opacities from Goyal et al. (2018), alongside Fe and ${{\rm{H}}}^{-}$. The model is fully described in Goyal (2019) and J. M. Goyal et al. (2020, in preparation). We simulate P–T profiles for each terminator by varying the recirculation factor (e.g., Fortney & Marley 2007)—a 1D proxy for advection due to winds.

Our self-consistent atmospheric profiles are shown in Figure 4. WASP-17b displays large compositional differences in TiO and VO for $P\lesssim {10}^{-2}$ bar due to metal oxide condensation on the cooler morning terminator. We therefore add an opaque cloud deck to the atmospheric model at 10 mbar, serving as a proxy for Ti and V condensates. At the higher temperature of WASP-12b, compositional differences instead arise from H₂ and H₂O dissociation in the upper atmosphere.

Zoom In Zoom Out Reset image size

Transmission spectra are computed for each planet following the same methodology as in Section 3.2. Opacity and mean molecular weight contributions from the altitude-dependent abundances of H₂, He, H, ${{\rm{H}}}^{-}$, Na, K, H₂O, TiO, and VO are considered. For WASP-12b we discount TiO and VO, as their large spectral signatures are inconsistent with current observations at optical wavelengths (Kreidberg et al. 2015; Sing et al. 2016). The resultant 2D spectra are convolved to a similar resolution and precision as literature spectra for each planet (Sing et al. 2016): 200 ppm/400 ppm for WASP-17b (STIS/WFC3) and 100 ppm/50 ppm for WASP-12b. Retrievals are conducted as previously, with the addition of a cloud pressure parameter, ${P}_{\mathrm{cloud}}$, for WASP-17b.

Retrieval results for each planet are shown in Figure 4. We obtain biased retrieved P–T profiles and abundances, despite 1D model spectra providing adequate fits to the 2D atmosphere spectra. Both planets have retrieved P–T profiles biased to colder temperatures, with more extreme biases for the ultra-hot Jupiter WASP-12b.

WASP-17b provides an example where retrieval biases are reasonably contained. Our retrieved photosphere temperature, ${T}_{10\mathrm{mbar}}={1265}_{-237}^{+273}$ K, is biased by ∼250 K below the true terminator-averaged temperature. This is consistent with reported literature temperatures (Pinhas et al. 2019; Welbanks et al. 2019). The retrieved abundances of TiO, VO, and H₂O agree with the true terminator-averaged values to 1σ. This demonstrates that abundances derived from transmission spectra with large error bars (>100 ppm) can still be considered accurate despite 2D terminator differences.

WASP-12b provides a cautionary tale for how 1D retrievals can infer erroneous atmospheric properties. The retrieved temperature, ${T}_{10\mathrm{mbar}}={1711}_{-236}^{+472}$ K, is biased by almost 1000 K below the true terminator-averaged temperature. This is consistent with the cold temperature derived by Kreidberg et al. (2015) to 1σ and that of Pinhas et al. (2019) and Welbanks et al. (2019) to 2σ (Table 1). The abundance posteriors⁹ display a bimodal solution, with a high-metallicity mode (maximum likelihood) and a roughly solar metallicity mode. We represent each solution separately in Figure 4, showing two 1σ error bars derived from separate retrievals with a cut imposed at log(${X}_{{{\rm{H}}}_{2}{\rm{O}}}$) = −2. The preferred mode has ${{\rm{H}}}^{-}$ and H₂O abundances biased by nearly two orders of magnitude: log(${X}_{{{\rm{H}}}^{-}}$) = $-{7.75}_{-0.41}^{+0.31}$ and log(${X}_{{{\rm{H}}}_{2}{\rm{O}}}$) = $-{1.41}_{-0.21}^{+0.17}$, echoing the lesson of Section 3 that 1D abundances for ultra-hot Jupiters must be carefully considered.

In the limit where $\tfrac{{\rm{\Delta }}T}{\bar{T}},{\rm{\Delta }}\mathrm{ln}X\to 0$, we have a uniform, 1D, atmosphere. We therefore expect ${T}_{1{\rm{D}}}\to \bar{T}$. In this limit, Equations (10) and (11) give $\tilde{f},\tilde{g}\to 0$, hence Equation (7) becomes

$$\begin{eqnarray}&&{T}_{1{\rm{D}}}=\bar{T}\left[\displaystyle \frac{{{\rm{\Psi }}}_{\lambda }^{{\prime} }}{{W}_{-1}\left({{\rm{\Psi }}}_{\lambda }^{{\prime} }{e}^{{{\rm{\Psi }}}_{\lambda }^{{\prime} }}\right)}\right]=\bar{T}\end{eqnarray} \tag{ B1 }$$

where ${{\rm{\Psi }}}_{\lambda }^{{\prime} }=-2\left[\gamma +\mathrm{ln}{\bar{\tau }}_{0,\lambda }\right]$ and the last equality uses $W({{xe}}^{x})=x$, following from the definition of the W function.

When ${\rm{\Delta }}\mathrm{ln}X=0$, we have $\tilde{g}=0$; hence Equation (7) becomes

$$\begin{eqnarray}&&{T}_{1{\rm{D}}}=\bar{T}\left[\displaystyle \frac{-2\left[\gamma +\mathrm{ln}{\bar{\tau }}_{0,\lambda }+\tilde{f}\left(\tfrac{{\rm{\Delta }}T}{\bar{T}}\right)\right]}{{W}_{-1}\left(-2\left[\gamma +\mathrm{ln}{\bar{\tau }}_{0,\lambda }+\tilde{f}\left(\tfrac{{\rm{\Delta }}T}{\bar{T}}\right)\right]{e}^{-2(\gamma +\mathrm{ln}{\bar{\tau }}_{0,\lambda })}\right)}\right].\end{eqnarray} \tag{ B2 }$$

Considering a hot Jupiter with $\bar{T}=1400\,{\rm{K}}$ and a typical morning–evening temperature difference ${T}_{{\rm{E}}}-{T}_{{\rm{M}}}\approx 200\,{\rm{K}}$ (Kataria et al. 2016), we have ${\rm{\Delta }}T/\bar{T}\approx 0.07$. As ${\rm{\Delta }}T/\bar{T}\ll 1$, we can Taylor expand Equation (10) to yield

$$\begin{eqnarray}&&\tilde{f}=-\displaystyle \frac{1}{4}\left[{\left(\displaystyle \frac{{\rm{\Delta }}T}{\bar{T}}\right)}^{2}+{ \mathcal O }{\left(\displaystyle \frac{{\rm{\Delta }}T}{\bar{T}}\right)}^{4}+...\,\right]\end{eqnarray} \tag{ B3 }$$

where the odd terms cancel due to symmetry. $\tilde{f}$ is then essentially a minor quadratic perturbing term in Equation (B2). As ${\bar{\tau }}_{0,\lambda }$ is defined in the deep atmosphere at $P={P}_{0}$ (Equation (4)), we have $\gamma +\mathrm{ln}{\bar{\tau }}_{0,\lambda }\gg \tilde{f}$ and hence Equation (B2) tends toward Equation (B1) and ${T}_{1{\rm{D}}}\approx \bar{T}$. Numerical exploration of Equation (B2) yields cooling biases of $\lesssim 5\,$ K, even for ultra-hot Jupiters. We conclude that pure temperature differences have little effect on retrieved temperatures.

To consider the effect of compositional differences, let us write Equation (7) with $\tilde{f}\approx 0$ (as shown in the last section)

$$\begin{eqnarray}&&{T}_{1{\rm{D}}}\approx \bar{T}\left[\displaystyle \frac{-2\left[\gamma +\mathrm{ln}{\bar{\tau }}_{0,\lambda }+\tilde{g}\left(\tfrac{{\rm{\Delta }}T}{\bar{T}},{\rm{\Delta }}\mathrm{ln}X\right)\right]}{{W}_{-1}\left(-2\left[\gamma +\mathrm{ln}{\bar{\tau }}_{0,\lambda }+\tilde{g}\left(\tfrac{{\rm{\Delta }}T}{\bar{T}},{\rm{\Delta }}\mathrm{ln}X\right)\right]{e}^{-2(\gamma +\mathrm{ln}{\bar{\tau }}_{0,\lambda })}\right)}\right].\end{eqnarray} \tag{ B4 }$$

Writing the denominator as ${W}_{-1}(y)$, the exponential factor implies $| y| \ll 0.1$. We can then employ an asymptotic expansion, valid for $-0.1\leqslant y\leqslant 0$, to write (Vazquez-Leal et al. 2019)

$$\begin{eqnarray}&&{W}_{-1}(y)\approx \mathrm{ln}(-y)-\mathrm{ln}(-\mathrm{ln}(-y))+\displaystyle \frac{\mathrm{ln}(-\mathrm{ln}(-y))}{\mathrm{ln}(-y)}.\end{eqnarray} \tag{ B5 }$$

Due to the logarithmic dependence of ${W}_{-1}(y)$ on y, while y itself varies roughly linearly in ${\rm{\Delta }}\mathrm{ln}X$, to a zeroth approximation we can take ${W}_{-1}(y)\approx -\mathrm{const}$. This allows a simpler functional form to be written

$$\begin{eqnarray}&&{T}_{1{\rm{D}}}\mathop{\propto }\limits_{\sim }\bar{T}\left[\gamma +\mathrm{ln}{\bar{\tau }}_{0,\lambda }+\tilde{g}\left(\displaystyle \frac{{\rm{\Delta }}T}{\bar{T}},{\rm{\Delta }}\mathrm{ln}X\right)\right].\end{eqnarray} \tag{ B6 }$$

Taking the limit ${e}^{{\rm{\Delta }}\mathrm{ln}X}\gg {e}^{-{\rm{\Delta }}\mathrm{ln}X}$, for which $\tilde{g}\to \mathrm{ln}2\,-\left(1-\tfrac{{\rm{\Delta }}T}{\bar{T}}\right){\rm{\Delta }}\mathrm{ln}X$, the asymptotic behavior of compositional differences is

$$\begin{eqnarray}&&{T}_{1{\rm{D}}}\mathop{\propto }\limits_{\sim }\bar{T}\left[-\left(1-\displaystyle \frac{{\rm{\Delta }}T}{\bar{T}}\right){\rm{\Delta }}\mathrm{ln}X+\left(\gamma +\mathrm{ln}{\bar{\tau }}_{0,\lambda }+\mathrm{ln}2\right)\right]\end{eqnarray} \tag{ B7 }$$

which explains the roughly linear decrease of T_1D with ${\rm{\Delta }}{\mathrm{log}}_{10}(X)$ shown in Figure 2. The factor of $\left(1-\tfrac{{\rm{\Delta }}T}{\bar{T}}\right)$ modulating the gradient also explains why the coldest values of T_1D occur for pure compositional differences.

Finally, we note that the condition for cooling biases to occur from Equation (B4) can be essentially reduced to $\tilde{g}\lt 0$. However, in the presence of both temperature and compositional differences, Equation (11) has a regime where $\tilde{g}\gt 0$ (the "wrinkle" in Figure 2). We can therefore define a "critical" mixing ratio difference, ${\rm{\Delta }}\mathrm{ln}{X}_{\mathrm{crit}}$, where cooling biases begin according to $\tilde{g}=0$, or

$$\begin{eqnarray}&&\left(\displaystyle \frac{{\rm{\Delta }}T}{\bar{T}}\right){\rm{\Delta }}\mathrm{ln}{X}_{\mathrm{crit}}=\mathrm{ln}\left[\cosh ({\rm{\Delta }}\mathrm{ln}{X}_{\mathrm{crit}})\right].\end{eqnarray} \tag{ B8 }$$

Numerically solving this equation results in ${\rm{\Delta }}\mathrm{ln}{X}_{\mathrm{crit}}\mathop{\propto }\limits_{\sim }\tfrac{{\rm{\Delta }}T}{\bar{T}}$ as long as $\tfrac{{\rm{\Delta }}T}{\bar{T}}\lesssim 0.7$ (with ${\rm{\Delta }}\mathrm{ln}{X}_{\mathrm{crit}}\to \infty$ as $\tfrac{{\rm{\Delta }}T}{\bar{T}}\to 1$). To obtain a quantitative handle on the scale of ${\rm{\Delta }}\mathrm{ln}{X}_{\mathrm{crit}}$, consider a conservative "extreme" temperature difference of ${T}_{{\rm{E}}}-{T}_{{\rm{M}}}\,=1000\,{\rm{K}}$ on an ultra-hot Jupiter with $\bar{T}=2500\,$ K. This case has $\tfrac{{\rm{\Delta }}T}{\bar{T}}=500/2500=0.2$ and a numerical solution of ${\rm{\Delta }}\mathrm{ln}{X}_{\mathrm{crit}}\,\approx 0.41$ (${\rm{\Delta }}{\mathrm{log}}_{10}{X}_{\mathrm{crit}}\approx 0.18$) or, equivalently, a difference of $\approx 50 \%$. Therefore even minor differences in mixing ratios between the terminators enter the cooling bias regime. As a rule of thumb, compositional differences satisfying ${\rm{\Delta }}{\mathrm{log}}_{10}X\gt 0.3$ (a factor of 2 difference) result in T_1D biased to hundreds of degrees colder than $\bar{T}$.