Clouds and Clarity: Revisiting Atmospheric Feature Trends in Neptune-size Exoplanets

Detailed characterization has long been the goal of exoplanet astrophysics, but only in recent years has this truly been possible. In most cases, exoplanet detection studies have limited ability to probe conditions on these planets directly, leaving exoplanet characterization as a more targeted and expensive second-echelon effort. Characterization efforts can take different forms, but among the most successful has been the atmospheric characterization of gaseous exoplanets, through transmission, emission, or direct spectroscopy. Currently, direct spectroscopy requires sufficiently high-contrast imaging of widely separated, self-luminous exoplanets, while transmission and emission spectroscopy require these planets to transit their host stars, but also require sufficiently warm and large planets to produce detectable atmospheric absorption or planet–star contrast. Given current observational capabilities, transmission spectroscopy has been one of the standout successes in atmospheric characterization, driven primarily by ground-based high-resolution spectrographs as well as space-based observatories like the workhorse Hubble Space Telescope (HST). With over 5500 exoplanets known, of which ∼4100 transit (NASA Exoplanet Science Institute 2020), we now have hundreds of planets characterized through transmission spectroscopy or spectrophotometry.

However, these successes are not without their limits. These instruments often have limited spectral ranges, necessitating the use of multiple observatories, filters, and observations to observe multiple transits of an exoplanet in order to assemble broad-band atmospheric transmission spectra, and many well-known transit spectra are markedly flat (e.g., GJ 1214b; Kreidberg et al. 2014), either due to opaque cloud layers, atmospheric hazes, or the muting of absorption features by high atmospheric metallicity, limiting our ability to determine the content of these atmospheres. A major exception has been the use of HST/Wide Field Camera 3 (WFC3) to observe water vapor absorption in dozens of exoplanetary atmospheres. H₂O has a well-known absorption feature near 1.4 μみゅーm, right in the middle of HST/WFC3's G141 grism. In clear atmospheres, this feature can be quite prominent, allowing for the direct comparison of the cloudiness or clarity of observed exoplanet transmission spectra. Where high-altitude clouds or hazes are present, observed spectra show either no water absorption or muted features, while clear atmospheres show strong water absorption.

Previous work (Stevenson 2016; Crossfield & Kreidberg 2017; Fu et al. 2017; Gao et al. 2020; Yu et al. 2021; Dymont et al. 2022; Estrela et al. 2022; Edwards et al. 2023) has attempted to find relationships between physical exoplanetary parameters (e.g., planet equilibrium temperature, T_eq or surface gravity g) and the strength of observed atmospheric features in transmission. Stevenson (2016) performed the first systematic study of this and introduced the method of examining clarity through the WFC3 water feature amplitude A_H (although focusing on a surface gravity-selected sample), while initial results from Crossfield & Kreidberg (2017) found a linear trend in A_H for a small sample of six warm Neptunes (R_p ∈ [2, 6] R_⊕, T_eq ∈ [500, 1100] K), and showed that hotter planets tend to have less aerosol obscuration than cooler planets. Fu et al. (2017) found a similar positive linear trend for a larger sample, which also included transmission spectra from hotter, more massive planets, and proposed a physical explanation based on cloud condensation instead of the haze formation hypothesis from Crossfield & Kreidberg (2017). Gao et al. (2020) focused on hot (800–2600 K) giant exoplanets, measuring the water feature amplitude compared to model predictions, and found that all but the hottest planets in their sample must be attenuated by clouds. Higher-order polynomial trends have been a popular model for this behavior. Yu et al. (2021) fit a parabolic trend to a sample of nine Neptune-sized planets, but were primarily conducting an analysis of haze evolution and did not elaborate on the implications of this trend. Dymont et al. (2022) examined a larger sample of 25 planets between 200 and 1000 K, but included rocky planets as small as 1.13 R_⊕ (GJ 1132b), and Jovian planets as large as 12.9 R_⊕ (WASP-67 b), and also found hints of the quadratic trend from Yu et al. (2021; without fitting any models) between T_eq and A_H , which depends entirely on K2-18 b and LHS 1140 b. Taking a Neptune-sized 16-planet subset of their full 70-planet analysis, Edwards et al. (2023) find a parabolic trend matching theoretical predictions from Kawashima et al. (2019), while Estrela et al. (2022) fit second- and fourth-order polynomials to 62 planets from 500 to 2500 K, finding two distinct populations: aerosol-free and partially cloudy/hazy planets, with minimal aerosol presence near 1500 K.

