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

Over the last decade, precise exoplanet transmission spectroscopy has revealed the atmospheres of dozens of exoplanets, driven largely by observatories like the Hubble Space Telescope. One major discovery has been the ubiquity of atmospheric aerosols, often blocking access to exoplanet chemical inventories. Tentative trends have been identified, showing that the clarity of planetary atmospheres may depend on equilibrium temperature. Previous work has often grouped dissimilar planets together in order to increase the statistical power of any trends, but it remains unclear from observed transmission spectra whether these planets exhibit the same atmospheric physics and chemistry. We present a reanalysis of a smaller, more physically similar sample of 15 exo-Neptune transmission spectra across a wide range of temperatures (200–1000 K). Using condensation cloud and hydrocarbon haze models, we find that the exo-Neptune population is best described by low cloud sedimentation efficiency (f sed ∼ 0.1) and high metallicity (100 × solar). There is an intrinsic scatter of ∼0.5 scale height, perhaps evidence of stochasticity in these planets’ formation processes. Observers should expect significant attenuation in transmission spectra of Neptune-size exoplanets, up to 6 scale heights for equilibrium temperatures between 500 and 800 K. With JWST's greater wavelength sensitivity, colder (<500 K) planets should be high-priority targets given their clearer atmospheres, and the need to distinguish between the “super-puffs” and more typical gas-dominated planets.

1. INTRODUCTION arXiv:2310.07714v2[astro-ph.EP] 18 Jan 2024 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 1214 b, 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 /WFC3 to observe water vapor absorption in dozens of exoplanetary atmospheres.H 2 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 highaltitude 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;Edwards et al. 2023;Estrela et al. 2022) has attempted to find relationships between physical exoplanetary parameters (e.g.planet equilib-rium 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 6 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 K -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 9 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 K and 1000 K, but included rocky planets as small as 1.13 R ⊕ (GJ 1132 b), 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 fourthorder 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.

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): and the depth of an absorption feature at a particular wavelength, measured in numbers of scale heights, is: 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 100 /H 1 ≈ 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 100 ≈ 1.33n 1 ).
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 which 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 it's 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 like 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 and their physical parameters are shown in Table 1. are not yet scheduled or approved for JWST observations.Selected targets have been labeled.

Planetary Parameters and Scale Heights
From literature values for planetary and system parameters (R p , M p , P orb , M * ), we recalculated the orbital semimajor axis a (according to Eq. 3), planet surface gravity g (according to Eq. 4), as well as planet equilibrium temperature T eq (according to Eq. 5) for all planets in our sample, assuming a planetary albedo of A = 0.3.
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 Eq. 1.Our calculated values are shown in Table 1.

Atmospheric Model Fitting and Retrievals
We used the open-source code petitRADTRANS (Mollière et al. 2019) to conduct isothermal, freechemistry atmospheric retrievals of the archival HST /WFC3 spectra of the planets in our sample.Our  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 fluxconserving resampling method in astropy's (Astropy Collaboration et al. 2013, 2018, 2022) specutils package.For each sampled spectrum i, we calculated ), as re-arranged 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 Hi 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 16 th , 50 th , and 84 th percentiles of the sampled values for n i .
The observed transmission spectra and atmospheric models for our planets are shown in Figure 2.

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's 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 2 + 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 σ fit = σ 2 A H + σ 2 i .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 4 chains for 3000 tuning steps and 2000 draws each.After all sampler runs converged (as evidenced by the Gelman-Rubin statistic R reaching a value of 1), we calculated χ 2 ν and BIC values for each model and found that the parabolic model A H = 2.3 × 10 −5 t 2 − 0.028t + 9.5 was preferred to either the constant or linear models, with χ 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 χ 2 ν = 23 and BIC=287, as well as increases the added scatter to 0.95 scale height.

MODEL COMPARISONS
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 2 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 2 H 2 , C 2 H 4 , C 2 H 6 , C 4 H 2 , 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.

RESULTS
We calculate χ 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.
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.
continuum between the slightly attenuated atmospheres which 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 Neptunesized 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 4 for H 2 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 secondorder trends fit to samples that extend to low temperatures (Yu et al. 2021;Edwards et al. 2023), showing significant increases in atmospheric clarity at cold (< 500K) and hotter (> 700K).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 ob-servational 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).

COMMUNITY PRIORITIES
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 (JWST Transiting Exoplanet Community Early Release Science Team et al. 2023;Rustamkulov et al. 2023;Alderson et al. 2023;Feinstein et al. 2023;Ahrer 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 ⊕ ). 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 < 500K 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 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 2 have been identified in the Hot Jupiter WASP-39 b (Rustamkulov et al. 2023;Alderson 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.

CONCLUSIONS
We have identified a coherent trend relating the clarity of Neptune-sized planetary atmospheres to their equilibrium temperatures, as well as modeling the intrin-sic scatter in planetary clarity measurements, across an ∼ 800K temperature range.These findings show behavior similar to previous analyses, again finding a parabolic trend with clearer atmospheres at cooler (< 500K) and hotter (> 700K), 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 coldtemperature 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.

Figure 1 .
Figure 1.This work's sample overlaid with known transiting planets.Hexagons indicate planets in our sample also being observed by JWST through Cycle 2, while triangles are not yet scheduled or approved for JWST observations.Selected targets have been labeled.

Figure 2 .
Figure 2. All HST /WFC3 G141 transmission spectra of our exo-Neptunes, scaled to show the relative AH values.The retrieved models are colored by the planetary equilibrium temperature, and the spectra are ordered by AH .

Figure 3 .
Figure 3. Retrieved spectral feature amplitude, AH , compared to clear, cloudy, and hazy atmosphere models from Morley et al. (2015).For the clear and cloudy models, the dotted lines indicate a 100× Solar metallicity atmosphere, and the dashed lines a 300× Solar atmosphere.All hazy models were calculated at 50× Solar metallicity, the dotted lines are 1% haze precursor conversion, and the dashed lines are 10%.The green open markers, shown for context, are the measured CH4-retrieved AH values for K2-18 b and TOI-270 d.The black open markers, shown for context, are the Kepler-51 planets and GJ 1214 b, which are excluded from the trend analysis.

Table 1 .
Planets in our sample and their relevant parameters.Scale heights were calculated assuming µ = 3.05 amu, and equilibrium temperatures calculated assuming an albedo of 0.3.* Value calculated as part of this analysis.† Knutson et al. (2014) ‡ Calculated using the Chen & Kipping (2017) mass-radius relationship retrieval framework fits for five parameters: planet mass, planet radius, isothermal atmospheric temperature, H 2 O abundance (log 10 X, where X is the H 2 O mass fraction), and gray cloud deck pressure (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 10 cloud deck pressure (in bars, U(−6, 2)) and log 10 water abundance (as mass fractions, U(−6, 0)) Our retrievals also included Rayleigh scattering from H 2 and He, as well as collisionally-inducedabsorption of H 2 -H 2 and H 2 -He as continuum opacities.

Table 2 .
Bayesian model fitting priors, fitted parameters, and model comparisons.