A unique hot Jupiter spectral sequence with evidence for compositional diversity

The emergent spectra of close-in, giant exoplanets (‘hot Jupiters’) are expected to be distinct from those of self-luminous objects with similar effective temperatures because hot Jupiters are primarily heated from above by their host stars rather than internally from the release of energy from their formation1. Theoretical models predict a continuum of dayside spectra for hot Jupiters as a function of irradiation level, with the coolest planets having absorption features in their spectra, intermediate-temperature planets having emission features due to thermal inversions and the hottest planets having blackbody-like spectra due to molecular dissociation and continuum opacity from the H− ion2–4. Absorption and emission features have been detected in the spectra of a number of individual hot Jupiters5,6, and population-level trends have been observed in photometric measurements7–15. However, there has been no unified, population-level study of the thermal emission spectra of hot Jupiters as there has been for cooler brown dwarfs16 and transmission spectra of hot Jupiters17. Here we show that hot Jupiter secondary eclipse spectra centred around a water absorption band at 1.4 μm follow a common trend in water feature strength with temperature. The observed trend is broadly consistent with model predictions for how the thermal structures of solar-composition planets vary with irradiation level, but is inconsistent with the predictions of self-consistent one-dimensional models for internally heated objects. This is particularly the case because models of internally heated objects show absorption features at temperatures above 2,000 K, whereas the observed hot Jupiters show emission features and featureless spectra. Nevertheless, the ensemble of planets exhibits some degree of scatter around the mean trend for solar-composition planets. The spread can be accounted for if the planets have modest variations in metallicity and/or elemental abundance ratios, which is expected from planet formation models18–21. A population study of near-infrared spectra of 19 hot giant planets shows a correlation between the strength of the 1.4 μm water band and temperature, which is broadly regulated by irradiation. However, the observed scatter around the mean is indicative of the effect of individual planetary formation pathways on the composition.