Here, we present a trend analysis of a homogeneous sample of Neptune-like exoplanets, in order to show that sample selection is important in order to identify atmospheric feature trends, especially as JWST operations continue and provide higher-quality, broader-wavelength characterization of planetary atmospheres than possible from HST alone.

2.1. Atmospheric Motivation

Typically, observed features are characterized by atmospheric scale heights, converting from relative units (planet-to-star radius ratio or transit depth) to physical units (kilometers), where the atmospheric scale height of a planet is (Deeg & Belmonte 2018)

$$\begin{eqnarray}&&H=\displaystyle \frac{{k}_{{\rm{B}}}{T}_{\mathrm{eq}}}{\mu g},\end{eqnarray} \tag{ 1 }$$

and the depth of an absorption feature at a particular wavelength, measured in numbers of scale heights, is

$$\begin{eqnarray}&&{\delta }_{\lambda }=\displaystyle \frac{{({R}_{{\rm{p}}}+{nH})}^{2}}{{\text{}}{R}_{{\rm{s}}}^{2}}-\displaystyle \frac{{\text{}}{R}_{{\rm{p}}}^{2}}{{\text{}}{R}_{{\rm{s}}}^{2}}\approx \displaystyle \frac{2{{nR}}_{{\rm{p}}}H}{{\text{}}{R}_{{\rm{s}}}^{2}}.\end{eqnarray} \tag{ 2 }$$

Here, K_B is Boltzmann's constant, T_eq the equilibrium temperature of the planet, μみゅー the atmospheric mean molecular weight, g the planet's surface gravity, R_p the planet radius, R_s the stellar radius, and n the number of scale heights.

However, this implies that the strength of an observed absorption feature may actually differ significantly depending on the assumed atmospheric metallicity. For example, two identical planets, at 1 × solar metallicity (μみゅー ≈ 2.3 amu) and 100 × solar metallicity (μみゅー ≈ 3.05 amu) will have H₁₀₀/H₁ ≈ 0.75, and so identical-depth absorption features in a 100 × solar metallicity atmosphere will be "stronger," i.e., correspond to more scale heights than a solar metallicity atmosphere (n₁₀₀ ≈ 1.33n₁).

Since atmospheric metallicity appears to vary across the exoplanetary mass range (Welbanks et al. 2019) and is expected from the solar system and from planetary formation models (Fortney et al. 2013), we must look at physically similar planets in a relatively narrow mass range to mitigate any mass–metallicity effect in order to have any hope of comparing their observed atmospheric features and finding underlying trends. We chose 15 planets for our analysis, including those from Crossfield & Kreidberg (2017) and others that have subsequently been observed with HST WFC3. We started with a specific radius range (R_pl ∈ [2, 6] R_⊕), in order to ensure we were not including rocky super-Earths, and included a few larger planets (HAT-P-26 b, due to its Neptune-like mass, and Kepler 51 b and d as, despite their low masses, they have well-known extended atmospheres and relatively low equilibrium temperatures). We also restricted our sample to temperatures <1000 K, in order to make sure we probe thermally consistent atmospheric physics. As such, we do not include planets such as 55 Cnc e, as it is nearly 1000 K hotter than any other planet in our sample, likely rocky, with no extended atmosphere, and has been shown not to have water vapor in its atmosphere. These planets are shown in Figure 1 and their physical parameters are shown in Table 1.

Table 1. Planets in Our Sample and Their Relevant Parameters

Planet (Spectrum References)	Rpl	Mpl	gpl a	Porb	R*	M*	Teff	Teq a	Ha	AH a	Parameter References
	(R⊕)	(M⊕)	(cm s−2)	(days)	(R⊙)	(M⊙)	(K)	(K)	(km)	(H)
K2-18 b (Benneke et al. 2019a)	2.61	8.63	1241	32.9400	0.44	0.50	3457	253	56	${3.23}_{-0.88}^{+0.98}$	Benneke et al. (2019a)
Kepler-51 d (Libby-Roberts et al. 2020)	9.46	5.7	62	130.1845	0.88	0.98	5670	332	1450	${0.13}_{-0.08}^{+0.24}$	Libby-Roberts et al. (2020)
TOI-270 d (Mikal-Evans et al. 2023)	2.13	4.78	1029	11.3796	0.38	0.39	3506	354	94	${3.50}_{-0.91}^{+1.06}$	Van Eylen et al. (2021)
Kepler-51 b (Libby-Roberts et al. 2020)	6.89	3.69	76	45.1542	0.88	0.98	5670	472	1691	${0.34}_{-0.25}^{+0.41}$	Libby-Roberts et al. (2020)
HD 3167c (Mikal-Evans et al. 2021)	3.01	9.8	1059	29.8454	0.86	0.86	5261	508	131	${1.76}_{-0.89}^{+0.90}$	Christiansen et al. (2017)
GJ 1214 b (Kreidberg et al. 2014)	2.74	8.17	1064	1.5804	0.21	0.18	3250	536	137	${0.09}_{-0.03}^{+0.07}$	Cloutier et al. (2021)
GJ 3470 b (Benneke et al. 2019b)	3.88	13.73	893	3.3366	0.48	0.51	3652	597	182	${1.03}_{-0.26}^{+0.20}$	Biddle et al. (2014)
GJ 436 b (Knutson et al. 2014)	3.96	21.7 b	1355	2.6439	0.44	0.45	3585	619	125	${0.58}_{-0.39}^{+0.32}$	Knutson et al. (2011)
GJ 9827 d (Roy et al. 2023)	2.07	5.2	1189	6.2015	0.61	0.61	4269	621	142	${3.21}_{-1.05}^{+0.99}$	Rodriguez et al. (2018)
TOI-674 b (Brande et al. 2022)	5.25	23.6	838	1.9771	0.42	0.42	3514	661	215	${1.05}_{-0.41}^{+0.46}$	Murgas et al. (2021)
HD 97658 b (Guo et al. 2020)	2.12	8.3	1809	9.4897	0.73	0.85	5212	681	103	${0.06}_{-0.05}^{+0.74}$	Ellis et al. (2021)
HAT-P-11 b (Fraine et al. 2014)	4.73	25.74	1133	4.8878	0.75	0.81	4780	796	193	${2.69}_{-0.64}^{+0.59}$	Bakos et al. (2010)
HD 106315c (Kreidberg et al. 2022)	4.35	15.2	787	21.0570	1.30	1.09	6327	811	281	${1.91}_{-0.87}^{+1.14}$	Barros et al. (2017)
HIP 41378 b (This work)	2.9	8.29 c	965	15.5712	1.40	1.15	6199	904	255	${2.64}_{-2.50}^{+2.62}$	Vanderburg et al. (2016)
HAT-P-26 b (Wakeford et al. 2017)	6.3	18.75	458	4.2345	0.79	0.82	5079	909	541	${2.97}_{-0.45}^{+0.40}$	Hartman et al. (2011)

Notes. Scale heights were calculated assuming μみゅー = 3.05 amu, and equilibrium temperatures calculated assuming an albedo of 0.3.

^a Value calculated as part of this analysis. ^b Knutson et al. (2014). ^c Calculated using the Chen & Kipping (2017) mass–radius relationship.

Download table as:  ASCIITypeset image

3.1. Planetary Parameters and Scale Heights

From the literature values for planetary and system parameters (R_p, M_p, P_orb, M_*), we recalculated the orbital semimajor axis a (according to Equation (3)), planet surface gravity g (according to Equation (4)), as well as planet equilibrium temperature T_eq (according to Equation (5)) for all planets in our sample, assuming a planetary albedo of A = 0.3.