We performed a statistical analysis of 19 hot Jupiter secondary eclipse spectra obtained with the Wide Field Camera 3 (WFC3) instrument on the Hubble Space Telescope (HST) using the G141 grism between 1.1 and 1.7 µm.This bandpass is primarily sensitive to water vapor in exoplanet atmospheres, and the largest molecular feature in this wavelength range is a water vapor absorption band centered at about 1.4 µm.Over the last decade a large sample of exoplanets have been observed using WFC3+G141 to understand their atmospheric water abundances 17,22 , and it has become an important tool in understanding exoplanet atmospheres.
We analyzed six new data sets following the data reduction procedure outlined in the Methods.We also performed a reanalysis of the spectrum of one planet, Kepler-13Ab.We combined these seven new analyses with twelve results from the literature to form a complete sample of planets observed with HST/WFC3+G141 in this wavelength region.Supplementary Table 1 contains detailed information on each of the twelve arXiv:2110.11272v1[astro-ph.EP] 21 Oct 2021 literature results we considered.The planets in this study have observed dayside temperatures in the HST/WFC3+G141 bandpass between 1450 − 3500 K and radii between 0.9 − 2.0 Jupiter radii.The full set of 19 spectra are shown in Figure 1.
Manjavacas et al. (2019) 16 presented the spectra of 10 of the hot Jupiters in this study but did not examine in detail the physical causes for the observed spectral features.Melville et  al. (2020)  23 similarly examined the spectra of a different subset of 12 hot Jupiters but only analyzed them in the context of models with fixed temperature-pressure (T-P) profiles, with no feedback between the T-P profile and the chemistry.Here we expand on these studies by doubling the sample of hot Jupiter secondary eclipse spectra and comparing the spectra to a grid of models with fully consistent T-P profiles to understand in detail what drives their feature strengths.Because our models combine a set of basic self-consistent assumptions which are expected to hold true for hot Jupiters (e.g., energy balance in the atmosphere and thermochemical equilibrium 24 ) with a complete set of relevant opacities, we can use them to create self-consistent predictions for hot Jupiter spectra, which can then be compared to the observed data.
Baxter et al. (2020) 13 presented an analysis of changes in Spitzer Space Telescope emission photometry with temperature and also examined a subset of planets observed with HST.However, because this study focused on broadband photometry, they were only able to give broad constraints on the hot Jupiter population, such as that high C/O ratios are disfavored.This study expands on that work by uniformly analyzing all HST thermal emission spectra and performing a more comprehensive analysis of their compositional diversity.
We created a grid of cloud-free irradiated 1D selfconsistent radiative-convective-thermochemical equilibrium models to compare to the dayside HST/WFC3+G141 thermal emission observations.These models were created using the Sc-CHIMERA framework 4,13,[25][26][27][28] which includes a broad array of opacity sources that are important for the temperature range explored here, including atomic and ionic opacities that are relevant at the high temperatures of ultra-hot Jupiters 24 .A full description of the models and complete list of opacities can be found in the Methods.
Figure 2 shows the T-P profiles and resultant secondary eclipse spectra for our fiducial model, which uses system parameters for a standard hot Jupiter (stellar effective temperature T eff = 5300 K, planetary gravity g = 10 m/s 2 , planetary metallicity M H = 0.0, planetary carbon-to-oxygen abundance ratio C O = 0.55, and planetary internal temperature T int = 150 K).Models at different temperatures were created by scaling the incident stellar flux to match the specified irradiation temperature.Figure 2 also shows the ratio of the absorption mean opacity (κ J ) to the Planck mean opacity (κ B ) as a function of equilibrium temperature at a pressure of 10 −2 bar, which is approximately the photospheric pressure in the HST/WFC3 bandpass (see the Methods for a full description of these opacities).This ratio describes the relative efficiency of stellar absorption vs. thermal re-radiation at that layer in the planet's atmosphere 29 .
In addition to the fiducial model grid, we created models with a variety of atmospheric/system parameters to see how individual parameters impact the 1D vertical structure and resulting population level trends observed in the dayside emission spectra.We examined models with stellar T eff = 3300 K, 4300 K, 6300 K, 7200 K, and 8200 K; planetary gravity, g = 1 m/s 2 and 100 m/s 2 ; metallicity, M H = −1.5 and 1.5; and C O = 0.01 and 0.85.We also included a model where the internal temperature varies with the planetary irradiation temperature to capture the internal entropy change that could be the cause of the hot Jupiter radii inflation 30 .Furthermore, we tested models in which the TiO and VO opacity were removed ad hoc until temperatures above 2000 K, 2500 K, or 3000 K in order to simulate cold-trapping in cooler regions of the atmosphere 3,31,32 (see the Methods section for a full description of these models).For all of these models, only one parameter was varied at a time while the other parameters were held fixed to the values in the fiducial model (e.g., a slice along a given parameter dimension).(c) Ratio of the absorption mean opacity (κ J ) to the Planck mean opacity (κ B ) as a function of equilibrium temperature (T eq ) in these models at a pressure of 10 −2 bar, assuming zero albedo and full redistribution.This ratio describes the relative efficiency of heating vs. cooling in the models 29 , and a ratio of κ J /κ B > 1 generally indicates the presence of a thermal inversion in the T-P profile.This panel is plotted with temperature on the y-axis for ease of comparison to Figure 3.
Our models predict three primary spectral regimes.At the lowest dayside temperatures (T day < 2100 K for the fiducial model), the models exhibit absorption features due to monotonically decreasing temperature profiles.At intermediate temperatures (2100 K < T day < 3000 K for the fiducial model), the modeled thermal structures exhibit a rising temperature with increasing altitude (decreasing pressure) due to the gas-phase onset of TiO and VO which push K J /K B > 1, in turn causing emission features.At the highest temperatures (T day > 3000 K for the fiducial model), the models still show strong thermal inversions (becoming stronger primarily due to the dissociation of water, an efficient coolant) but the resulting secondary eclipse spectra are relatively featureless because of a combination of high-temperature effects such as molecular dissociation and the onset of H − opacity, which cause all the WFC3+G141 wavelengths to probe the same altitude/pressure level, hence brightness temperature 3,4,[24][25][26] .The exact temperatures of the transitions between these regimes, as well as the strength of absorption and emission features present in the models, depend on the parameters of each set of models.
For both the models and the population of 19 observed hot Jupiters, we examined the degree of absorption or emission observed in the water feature at 1.4 µm, the primary feature in the HST/WFC3+G141 bandpass.We quantified their deviation from a blackbody using an HST water feature strength metric, which is illustrated in Supplementary Figure 1.For each data set, we first fit a blackbody to the two "out-ofband" regions of the spectrum, which have wavelengths of 1.22 − 1.33 µm and 1.53 − 1.61 µm and are defined based on where the models in Figure 2 show minimal water opacity.The temperature of this blackbody is referred to throughout this paper as the observed dayside temperature (T day ) in this bandpass.The water feature strength is then defined as where F B,in and F obs,in are the flux of the fitted blackbody and the observed data, respectively, in the "in-band" region shown in Supplementary Figure 1.The "in-band" wavelength region extends from 1.35 − 1.48 µm and captures the center of the primary water band observed in the HST/WFC3 bandpass.The shaded regions in Figure 1 show the extent of the "outof-band" and "in-band" regions.From this definition, S H 2 O will have a positive value when a water feature is observed in absorption, a negative value when a feature is observed in emission, and a value of zero if a blackbody is observed.We note that we use S H 2 O here instead of the traditional infrared Jand H-bands because the J-and H-bands exclude the strongest part of the water band at ≈ 1.4 µm.Therefore, the S H 2 O index we define gives us greater sensitivity to weak water features that may only produce significant deviation from a blackbody at the center of this absorption band.Figures 3 and 4 show the observed HST water feature strengths for the sample of 19 hot Jupiter emission spectra compared to those of the models.Supplementary Table 2 lists the value of S H 2 O for each planet.Figure 3 shows that the observed HST/WFC3 feature strengths generally fall within the region of parameter space spanned by the models, with almost all of the planets fully within the predicted spread of the models.The models considered here assume elemental abundance ratios that fall within the range of commonly expected outcomes from planet formation models [18][19][20][21] .We find that varying parameters in these simple models can explain the observed hot Jupiter population without having to appeal to less likely outcomes of planet formation (e.g., C/O> 1 18,19,21 ) or exotic chemistry.
In order to demonstrate the difference between models of hot Jupiter atmospheres (which are primarily irradiated from above by their host stars) and self-luminous atmospheres (which are primarily heated from below by the object's interior), we created a separate model grid of self-luminous, cloudfree models using the same Sc-CHIMERA code setup.These models used an identical set of parameters to the hot Jupiter models, with the exception of irradiation from within the body instead of from an exterior star.The fiducial self-luminous models had g = 1000 m/s 2 , M H = 0.0, and C O = 0.55.We also created models with metallicities of M H = −1.0 and 1.0.Additionally, while a gravity of g = 1000 m/s 2 is typical for a brown dwarf, we created grids with g = 100 m/s 2 and g = 10 m/s 2 for direct comparison to lower-gravity hot Jupiters.
We used Equation (1) to derive water feature strengths for the grid of self-luminous models.Figure 3 shows these self-luminous water feature strengths compared to those from the hot Jupiter models and hot Jupiter observations, as well as water feature strengths derived from the HST/WFC3 brown dwarf spectra presented in Manjavacas et al. (2019) 16 .The self-luminous models generally show very distinct spectra from the hot Jupiter models.In particular, the self-luminous models consistently show negative values of S H 2 O indicating absorption features across the full range of temperatures modeled.This is because the atmospheric thermal inversions which produce emission features can only appear in atmospheres primarily heated from above (e.g., 33 ).Additionally, at temperatures below 2000 K where both sets of models show absorption features, the self-luminous models show consistently deeper features than the hot Jupiter models.We find this is due to the self-luminous models generally showing steeper T-P profiles, and therefore deeper absorption features, than the hot Jupiter models.We find that the hot Jupiter data are discrepant from all of the self-luminous models at ≥ 10σ significance.However, at the low temperatures of observed brown dwarfs (photospheric T < 2000 K), we find that the brown dwarf and hot Jupiter observations show similar water feature strengths.This is likely due to clouds, which can act to mute water feature strengths in both hot Jupiters and brown dwarfs at T < 2000 K (see Supplementary Figure 4).This muting of water features makes the brown dwarf population appear consistent with the cloud-free hot Jupiter models between 1 and 2 microns as shown in Figure 3.This is because the effect of clouds in muting brown dwarf absorption features is degenerate with the steepness of their T-P profiles, as shown previously by Burningham et al. (2017)  34 .Therefore, a cloudy brown dwarf with a steep T-P profile can appear to have a similar water feature strength as a model for a cloudfree hot Jupiter with a less steep T-P profile.However, we emphasize that at temperatures T > 2000 K, the population of observed hot Jupiters show emission features and featureless spectra (S H 2 O ≤ 0), which are consistent with our models of hot Jupiter atmospheres but inconsistent with the absorption features shown by self-luminous objects.We also note that the scatter in the observed brown dwarf population is likely due to variation in gravity and modest variation in metallicity, as changes in the gravity can influence the cloud cover.Brown dwarfs span a much wider range of gravities than hot Jupiters but are generally expected to have much smaller differences in atmospheric composition, particularly C/O ratio 35 .Although the observed population of hot Jupiter emission spectra generally matches the trends in our hot Jupiter model predictions, we find that no single model track is the best fit for all 19 of the observations.When taking each model track individually, the data are discrepant from each track at ≥ 2σ significance.The fact that different data sets are best fit by models with different parameters suggests there may be one or more parameters varying between the planets.To determine which parameters can most easily explain the scatter in the observed data, we examined the water feature strength variation we could achieve through changing each of our model parameters individually.Figure 4 shows water strengths for each individual model we examined.We found that the stellar effective temperature, planet gravity, and extent of internal heating had relatively small impacts on the predicted water feature strengths throughout the range of temperatures of the hot Jupiter population.The relatively small impact of the gravity on hot Jupiter water feature strengths is notably different than the large impact of gravity on the feature strengths of observed brown dwarfs.This difference is because the observed hot Jupiters generally show a much smaller span of gravities than the wide range represented in the population of observed brown dwarfs.Additionally, the models with TiO/VO opacity removed at different temperatures could only account for some of the scatter at intermediate temperatures and could not explain scatter at the highest or lowest temperatures, where we have observed the most precise secondary eclipse spectra.However, changing the atmospheric metallicity and C/O ratio had a significant impact on the predicted HST/WFC3 water feature strengths.We found the observed scatter could be explained if the planets have atmospheric metallicites between 0.03 − 30x solar and C/O ratios between 0.01-0.85(0.02 − 1.5x solar).With the current data we are unable to compare each planet's specific atmospheric composition to this prediction, as even the most detailed HST  16 .
secondary eclipse observations only constrain the metallicity to within 0.5 dex and often cannot constrain the C/O ratio, or can only place an upper limit (e.g., 4,6,26 ).However, such variation is expected from planet formation models 18,19 and has been suggested by a handful of transmission spectra studies (e.g., 17,21 ).The scatter we observed in emission spectra lends further support to the concept of compositional diversity among hot Jupiters.
Our hypothesis that hot Jupiters show compositional diversity can be tested through high-precision observations that cover more of the key O-and C-bearing molecules than are included in existing data sets (e.g, H 2 O, CO, CO 2 , and CH 4 ).Such observations will be possible with the upcoming James Webb Space Telescope (JWST) 36 and stabilized, highresolution spectrographs on large ground-based telescopes that have broad wavelength coverage 37 .Simultaneous detection of multiple molecules would lead to more precise constraints on metallicities and carbon-to-oxygen abundance ratios (and additional elemental ratios including nitrogen, etc.) 38 .Beyond testing our hypothesis, more precise compositional constraints on exoplanet atmospheres would inform our understanding of the formation and evolution processes that have produced the diverse planetary systems revealed over the last 25 years.In each case all other parameters are held fixed at the fiducial model values.The error bars include uncertainties in the stellar effective temperature.We found that changing the stellar temperature, planetary gravity, and internal heating in our models had little impact on the derived water feature strengths, and changing the TiO/VO only had an impact at intermediate temperatures, but changing the atmospheric C/O ratio and metallicity can explain the diversity of observed secondary eclipse spectra.

New Observations and Data Reduction
We reduced and analyzed HST/WFC3+G141 spectra of six planets.At the time this study was begun, these were all of the remaining secondary eclipse data sets in the HST archive that had not been published yet.Since we began this project, results for three planets have been published 6,[39][40][41] .In all of these cases, our reductions produced spectra consistent with the published results.Supplementary Table 3 lists the details of these observations, which included single eclipses of HAT-P-41b, KELT-7b, WASP-74b, WASP-76b, and WASP-79b, as well as five eclipses of WASP-121b.
We reduced the data using the data reduction pipeline described in Kreidberg et al. (2014)  42 .We used an optimal extraction procedure 43 and masked cosmic rays.To subtract the background out of each frame, we visually inspected the images to find a clear background spot on the detector and subtracted the median of this background area.The uncertainties on the measurements were determined by adding in quadrature the photon noise, read noise, and median absolute deviation of the background.
Following standard procedure for HST/WFC3 eclipse observations, we discarded the first orbit of each visit.The spectra of each planet were binned into 14 channels at a resolution R ≈ 30 − 40.We also created a broadband white light curve for each planet by summing the spectra over the entire wavelength range.
We fit both the white light curves and spectroscopic light curves with a model in the form where M(t) is the modeled flux, E(t) is an eclipse model found using batman 44 , and the rest of the equation is a systematics model based on Berta et al. (2012) 45 .In this systematics model, c is a normalization constant, s is a scaling factor to account for an offset in normalization between scan directions, v is a visit-long linear trend, t vis is the time since the beginning of the visit, r 1 and r 2 are the amplitude and time constant of an orbit-long exponential ramp, respectively, and t orb is the time since the beginning of the orbit.For the white light curves, the free parameters in the eclipse model were the mid-eclipse time T 0 and the planet-to-star flux ratio F p /F s .For the spectroscopic light curves, the mid-eclipse time T 0 was fixed to the best-fit value from the white light curve and the only free parameter in the eclipse model was the planet-to-star flux ratio F p /F s .The single eclipses observed for most of these planets had poor coverage of ingress and egress, so they could not constrain parameters such as the inclination to the level of precision provided by previous observations.Therefore, following best practices from previous studies (e.g., 26,[39][40][41]46 ), all other eclipse parameters were fixed to the literature values listed in Supplementary Table 4. Fo the systematics model, c, v, and s were allowed to vary between visits, while r 1 and r 2 were fixed to the same values for all visits.Four of the data sets (for HAT-P-41b, WASP-74b, WASP-79b, and WASP-121b) only used forward scanning instead of bi-directional scanning, so for these observations we fixed s = 1.The first secondary eclipse observation for WASP-121b occurred two years before the other four observations and showed significant differences in the ramp shape, so we allowed this first eclipse to be fit with different values of r 1 and r 2 than the other four visits.
The data sets for WASP-76b and WASP-79b showed additional correlated noise after applying this systematic model, so for these two data sets we tested adding an additional quadratic term to the visit-long trend.While adding this additional term was able to correct for the correlated noise, it introduced strong degeneracies between the fit parameters and the planet-to-star flux ratio.In order to avoid these degeneracies, we fit for only a linear visit-long trend in our final fit and used the divide-white method to correct for the additional noise 42 .
We estimated the parameters with a Markov Chain Monte Carlo (MCMC) fit using the emcee package for Python 47 .The final secondary eclipse spectra for all of the planets are shown in Supplementary Figure 2, and the planet-to-star flux ratio in each wavelength bin is listed in Supplementary Table 5.The white light curves had reduced chi-squared values between 1.9 < χ 2 ν < 15.2.The spectroscopic light curves generally achieved photon-limited precision, with ≈ 90 % of the light curves having reduced chi-squared values between 0.7 < χ 2 ν < 2.0.However, occasional individual spectroscopic light curves had higher reduced chi-squared values between 2.0 < χ 2 ν < 3.4.Therefore, before fitting each spectroscopic light curve we rescaled the uncertainties by a constant factor such that each light curve had χ 2 ν = 1 to give more conservative error estimates.
WASP-76 has a companion star whose spectrum is blended with that of the primary star in the WFC3 data.We corrected for the presence of this companion star using the following equation: where F corr is the corrected planet-to-star flux ratio in a given bandpass, F obs is the observed flux ratio in that bandpass including the companion star contamination, F B is the flux of the companion star in that bandpass, and F A is the flux of the primary star in that bandpass.We used ATLAS models 48 with temperatures of 6250 K and 4824 K to represent the primary star and the companion star, respectively 49 .

Reanalysis of Kepler-13Ab
In addition to the six new data reductions described above, we performed a reanalysis of the emission spectrum of Kepler-13Ab.The details for the two observed secondary eclipses of Kepler-13Ab are listed in Supplementary Table 3.These data were also reduced using the data reduction pipeline described in Kreidberg et al. ( 2014) 42 , and we again discarded the first orbit of each visit.We additionally discarded 14 of the 1008 observed spectra because they showed anomalously low fluxes in the broadband white light curve compared to the rest of the spectra.The spectrum was binned into 14 channels at a resolution of R ≈ 30 − 40.
The spectrum of Kepler-13Ab was observed in stare mode.Stare mode observations commonly show one or more of three types of systematics: a visit-long trend, an L-shaped hook trend over an individual orbit, and thermal breathing 50 .We tested including all of these components in our fit and found that only a visit-long trend was necessary to explain the systematics.Therefore, we fit both the white light curve and the spectroscopic light curves with a model in the form Following our method for the other data sets, the free parameters in the white light curve fit were c, v, T 0 , and F p /F s .For the spectroscopic light curves, the free parameters were c, v, and F p /F s , and T 0 was fixed to the best-fit value from the white light curve.All other eclipse parameters were fixed to the literature values listed in Supplementary Table 4.The parameters c and v were allowed to vary between visits.We estimated the parameters with a Markov Chain Monte Carlo (MCMC) fit using the emcee package for Python 47 .The final secondary eclipse spectrum for Kepler-13Ab is shown in Supplementary Figure 2, and the planet-to-star flux ratio in each wavelength bin is listed in Supplementary Table 5.The white light curve had a reduced chi-squared of χ 2 ν = 1.33.The spectroscopic light curves generally achieved photon-limited precision and had reduced chi-squared values between 0.94 < χ 2 ν < 1.08.

Observed Dayside Temperatures and Water Feature Strengths
As described in the main text, we determined the water feature strength (S H 2 O ) of each observed exoplanet using Equation (1).In order to ensure that S H 2 O would produce the same value for identical planets orbiting different stars, we first subtracted out the stellar signal from the observed data.We used ATLAS models 48 to calculate the stellar flux in our defined out-ofband and in-band regions.We then calculated the planetary flux in each of these regions using the equation where D is the observed secondary eclipse depth in each bandpass, F s is the stellar flux from the ATLAS model, and R p and R * are the planetary and stellar radius, respectively.We measured the dayside temperature of each observed planet by fitting a blackbody to the "out-of-band" regions indicated in Supplementary Figure 1.Similar to previous studies 51 , we found a linear relationship between this observed dayside temperature and the planetary irradiation temperature given by T day = 0.807 +0.008 −0.004 T irr + 71 +25 −8 , where T irr = T e f f R * /a is the irradiation temperature and a is the semi-major axis.

Model Grid
We created a new grid of self-consistent, 1D hot Jupiter models to compare their emission spectra to the population of observed planets.These models were generated using the Sc-CHIMERA code (validated against established brown dwarf models 52 and analytic models 27 ) assuming cloudfree, radiative-convective-thermochemical equilibrium atmospheres.The models' assumption of chemical equilibrium is likely a good approximation for the highly irradiated planets that make up the majority of our observed population 53 .
A two stream source function technique 54 is employed to solve for the planetary thermal fluxes at each atmospheric level (under the hemispheric mean approximation).We modeled the stellar flux via a standard two stream approximation (for both direct and diffuse fluxes, under the quadrature approximation) assuming cosine incident angle of 0.5, utilizing the PHOENIX models for the stellar spectra 55 .A Newton-Raphson iteration 56 is used to determine the temperature at each model layer which ensures zero net flux divergence.We include absorption cross-sections from 0.  [57][58][59][60][61] , H − bound-free and free-free 62,63 , and H 2 /He Rayleigh scattering, and additional UV opacities for CO, SiO, and H 2 59 .Pre-computed cross-sections were converted into correlated-K coefficients at a spectral resolution of 250 using a 10 point double Gauss quadrature (with half covering the top 5% of the correlated-K cumulative distribution function) with mixed-gas optical depths computed using the random-overlap resort-rebin framework (e.g., 64,65 ).Thermochemical equilibrium molecular abundances were computed using the NASA CEA Gibbs free energy minimization code 66 combined with elemental-rain out due to condensate formation (all major Si, Fe, Mg, Ca, Al, Na, and K bearing condensates are included) given the Lodders et al. (2009)  67 elemental abundances.
We parameterized the model atmospheres with a set of five parameters: the stellar effective temeprature (T e f f ), the planetary gravity (g), the planetary metallicity ( M H ), the planetary carbon-to-oxygen ratio ( C O ), and the planetary internal temperature (T int ).Our fiducial models had the following parameter values: T e f f = 5300 K, g = 10 m/s 2 , M H = 0.0, C O = 0.55, T int = 150 K. Models at different irradiation temperatures were created by re-scaling the incident stellar spectrum / (the PHOENIX model for a given stellar effective temperature) by the ratio of the desired irradiation temperature to the bolometric temperature of a planet at 0.05 AU around a 1 solar radius star.We created hot Jupiter models with irradiation temperatures between 500 − 3600 K, with step sizes of 50 − 200 K.Following Lothringer & Barman (2019) 29 , we calculate the absorption mean opacity κ J and the Planck mean opacity κ B at a pressure of 10 −2 bar as a function of equilibrium temperature for our fiducial models.The absorption mean opacity at a given pressure P is given by where κ λ is the monochromatic true absorption coefficient, J λ is the mean intensity at a given wavelength, and T is the local temperature in the planet's atmosphere.The Planck mean opacity is given by where B λ (T ) is the Planck function.The absorption mean opacity represents the efficiency with which the atmosphere can absorb photons, while the Planck mean opacity represents the efficiency with which the atmosphere can emit photons 29 .Therefore, the ratio κ J /κ B describes the relative efficiency of stellar absorption vs. thermal re-radiation, and a ratio κ J /κ B > 1 generally indicates the presence of a thermal inversion in the T-P profile.The hot Jupiters in this study can generally be thought of as emitting most of their radiation at near-infrared wavelengths, whereas incoming starlight from their host stars peaks at visible wavelengths.Therefore, increasing the amount of molecules such as TiO that are optically active at visible wavelengths will increase κ J , and increasing the amount of molecules such as H 2 O that are optically active at near-infrared wavelengths will increase κ B .We also created subset grids as a function of irradiation temperature where a single parameter dimension was varied while all other parameters were held fixed to their fiducial model values (no cross-variance).We examined models with a stellar T eff = 3300 K, 4300 K, 6300 K, 7200 K, and 8200 K; g = 1 m/s 2 and 100 m/s 2 ; M H = −1.5 and 1.5; and C O = 0.01 and 0.85.For models with different metallicities, elemental abundance ratios were held constant while the overall metallicity was re-scaled relative to H. We also created a model grid where the internal temperature varies with the planetary irradiation temperature following Equation 3 in ref ( 30 ).Individual model tracks with irradiation temperature for each of these variations are shown in Figure 4.
Opacity from gaseous TiO/VO is theorized to be a driving force behind the transition between uninverted hot Jupiter atmospheres with monotonically decreasing T-P profiles and atmospheres containing thermal inversions 2 .Some previous observations of hot Jupiters have suggested that vapor TiO/VO may not be present in high-temperature atmospheres if it is condensed in cooler parts of the atmosphere (e.g., 32 ).This process, known as cold-trapping, effectively works to remove TiO/VO from places in the atmosphere where vaporized TiO/VO would be expected to be present in equilibrium.In order to study the impact of potential cold-trapping, we created models where the TiO and VO opacities are artificially set to zero until a given temperature threshold.We tested models where TiO/VO opacity is zeroed out for temperatures below 2000 K, 2500 K, and 3000 K.These tracks are also shown in Figure 4. Supplementary Figure 3 shows the best-fit model from this complete grid for each individual data set, and Supplementary Table 6 lists the reduced chi-squared values for these best-fit models.We find that the model grid is generally able to produce good fits to the data, with the best fits to all but two data sets having reduced chi-squared values below 2.6.However, different data sets are best fit by models with different values for the atmospheric metallicity and C/O ratio, which suggests their atmospheres may have diverse compositions.
Recent studies have suggested clouds may have an impact on the strength of molecular absorption features observed in thermal emission (e.g., 11,68 ).To test the impact the presence of clouds would have on the trends in our models, we created two cloudy models.We used the cloud model of Ackerman & Marley (2001) 69 , as implemented by Mai & Line (2019) 70 .Both models had a constant vertical mixing strength of 10 8 cm 2 /s using the Zahnle et al. (2016) 71 timescale prescription.We tested models with sedimentation efficiencies of f sed = 0.1 and 1.0.These models are shown compared to the fiducial model in Supplementary Figure 4. Adding clouds acts to weaken the water feature strengths below a dayside temperature of about 2000 K, with a smaller f sed leading to more effective weakening.While clouds may provide a potential explanation for the weak water feature strength of HD 189733b, the lowest-temperature hot Jupiter in our population study, we find that including clouds can not generally explain the scatter we see in water feature strengths and has no impact on the feature strengths above T day = 2000 K. Our results agree with those from general circulation models, which also show that clouds have little to no impact at temperatures above ≈ 2000 K 72,73 .

Self-Luminous Object Models
In order to demonstrate the difference between models of hot Jupiter atmospheres (which are primarily irradiated from above by their host stars) and self-luminous atmospheres (which are primarily heated from below by the object's interior), we created a separate model grid of cloud-free selfluminous object models using the same Sc-CHIMERA code setup.We parameterized the self-luminous model atmospheres with a set of four parameters: the effective temperature (T e f f ,bd ), the gravity (g), the metallicity ( M H ), and the carbon-to-oxygen ratio ( C O ).These models thus used an identi-/ cal set of parameters to the hot Jupiter models, with the exception of irradiation from within the self-luminous body instead of from an exterior star.The fiducial self-luminous models had g = 1000 m/s 2 , M H = 0.0, and C O = 0.55.We created models with effetive temperatures between 1000 − 2800 K with a step size of 200 K.We also created grids with metallicities of M H = −1.0 and 1.0.Additionally, while a gravity of g = 1000 m/s 2 is typical for a brown dwarf, we created grids with g = 100 m/s 2 and g = 10 m/s 2 for direct comparison to lower-gravity hot Jupiters.We note that the self-luminous models are cloud-free and therefore likely overpredict water feature strengths at temperatures below ≈ 2000 K.

Figure 1 .
Figure 1.Secondary eclipse spectra showing the planet-to-star flux ratio (F p /F s ) as a function of wavelength for all 19 hot Jupiters considered in this study.Data sets are colored by dayside temperature, which is measured as described in the Methods and shown by the colorbar.Solid lines indicate interpolations from our solar composition fiducial model grid (see the Methods section for a description), while dashed lines indicate best-fit blackbodies.Shaded regions indicate the "out-of-band" and "in-band" regions used to calculate the water feature strength (S H 2 O ) for each observed secondary eclipse spectrum.Note that for several data sets, the error bars are smaller than the point size.

Figure 2 .
Figure 2. (a) Temperature-pressure (T-P) profiles and (b) resulting dayside planet fluxes for the fiducial model grid, which covers approximately the same range of temperatures as spanned by the observations.The full model specifications are detailed in the Methods.The fiducial model uses a 5300 K stellar effective temperature, a solar composition planetary atmosphere ( M H = 0.0 and C O = 0.55), a planetary gravity of 10 m/s 2 , and a planet internal temperature of 150 K. Blue and yellow lines show models with the coolest and warmest irradiation temperatures, respectively.The model grid spans a range of irradiation temperatures between 500 − 3600 K, with step sizes of 50 − 200 K.For clarity, only every other model in the grid is shown here.Grey shaded bands indicate the "out-of-band" and "in-band" regions used to calculate the water feature strength (S H 2 O ) for each model.(c)Ratio of the absorption mean opacity (κ J ) to the Planck mean opacity (κ B ) as a function of equilibrium temperature (T eq ) in these models at a pressure of 10 −2 bar, assuming zero albedo and full redistribution.This ratio describes the relative efficiency of heating vs. cooling in the models29 , and a ratio of κ J /κ B > 1 generally indicates the presence of a thermal inversion in the T-P profile.This panel is plotted with temperature on the y-axis for ease of comparison to Figure3.

Figure 3 .
Figure 3. HST water feature strength diagram comparing observed secondary eclipse spectra to the model predictions in Figure 2. The y-axis shows the temperature of a blackbody fit to the "out-of-band" regions defined in Supplementary Figure 1, which is the observed dayside temperature T day .The x-axis shows the strength of the observed feature in the water band at 1.4 µm compared to this blackbody, as defined by Equation (1).Featureless, blackbody-like spectra have S H 2 O = 0 and absorption/emission features have positive/negative values of S H 2 O , respectively.The gray line and points show the fiducial hot Jupiter models pictured in Figure 2. The light gray shaded region shows the full range of hot Jupiter model predictions assuming different values for the stellar effective temperature; the temperature where TiO opacity becomes important; and the planet gravity, C/O ratio, metallicity, and internal heat.Similarly, the brown line and points show the fiducial self-luminous object models, and the tan shaded region shows the full range of self-luminous models assuming different values for the planet gravity and metallicity.Colored points with 1σ error bars show all planets with HST/WFC3 spectra, and boxes around planet names indicate new data reductions in this publication.The color scale indicates the planetary equilibrium temperature.The error bars include uncertainties in the stellar effective temperature.Brown points show isolated brown dwarf spectra observed with HST/WFC3 from Manjavacas et al. (2019) 16 .

Figure 4 .
Figure 4. Diagrams illustrating the change in HST water feature strength from models with different parameters.All diagrams show the observed hot Jupiter data as black points with 1σ error bars, while the lines show tracks for models with varying C/O ratio (a), metallicity (b), stellar temperature (c), gravity (d), internal heating (e), and the temperature to which TiO/VO opacity were ignored (f).In each case all other parameters are held fixed at the fiducial model values.The error bars include uncertainties in the stellar effective temperature.We found that changing the stellar temperature, planetary gravity, and internal heating in our models had little impact on the derived water feature strengths, and changing the TiO/VO only had an impact at intermediate temperatures, but changing the atmospheric C/O ratio and metallicity can explain the diversity of observed secondary eclipse spectra.