$$\begin{eqnarray}&&a={\left({P}_{\mathrm{orb}}^{2}\displaystyle \frac{{{GM}}_{* }}{4{\pi }^{2}}\right)}^{1/3}\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&g=\displaystyle \frac{{{GM}}_{{\rm{p}}}}{{\text{}}{R}_{{\rm{p}}}^{2}}\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&{T}_{\mathrm{eq}}={{\rm{T}}}_{\mathrm{eff}}\sqrt{\displaystyle \frac{{\text{}}{R}_{* }}{2a}}{(1-A)}^{1/4}.\end{eqnarray} \tag{ 5 }$$

To calculate the scale heights of the atmospheres of these planets, we assumed μみゅー = 3.05 amu (as these Neptune-sized worlds are more likely to be nearer 100 × solar metallicity than 1 × solar; Welbanks et al. 2019), and calculated H according to Equation (1). Our calculated values are shown in Table 1.

3.2. Atmospheric Model Fitting and Retrievals

We used the open-source code petitRADTRANS (Mollière et al. 2019) to conduct isothermal, free-chemistry atmospheric retrievals of the archival HST/WFC3 spectra of the planets in our sample. Our retrieval framework fits for five parameters: planet radius, planet mass, isothermal atmospheric temperature, H₂O abundance (${\mathrm{log}}_{10}X$, where X is the H₂O mass fraction), and gray cloud deck pressure (${\mathrm{log}}_{10}P$, in bars). We placed Gaussian priors on planetary mass and radius based on the literature values and our calculated gravity distributions, a uniform prior on isothermal atmospheric temperature of T_eq ± 200 K, and placed uniform priors on the log₁₀ cloud deck pressure (in bars, ${ \mathcal U }(-6,2)$) and log₁₀ water abundance (as mass fractions, ${ \mathcal U }(-6,0)$). Our retrievals also included Rayleigh scattering from H₂ and He, as well as collisionally induced absorption of H₂–H₂ and H₂–He as continuum opacities. We fixed the reference pressure of our retrievals at 0.01 bar, and calculated model spectra at R = 1000.

For each planet, to measure A_H , we sampled 100 random model spectra from our retrievals and rebinned them to a resolution of R = 250, using the flux-conserving resampling method in astropy's (Astropy Collaboration et al. 2013, 2018, 2022) specutils package. For each sampled spectrum i, we calculated ${{\rm{A}}}_{{{\rm{H}}}_{{\rm{i}}}}={{\rm{n}}}_{{\rm{i}}}=({\delta }_{1.4,{\rm{i}}}-{\delta }_{1.25,{\rm{i}}})({\text{}}{{\rm{R}}}_{* }^{2}/2{{\rm{H}}}_{{\rm{i}}}{\text{}}{{\rm{R}}}_{{\rm{p}},i})$, as rearranged from Equation (2), where δでるた_{λらむだ,i} is the transit depth of the sampled spectrum i at wavelength λらむだ, R_* the stellar radius for the system in question, R_p,i the retrieved radius of the planet for sampled spectrum i, H_i the atmospheric scale height (calculated with the retrieved mass and radius values for each sampled spectrum), and ${n}_{i}={A}_{{H}_{i}}$ the amplitude of the feature in numbers of scale heights H_i . The final A_H value and its confidence interval for a particular planet was calculated from the 16th, 50th, and 84th percentiles of the sampled values for n_i .

Zoom In Zoom Out Reset image size

The observed transmission spectra and atmospheric models for our planets are shown in Figure 2.

Zoom In Zoom Out Reset image size

3.3. Bayesian Trend Modeling

We searched for trends relating A_H to various planetary parameters in our sample, such as equilibrium temperature, planet surface gravity, and stellar metallicity, but found no hints of useful relationships between any other than A_H and equilibrium temperature. We then worked to identify any possible temperature dependent trend in a hierarchical Bayesian framework, in order to both distinguish between possible polynomial models as well as quantify the inherent scatter in the data. We also exclude two planets in our sample, Kepler-51 b and d, from the trend fitting, as these so-called "super-puffs" are likely not in the same physical category as the rest of the sample due to their anomalously low densities and surface gravities. It is unclear whether these, along with GJ 1214 b and its exceptionally flat spectrum, are part of the same population as the other cold planets that have strong absorption features. Perhaps below a certain T_eq, the clarity of sub-Neptune atmospheres is essentially random with respect to T_eq, or perhaps something splits these planets into clear and cloudy populations. Studies of GJ 1214 b's atmosphere with JWST continue to indicate it too may be unlike its peers, with especially high (>1000 × solar) metallicity (Gao et al. 2023). As such, we also exclude it from the final trend fitting.

We aim to distinguish between several trend models: one with zero correlation between A_H and temperature (A_H = some constant c), a linear model (A_H = bt + c) and one with a parabolic trend (A_H = at² + bt + c), where t is the equilibrium temperature. We also included an extra parameter σしぐま_i to account for unknown scatter (whether a result of stochastic astrophysical processes or various modeling insufficiencies in our polynomial trends) within our sample, added in quadrature to the measured A_H uncertainties as ${\sigma }_{\mathrm{fit}}\,=\sqrt{{\sigma }_{{A}_{H}}^{2}+{\sigma }_{i}^{2}}$.

We used broad uniform priors on our model parameters, and used the No U-Turn Sampler available in exoplanet (Foreman-Mackey et al. 2021) for model fitting and posterior sampling, running four chains for 3000 tuning steps and 2000 draws each. After all sampler runs converged (as evidenced by the Gelman-Rubin statistic $\hat{R}$ reaching a value of 1), we calculated ${\chi }_{\nu }^{2}$ and Bayesian Information Criterion (BIC) values for each model and found that the parabolic model

$$\begin{eqnarray*}&&{A}_{H}=2.3\times {10}^{-5}{t}^{2}-0.028t+9.5\end{eqnarray*}$$

was preferred to either the constant or linear models, with ${\chi }_{\nu }^{2}=2.0$ and BIC=25.1. The best-fit parabolic model also includes 0.56 scale height of added scatter. Our priors, fitted parameters, and model comparisons are shown in Table 2. Including the Kepler-51 planets and GJ 1214 b results in a worse fit to the data, with ${\chi }_{\nu }^{2}=23$ and BIC = 287, as well as increases the added scatter to 0.95 scale height.

Table 2. Bayesian Model Fitting Priors, Fitted Parameters, and Model Comparisons

Model	Parameter	Prior	Fitted Value
Constant, ${\chi }_{\nu }^{2}=5.9$, BIC = 67.9
	c	${ \mathcal U }(-50,50)$	${1.8}_{-0.27}^{+0.29}$
	σしぐまi	${ \mathcal U }(0,10)$	${1.1}_{-0.20}^{+0.24}$
Linear, ${\chi }_{\nu }^{2}=6.8$, BIC = 73.5
	b	${ \mathcal U }(-50,50)$	$-{7.1}_{-15}^{+15}\times {10}^{-3}$
	c	${ \mathcal U }(-50,50)$	${2.3}_{-1.0}^{+1.0}$
	σしぐまi	${ \mathcal U }(0,10)$	${1.1}_{-0.20}^{+0.25}$
Parabola, ${\chi }_{\nu }^{2}=2.0$, BIC = 25.1
	a	${ \mathcal U }(-50,50)$	${2.3}_{-0.47}^{+0.45}\times {10}^{-5}$
	b	${ \mathcal U }(-50,50)$	$-{0.028}_{-0.0055}^{+0.0058}$
	c	${ \mathcal U }(-50,50)$	${9.5}_{-1.8}^{+1.7}$
	σしぐまi	${ \mathcal U }(0,10)$	${0.56}_{-0.15}^{+0.20}$

Download table as:  ASCIITypeset image

After identifying the best-fit parabolic trend, we compared our findings to a grid of sub-Neptune atmospheric models from Morley et al. (2015). This grid modeled the atmosphere of GJ 1214 b in clear, cloudy, and hazy configurations, across various metallicities (from 50 × to 1000 × solar) and cloud/haze parameters. The grid was also calculated at varying insolations (compared to GJ 1214 b's actual insolation), corresponding to ∼450 K at the low end to 1100 K on the high end. As with our retrievals, we resampled the model spectra to R = 250, and calculated A_H using the same GJ 1214 b planetary parameters (but recalculating H according to the model T_eq). The equilibrium cloud models include salt and sulfide species (ZnS, KCl, and Na₂S) calculated at varying f_sed and metallicity, and are more fully described in Section 2.3 of Morley et al. (2015). The haze models were all calculated at 50 × solar metallicity, included soots with five precursor species (C₂H₂, C₂H₄, C₂H₆, C₄H₂, and HCN), calculated at varying particle sizes and precursor conversion fractions, and are more fully described in Section 2.4 of Morley et al. (2015). Using the 100 × and 300 × solar cloudy models and the 1% and 10% haze precursor conversion fractions to bracket reasonable parameter ranges for our sample, we clearly see in Figure 3 that almost all of our sample planet spectra strongly attenuated by atmospheric aerosols compared to the clear atmosphere models.

Zoom In Zoom Out Reset image size

We calculate ${\chi }_{\nu }^{2}$ values for our planets with T_eq > 450 K compared to the various contours for our clear, cloudy, and hazy models to find which set of models is most consistent with our observed A_H values. These model comparisons are shown in Table 3. We find that all our planets are likely to be attenuated by atmospheric aerosols (sometimes very significantly, up to ∼6 scale heights). This attenuation is strongest between ∼500 and ∼700 K, corresponding to a range of especially efficient aerosol production identified by theoretical models (Morley et al. 2015; Gao et al. 2020). The attenuation is also stronger for the 100 × solar models compared to the 300 × solar models.

Table 3. Atmospheric Models and Model Comparisons

Model	${\chi }_{\nu }^{2}$
Clear Models
100 × solar	149
300 × solar	32
Cloudy Models
fsed = 1, 100 × solar	62
fsed = 1, 300 × solar	13
fsed = 0.1, 100 × solar	2.9
fsed = 0.1, 300 × solar	4.7
fsed = 0.01, 100 × solar	3.8
fsed = 0.01, 300 × solar	9.0
Hazy Models
r = 1 μみゅーm, 1% conversion	102
r = 1 μみゅーm, 10% conversion	8.6
r = 0.1 μみゅーm, 1% conversion	8.1
r = 0.1 μみゅーm, 10% conversion	12.6
r = 0.01 μみゅーm, 1% conversion	7.4
r = 0.01 μみゅーm, 10% conversion	12.9

Download table as:  ASCIITypeset image

Even so, we note some interesting points. First, with respect to the cloudy models, the majority of our planets have low sedimentation efficiency, f_sed = 0.1 or lower, and the entire sample is most consistent with f_sed = 0.01 and 100 × solar metallicity. With respect to the hazy models, our sample is most consistent with small haze particles (r = 0.01 μみゅーm) and small (1%) precursor conversion rates. Overall, the cloudy models fit better than the hazy models, although the sparsity of the grid means this is a relatively low-significance finding. More detailed modeling efforts may clarify this, as discussed in Section 6.

Second, while they were not included in the trend fitting, the Kepler-51 planets cannot be totally ignored. While these "super-puffs" clearly stand out as exceptionally low-density planets, other temperate, puffy sub-Neptunes have been discovered since (e.g., TOI-700 c, WASP-193 b), which will also be amenable to atmospheric transmission spectroscopy. Observing as many cold planets as possible is sorely needed to determine whether our observed trend holds. These observations would identify whether cold sub-Neptunes exist in a continuum between the slightly attenuated atmospheres that follow our trend and the especially cloudy and hazy planets like the Kepler-51 planets and GJ 1214 b. If these are outliers among the broader population, it bodes well for future studies of habitable zone systems. Perhaps there is a range of temperatures where cloud rainout is efficient (especially for heavier species like ZnS and KCl), while low insolation means haze production is not, giving especially clear views into these temperate planets' atmospheres.

Finally, we find that, among our sample of Neptune-sized exoplanets, there is an intrinsic scatter in the clarity of these atmospheres of ∼0.5 scale height. This scatter may be evidence that our trend model choices are insufficient, or possibly be due to some randomness in the outcomes of planet formation and evolutionary processes. Perhaps the trend identified is analogous to a "main sequence" of planetary atmospheres, and (taking into account the super-puffs) future observations with JWST will be able to confirm or reject our identification of this scatter.

We also note that, while we used water vapor as the major opacity source in this analysis, methane absorption has significant overlap with water in near-IR bands, and may dominate compared to water vapor in certain exoplanetary atmospheres <600 K (Bézard et al. 2022). In fact, this has been identified in K2-18 b, previously thought to have water vapor in its atmosphere from HST data (Benneke et al. 2019a), but now shown instead to have methane from broad-wavelength JWST observations (Madhusudhan et al. 2023). With this in mind, we also ran methane-only retrievals for K2-18 b and TOI-270 d to compare to our water-only retrievals and found the methane-calculated A_H values to be consistent within the 1σしぐま uncertainties. The substitution of CH₄ for H₂O in these coolest planets does not change our conclusions. These are shown in Figure 3 as open green markers.

These results echo previous analyses: our model comparisons provide evidence for a cloud-based mechanism shaping this trend (predicted by Fu et al. 2017) over a haze-driven trend (predicted by Crossfield & Kreidberg 2017 and Yu et al. 2021). Qualitatively, our quadratic trend behaves similarly to previous second-order trends fit to samples that extend to low temperatures (Yu et al. 2021; Edwards et al. 2023), showing significant increases in atmospheric clarity at colder (<500 K) and hotter (>700 K) temperatures. The major differences between these trends are primarily due to sampling, as the inclusion of particular planets (like Kepler-51 d in Yu et al. 2021 and Edwards et al. 2023) can shift the parabola minimum to colder equilibrium temperatures. Robustness against planetary samples is encouraging, as it appears the observational predictions made by these trends will be useful for a broad range of future spectroscopic studies. Most of our measured A_H values are also consistent to previous analyses within error, although our method differs to previous analyses due to our retrieval framework. This confirms that, as an observational metric, A_H is robust against methodological changes (even if the physical insights gained about any particular planet transmission spectrum from A_H alone appear to be limited).

Despite the proliferation of HST WFC3 transmission spectra in recent years, significant trends relating the clarity of these atmospheres to planetary and system parameters still evade us. WFC3, while especially sensitive to water vapor and uniform cloud absorption within a narrow spectral range, cannot characterize the full panoply of absorbers present in exoplanet atmospheres (perhaps most significantly shown by K2-18 b's JWST methane detection in Madhusudhan et al. 2023 over the HST water detection from Benneke et al. 2019a), nor precisely constrain physical cloud models. As already shown from the first year of JWST analyses, exoplanetary atmospheres are complex, dynamic, and contain clear evidence of atmospheric physics beyond anything HST could provide (Ahrer et al. 2023; Alderson et al. 2023; Feinstein et al. 2023; JWST Transiting Exoplanet Community Early Release Science Team et al. 2023; Rustamkulov et al. 2023; Tsai et al. 2023).

Among the scheduled and completed JWST programs (through Cycle 2), there are 106 unique exoplanet systems. Of the planets in these systems, 36 fit our sample radius criteria (R_pl ∈ [2, 6] R_⊕). In total, 8 of the planets we analyzed in this work are on this list as well, although these alone are not sufficient to adequately sample the range of equilibrium temperatures. We need to add more exoplanets to the JWST observational sample, especially temperate planets that are large enough to sustain extended atmospheres, so we can confirm the <500 K leg of the trend observed in this work.

We also need significantly improved model grids for sub-Neptune atmospheres. Aerosol production is a complex process, and low-resolution observations have historically been unable to reveal much about these. New JWST observations, with higher resolution and a much broader wavelength range, have already identified signatures of wavelength-dependent silicate cloud absorption in the planetary-mass object VHS 1256 b (Miles et al. 2023). Haze precursors like SO₂ have been identified in the Hot Jupiter WASP-39 b (Alderson et al. 2023; Rustamkulov et al. 2023; Tsai et al. 2023), and GJ 1214 b's high-altitude aerosols have been found to be very reflective (A = 0.51; Kempton et al. 2023). Extending modeling efforts to colder planets, across a wide range of metallicities, cloud species, particle sizes, atmospheric mixing, and other parameters, should enable us to constrain these processes in these planets' atmospheres.

We have identified a coherent trend relating the clarity of Neptune-sized planetary atmospheres to their equilibrium temperatures, as well as modeling the intrinsic scatter in planetary clarity measurements, across an ∼800 K temperature range. These findings show behavior similar to previous analyses, again finding a parabolic trend with clearer atmospheres at cooler (<500 K) and hotter (>700 K) temperatures, with a minimum between 500–700 K. We reproduce this behavior for our restricted sample of Neptune-like planets compared to more diverse samples in previous analyses, implying the observational implications of this work should be applicable to many future HST and JWST targets. Compared to models, all observed planetary atmospheres are consistent with attenuation by at least some cloud cover as a result of vertically extended clouds or hazes with small particles (f_sed = 0.1, r = 0.01 μみゅーm), and have relatively high metallicity (∼100 × solar). Confirmation or rejection of cold-temperature theoretical predictions still eludes us, as we need updated self-consistent modeling efforts across a range of atmospheric parameters for these cold planets. JWST observations in the next few cycles should both continue to test our identified trend and help identify the atmospheric diversity that may be present among cold (T_eq < 500 K) planets.

We thank Luis Welbanks for helpful discussions and the anonymous referee whose comments significantly improved the work.

Facility: Exoplanet Archive -

Software: Astropy (Astropy Collaboration et al. 2013, 2018, 2022), exoplanet (Foreman-Mackey et al. 2021).

HIP 41378 b was included as a target in HST program GO-15333, which also observed HD 3167c, GJ 9827 d, TOI-674 b, and HD 106315 c.

We observed two transits of HIP 41378 b with HST's WFC3 instrument, using the G141 grism, before dropping this planet in favor of TOI-674 b (Brande et al. 2022) due to scheduling difficulties. We followed standard data reduction procedures for HST/WFC3 time series data, as documented in Kreidberg et al. (2022). Following common practice, we fit the light curves simultaneously with a transit model and an instrument systematics model. Our transit model included variable R_p/R_*, and other system parameters were fixed to the best-fit values from Berardo et al. (2019). The limb-darkening parameters were fixed on predictions from a PHOENIX stellar model at the effective temperature of the host star. The instrument systematics model consisted of a visit-long linear slope, a scaling factor between upstream and downstream spatial scans, and exponential ramps fit to each orbit (following the model ramp parameterization from Kreidberg et al. 2022). The observed transmission spectrum is given in Table 4. The raw data for HIP 41378 b is available at the Mikulski Archive for Space Telescopes at the Space Telescope Science Institute at the following doi:10.17909/ryg6-3c46.

Table 4. HIP 41378 b HST/WFC3 G141 Transmission Spectrum

Wavelength	Transit Depth ${({R}_{{\rm{p}}}/{R}_{* })}^{2}$
(μみゅーm)	(ppm)
1.175	347 ± 32.7
1.225	273 ± 32.2
1.275	244 ± 35.5
1.325	326 ± 28.6
1.375	302 ± 32.7
1.425	285 ± 33.2
1.475	320 ± 34.9
1.525	296 ± 34.7
1.575	283 ± 35.6
1.625	294 ± 41.1

Download table as:  ASCIITypeset image