Five key exoplanet questions answered via the analysis of 25 hot Jupiter atmospheres in eclipse

Population studies of exoplanets are key to unlocking their statistical properties. So far the inferred properties have been mostly limited to planetary, orbital and stellar parameters extracted from, e.g., Kepler, radial velocity, and GAIA data. More recently an increasing number of exoplanet atmospheres have been observed in detail from space and the ground. Generally, however, these atmospheric studies have focused on individual planets, with the exception of a couple of works which have detected the presence of water vapor and clouds in populations of gaseous planets via transmission spectroscopy. Here, using a suite of retrieval tools, we analyse spectroscopic and photometric data of 25 hot Jupiters, obtained with the Hubble and Spitzer Space Telescopes via the eclipse technique. By applying the tools uniformly across the entire set of 25 planets, we extract robust trends in the thermal structure and chemical properties of hot Jupiters not obtained in past studies. With the recent launch of JWST and the upcoming missions Twinkle, and Ariel, population based studies of exoplanet atmospheres, such as the one presented here, will be a key approach to understanding planet characteristics, formation, and evolution in our galaxy.


INTRODUCTION
More than 4700 exoplanets are currently known and for about 80 of these we have atmospheric data recorded with the Hubble Space Telescope (HST), Spitzer Space Telescope and ground-based instruments. To date, most data-oriented atmospheric analyses have focused on individual targets and only a very limited number of studies are available on populations of exo-atmospheres observed with the transit or eclipse technique (Sing et al. 2016a;Tsiaras et al. 2018;Pinhas et al.  . Best fit spectra of the 25 planets considered in this study for the HST+Spitzer observations in eclipse. Individual analyses and additional retrievals, including a simpler Blackbody fit, can be found in the Appendix. The planets are ordered in increasing atmospheric temperatures, which are traced by the colours (from blue to red). The retrieved atmospheric temperature is obtained by weighting the retrieved thermal profile by the contribution function. Unlabeled x-axes range from 1.1 µm to 2 µm.
highlighted the potential impact of metal hydrides and metal oxides in regulating the thermal properties of the day-side of hot Jupiters (Burrows 2014;Spiegel et al. 2009;Hubeny & Burrows 2009;Nugroho et al. 2017). Theoretical predictions (Lodders 2002;Fortney et al. 2008) have suggested that, in analogy to Brown dwarfs, planets hotter than ∼ 1700 K might have metal hydrides and metal oxides in the gaseous form: these molecular species are excellent absorbers in the optical /near infrared and might cause a hot layer in the atmosphere, and therefore a thermal inversion. By contrast, planets colder than ∼ 1700 K might not exhibit such thermal inversions, as metal hydrides and metal oxides have been sequestered into condensates.
2) Are the eclipse spectra of the hottest atmospheres consistent with blackbody curves? Recent theoretical studies (Bell & Cowan 2018;Lothringer et al. 2018;Parmentier et al. 2018) have suggested H − being an important opacity source, generated by the thermal dissociation of H 2 and H 2 O at very hot temperatures (T > 2500 K). Eclipse spectra of very hot atmospheres would therefore resemble blackbody curves due to the continuous shape of H − emission, and the absence of other molecules which have been dissociated. Observational evidence in favor of those predictions were obtained by studies of individual planets with HST and Spitzer (Mansfield et al. 2018;Arcangeli et al. 2019;Kreidberg et al. 2018), and more recently in the HST population analysis from Mansfield et al. (2021). However, many other works (Haynes et al. 2015;Mikal-Evans et al. 2020; have found spectral signatures in similarly hot atmospheres. 3) What is the dayside-terminator contrast in exoplanet atmospheres? Planetary atmospheres present vertical and horizontal inhomogeneities: this is particularly true for tidally locked planets, for which large day-night thermal gradients, sometimes in excess of 1000 K, and asymmetric circulation patterns have been predicted (Cho et al. 2003;Showman et al. 2010;Cowan & Agol 2011;Roth et al. 2021). Recent simulations using pseudo-spectral methods (Skinner & Cho 2021a,b;Cho et al. 2021) have revealed the complex and turbulent nature of these atmospheres, displaying highly dynamic small and large scale storms that develop and evolve in time. Some studies (Tan & Komacek 2019; Roth et al. 2021) have hypothesised that optical absorbers and thermal dissociation/recombination processes would impact the atmospheric dynamics. Namely, if H 2 dissociates at the day-side and recombines at the night-side, this would increase the energy transport, thus reducing the day-night temperature contrast. Mansfield et al. (2020a) provided observational evidence of these effects in the Spitzer phase-curve of KELT-9 b. On the other hand, optical absorbers such as TiO, VO and FeH are more thermally stable, so they would not contribute much to this effect. The comparison between eclipse and transit observations might constrain these dynamical processes.

4)
Are metallicity and C/O viable observables to understand planet formation? Hot Jupiters are thought to form via a three step process (Mizuno 1980;Bodenheimer & Pollack 1986;Ikoma et al. 2000): solid core accretion, runaway gas accretion and migration. The core accretion is believed to occur in the outer regions of a protoplanetary disk, where the abundance of solid materials leads to a rapid growth of a planetary core before disk gas dispersal. If that is the case, the composition of giant exoplanets, which is a direct outcome of these planetary formation processes, is predicted to be sub-stellar in heavy elements such as C, O, and refractory elements, because most of the heavy elements would be sequestered in the cores. The inferred super-stellar bulk metallicities of Jupiter, Saturn and warm exo-Jupiters (Saumon & Guillot 2004;Miller & Fortney 2011;Thorngren et al. 2016;Welbanks et al. 2019), however, suggest a more complicated picture, where heavy elements are also captured through some other processes, probably, during gas accretion and/or migration (Hasegawa et al. 2018;Shibata & Ikoma 2019;Shibata et al. 2020;Turrini et al. 2021). Also, while difficult to determine today, the C/O ratio of exoplanets may place constraints on where and how the giant planets collect gas and solids in evolving protoplanetary disks (Madhusudhan et al. 2017;Mordasini et al. 2016;Brewer et al. 2017;Booth et al. 2017;Eistrup et al. 2018;Cridland et al. 2019). Constraints on those two parameters would significantly improve our understanding of planetary formation.

5)
Can refractory elements help us understand exoplanet formation? In addition to the metallicity and C/O ratio, other elemental ratios, such as N/O, S/O (Turrini et al. 2021) or even refractory elements (Lothringer et al. 2020), may help constrain planet formation scenarios. Their potential, however, remains unexplored by observational studies, as their tracers are more difficult to detect. While HST is not particularly sensitive to N and S bearing species, refractory elements such as TiO, VO and FeH have been detected previously in eclipse spectra. Day-side Temperature (K) Figure 2. Correlations between the retrieved abundances of optical absorbers (TiO, VO, FeH and H − ) and the thermal gradient in the atmosphere. The optical absorber abundances are estimated using the weighted average of the individually retrieved abundances for the detected molecules. If no optical absorber is detected, only a upper limit is available. The thermal gradient is the difference between atmospheric temperatures at 5% and 95% of the atmospheric contribution function. Colors on the datapoints represent averaged retrieved temperatures. Overall, Planets with inverted thermal profiles have temperatures above 2000K and possess optical absorbers (shaded orange region). Planets without thermal inversions have temperatures below 2500K and do not possess optical absorbers (shaded green region).
observational data to spectra as 'data reduction', to distinguish this process from the analysis/interpretation of the atmospheric spectrum. Only two exoplanets, Kepler-13A b and WASP-33 b, were taken as is from the literature. WASP-33 b orbits a pulsating δ-Scuti star and Kepler-13A b is part of a triple star system. Reduction of the data for these particular targets requires a particularly careful modelling of the host stars, which is not explored in this study. For our retrieval analyses we have used Alfnoor (Changeat et al. 2020a), a tool that extends the atmospheric retrieval capabilities of TauREx3 (Al- Refaie et al. 2021b) to populations of atmospheres. More technical details on our tools and methods are reported in Materials and Methods.
For each planet, we tested predefined scenarios, varying the molecular species included, the model assumptions, and whether Spitzer data are included. As potential biases in retrieval studies can arise from combining HST and Spitzer observations (Yip et al. 2020;Changeat et al. 2020b), we assessed the robustness of our results against possible biases by artificially modifying the Spitzer data and performing additional retrievals. More specifically, we repeated our analysis with Spitzer data offset by +100ppm, -100ppm, and doubled uncertainties.
The free retrieval runs assumed constant with altitude abundances and included molecular species such as H 2 O, CO, CO 2 , CH 4 in a 'reduced' run. Refractory species (TiO, VO and FeH) and H − were added in a 'full' run. Please note that H − absorbing properties are traced by the e − abundance (John 1988). Atomic and ionic species are not considered here as they do not absorb in the HST and Spitzer wavelength regions. In the rest of the paper, we use the term 'refractory' to refer to TiO, VO and FeH, and we use the term 'optical absorber' to refer to TiO, VO, FeH and H − . We have also attempted equilibrium chemistry retrievals, using the GGChem code (Woitke et al. 2018) and including all the supported species, on our entire population. We used the Bayesian Evidence to compare the results with the free runs. In all models, the temperature profile is determined using a heuristic N-point profile. Since the thermal profile varies with altitude, the reported value is the mean atmospheric temperature obtained in the retrievals, weighted by the atmospheric contribution function. It characterizes the thermal conditions of the atmosphere in the region probed by the observations, Note that this is not equal to the Blackbody temperature or the equilibrium temperature. For reference, best-fit Blackbody temperatures are provided in the Appendix D, which analyses the planets individually.
For a subset of 17 planets we have also analysed the G141 transit observations using Iraclis, to include additional information about the terminator region of those planets in our study.

RESULTS
The HST eclipse observations reduced with our Iraclis pipeline are shown in Figure 1, alongside the Spitzer photometric data. A summary of the observations considered in this paper and relative references are given in Table B1 of the Appendix. We also show in Figure 1 the best-fit spectra obtained with Alfnoor. The cooler planets in the sample, such as WASP-43 b, may exhibit H 2 O features in absorption (around 1.4 µm) and a steep increase in the spectrum towards longer wavelengths. Hotter planets, such as WASP-121 b, show emission features from optical absorbers (from 1.1 to 1.4 µm) and their near-Infrared spectra are more leveled.
While a detailed analysis of each planet is provided in Appendix D, we summarise our key findings in Table  1. Overall, after inspecting the retrieval results of our selected sample of 25 planets, we find that the HST-only spectra are difficult to interpret using free retrievals, as these may also converge to unphysical solutions. This behaviour can be explained by the narrow wavelength range covered by HST (see Appendix C), which leads to large degeneracies and does not support the complexity of an emission model where both the thermal profile and atmospheric chemistry have to be disentangled. When Spitzer IRAC observations at 3.6 µm and 4.5 µm are included, the thermal profiles can be retrieved more easily. In this case, a free chemical retrieval which includes only H 2 O, CO, CO 2 , CH 4 (reduced model) does not fit well the entire population of exoplanets. To improve our fit, we have to include other plausible absorbers such as refractory species (TiO, VO and FeH) and H − (full model). These species are referred as optical absorbers throughout this work. When Spitzer data is considered, the full model obtains the highest Bayesian Evidence (Jeffreys 1961), ln(E), compared to a simple blackbody fit or a reduced model. The full model is rejected in only one case, i.e. HD 209458 b (this case is developed further in Appendix B and in the planet's section). Equilibrium chemistry is expected to be a relatively valid assumption to describe planets in the 1000 to 2000 K range. However, atmospheric retrievals following equilibrium chemistry, while not strongly rejected by the Bayesian Evidence (e.g ∆ ln(E) < 3), are almost never favoured as best solutions. Possibly, subtle departures from equilibrium chemistry may be better captured by the free retrievals. Another possibility could be that the elemental abundances of refractory elements (Ti, V and Fe), for which their ratios (Ti/O, V/O and Fe/O) remain solar in our equilibrium tests, are enhanced (see key question 5 bellow).
Combining HST and Spitzer observations should be done carefully as it can introduce biases in retrievals (Yip et al. 2021). We tested the robustness of our conclusions by including offsets in the Spitzer data, or artificially increasing the photometric uncertainties. We find that while planets can be affected individually, our results on the entire population remain unchanged.
The transit observations and their associated fits are reported in Figure A1 of Appendix A. At the terminator, water vapour is recovered most of the times, independently from the planet's atmospheric temperature. There, atmospheric clouds are also detected in many planets. Concerning the key open questions listed in the Introduction, focusing more on the results of the full model with HST+Spitzer, we find: 1) Refractory molecules and H − correlate with thermal inversions.
We show in Figure 2 the weighted mean abundances of optical absorbers (TiO, VO, FeH and H − ) as a function of the retrieved thermal gradient in the atmospheres of all the planets. Maps displaying the retrieved abundances for each individual species as a function of the atmospheric temperature are available in the Appendix (Figures A3 to A6). We see in Figure 2 a clear correlation between the retrieved abundances of optical absorbers and the retrieved vertical thermal gradients: planets with positive thermal gradients show features associated with optical absorbers. While correlation does not imply causation, refractory elements (here TiO, VO and FeH) and H − are efficient absorbers of stellar light. In the hottest atmospheres (T > 2000 K), refractory elements are thermally stable and can remain in the gas phase. They could therefore provide a natural explanation for thermal inversions. At even higher temperatures, H 2 is thermally dissociated, leading to abundant H − opacity. H − also absorbs in the visible and could be another contributor to the retrieved thermally inverted profiles.
2a) Spectra of the hottest exoplanets do not resemble blackbodies. Figure 1 demonstrates that the hottest planets in our sample show features associated with molecular species. This is also confirmed by the Bayesian Evidence of retrievals including absorbing species as opposed to featureless blackbody curves as shown in the complementary figure A2 of the Appendix). The presence of refractory elements lead to spectral features that make the spectra deviate from the pure blackbody.
2b) Dissociation processes do occur in the hottest atmospheres. Day-side Temperature (K) Figure 3. Retrieved abundances of e − versus H2O at the day-side in our sample of 25 hot Jupiters recovered in the HST+Spitzer full runs. The colours on the datapoints indicate the retrieved atmospheric temperatures weighted by the contribution functions with legend in the colobar. The shaded green region indicates the predicted abundances from equilibrium chemistry at solar metallicity (M=1) and C/O (C/O=0.55) between 1 bar and 0.01 bar, which is around the region probed by the observations. In orange, the metallicity is increased to 100 times solar. In red, the metallicity is decreased to 0.1 times solar. In blue, the metallicity is solar and the C/O is increased to 1. The planets separate in three regimes: A solar hot Jupiter regime where water vapour is detected in moderately hot atmospheres associated with a decreasing thermal profile (dashed green); An ultra-hot Jupiter regime where thermal inversions and H − emission are detected in very hot atmospheres (dashed orange); A high C/O hot Jupiter regime where temperatures remain moderate, but where water is not detected (dashed blue).
Our results also confirm the presence of dissociation processes in the hottest atmospheres: Figure 3 indicates the apparition of H − opacity from the dissociation of H 2 and H 2 O. The trend appears to be driven by the retrieved atmospheric temperature: H 2 O is detected in the spectra of cooler atmospheres (high C/O and solar HJ) and H − becomes more dominant for atmospheres hotter than 2500 K (ultra HJ). As demonstrated by the different shaded regions of this figure, the metallicity and C/O ratio impacts the transition from non-dissociated to thermally dissociated regimes. Additional figures are available in the Appendix A, showing the temperature dependence for each of the investigated molecules.
3) The day-terminator contrast of exoplanets does not appear to correlate with temperature and/or dissociation.
In Figure A7 we compare the atmospheric temperature on the day-side, as retrieved from eclipse spectra, and the atmospheric temperature at the terminator, as retrieved from transit spectra, to the estimated equilibrium temperature. Overall, we find that the day-side is about 20% hotter than the equilibrium temperature, while the terminator is about 30% cooler as also expected by theoretical studies (Cho et al. 2003;Tan & Komacek 2019;Parmentier & Crossfield 2017). We observe some scatter across the population in the recovered temperature (see Figure A7 in the Appendix), which could be explained by uncertainties in the observation, intrinsic variability in the population or atmospheric temporal  . Day-terminator thermal gradient as a function of the e − abundance recovered using the free retrievals. Only planets that have both eclipse and transit spectra and for which the terminator temperature is constrained are shown in this plot.
variability (Cho et al. 2003(Cho et al. , 2015Komacek & Showman 2020;Cho et al. 2021;Skinner & Cho 2021a,b). Our retrieved temperature at the terminator is compatible with the predicted temperatures from General Circulation Models (GCMs) at the night-side. This confirms potential biases of transit observations towards lower retrieved temperatures, most likely due to 3D effects or the presence of clouds (Caldas et al. 2019;MacDonald et al. 2020;Skaf et al. 2020;. The presence of H − on the day-side does not appear to be correlated with an increased of the day-terminator thermal gradient, as seen in Figure 4. For instance, for CoRoT-1 b, the terminator is found hotter than the day-side. While our results have large uncertainties, our study cannot confirm that the dissociation and recombination of H 2 have large effects on day-terminator contrast (Tan & Komacek 2019;Roth et al. 2021), as demonstrated in other observational studies. To our knowledge, the only planet for which such effect has been reported is KELT-9 b (Mansfield et al. 2020a), but given the fact that we do not have transit observations for this planet, we are unable to verify this claim.

4)
Planets have solar to sub-solar water vapour abundances, which suggest two reservoirs for planet formation.
In most studies, metallicity and C/O ratio of exoplanets are inferred from the detection of water vapour and C-bearing species, or from chemical equilibrium. Here, we demonstrate that for the coolest planets, water vapour -if detected -can be used as a proxy for metallicity, allowing us to estimate the O/H ratio. This is shown in Figure A8, where the estimated metallicity from O/H is consistent with the equilibrium retrievals (see the inset). For these planets, the chemistry of H 2 O is relatively well understood and should be well described by chemical equilibrium. For the hotter planets (T > 2500 K ), water vapour dissociates, as shown in point 2), and H 2 O becomes a poor indicator of the planet metallicity. In many cases, inferring metallicity and C/O ratio from water only is degenerate. We attempted to constrain the abundances of Vertical Thermal gradient (K) Figure 5. Averaged refractory element (TiO, VO and FeH) abundances recovered in the HST+Spitzer full runs, weighted by their standard deviation, and plotted against the retrieved atmospheric temperature on the day-side. Green area: predictions from a solar composition and equilibrium chemistry between 1 bar and 0.01 bar. Orange area: the metallicity is assumed to be 100 times solar. Red area: the metallicity is assumed to be 0.1 times solar. Blue area: The metallicity is solar and C/O is assumed equal 1.
the carbon-bearing species and estimate C/O from our free retrievals. At the equilibrium, CH 4 is predicted to be scarce in hot atmospheres, while CO should be the main carbon-carrier. We detected evidence of CO only in HD 209458 b. Instead, in a few instances, such as in WASP-43 b or HAT-P-7 b, we found evidence for high abundances of CO 2 . These detections are likely driven by the stronger absorption of CO 2 in the 4.5 µm Spitzer channels and its additional absorption at 1.5 µm. We note that, except for HAT-P-7 b, performing retrievals without CO 2 in the cases where it is detected leads to detection of CO instead, albeit with a slightly lower Bayesian evidence (∆ ln(E) < 5). As such, we believe the recovered abundance for CO 2 are not reliable and might here be overestimated.
While directly inferring the C/O from C-bearing species remains difficult, the abundance of water vapour also affects this parameter, which is shown in Figure 3 by the blue shaded region. In this figure, three regimes are identified. The first regime (orange dashed line) includes the hottest planets for which H − is detected. We label those as ultra HJ. For cooler planets two reservoirs exist. The first reservoir (green dashed line) contains planets for which water vapour is detected. The recovered abundances are roughly consistent with solar predictions (see also Figure A3). There are exceptions, such as HD 189733 b and HD 209458 b, for which the data allows us to detect much lower water abundances. Note also that WASP-18 b and CoRoT-1 b have water abundances that are consistent with solar, but those planets also present a particularly high dissociation. The third reservoir (blue dashed line) encompasses planets which do not show the spectral feature of H 2 O and H − . The non-detection of water in those planets could be due to either a depletion of the molecule, or other molecules/clouds masking the water signal. If these atmospheres are indeed depleted of H 2 O, as also suggested in previous studies Welbanks et al. 2019;Pinhas et al. 2019), this provides strong constraints on metallicity and C/O ratio. The depletion could be explained by either an overall subsolar metallicity, or alternatively a high C/O (see blue region in Figure 3). When this information is combined with point 5), we argue that the latter is more probable. At present, the bulk metallicities and C/O ratios of exoplanets are difficult to constrain from the HST and Spitzer data only, partly due to the limited amount of tracers that our observations are sensitive to. Future telescopes will have the potential to characterise the metallicity and C/O for exoplanets very accurately and confirm these predictions.

5) A contradiction exists between the abundances of volatiles and refractory elements.
For the hotter planets, we can constrain refractory elements from the abundance of metal oxides and hydrides. As their individual abundances could be unreliable due to their overlapping features, we show in Figure 5 their averaged retrieved abundances weighted by standard deviation. Shaded regions demonstrate that refractory species (TiO, VO and FeH) are expected to condensate at the cooler temperatures (T < 1300K) and thermally dissociate for temperatures higher than 3500K. For intermediate temperatures, they can exist in a gaseous form and be detectable. For planets hotter than 2500 K, at least one metal oxide/hydride is present and, interestingly, the recovered weighted abundance is higher than the one predicted from equilibrium chemistry for solar metallicity values. However, it is compatible with 100 times solar. While we cannot entirely exclude systematic errors from missing molecules in our retrievals or other biases arising from model assumptions, this result contradicts our estimates of the metallicity from water vapour, which suggested solar to sub-solar water abundances (see Figure 3). If confirmed, our result suggests that the water depletion identified here and found in previous studies, would come from a depletion in oxygen rather than an overall sub-solar metallicity, which would have strong implications for planetary formation. We note that a recent study (Welbanks et al. 2019) reached the same conclusion when using Alkali as secondary tracers of metallicity, which were also found to have super-solar abundances.

DISCUSSION AND CONCLUSION
We have presented here the first retrieval population study of exoplanet atmospheres observed in eclipse with the Hubble and Spitzer Space Telescopes. Our sample includes 25 hot gaseous planets. When combining the HST data with the available Spitzer data, the trends in our population are stable to changes of ±100 ppm in the Spitzer data, and an increase of observed noise by a factor two. Overall, we have found that: • the coolest planets in our sample (T < 2000 K) have non-inverted thermal profiles with signatures of water absorption.
• The hottest planets in our sample (T > 2000 K) have inverted thermal profiles with signatures from thermal dissociation (H − ) and refractory species (TiO, VO or FeH). Their spectra in the HST wavelength range are not consistent with a simple blackbody emission.
• The dayside-terminator thermal gradient is not found to be correlated with equilibrium temperature or H − opacity.
• Metallicity and C/O are difficult quantities to constrain from free retrievals of current data. Our results suggest water abundance being solar to sub-solar in the sample analyses. If confirmed, this result would have important implications on planet formation.
• Metal oxides and hydrides are found in excess of solar abundances in the hotter planets, thus contrasting with our results on the water abundances. If confirmed, this result would inform our current understanding of giant planets' formation.
• For a number of planets in our sample, equilibrium chemistry retrievals are not the preferred solution. While we cannot strongly reject equilibrium chemistry, this suggests that disequilibrium mechanisms might be important and highlight the importance of carrying unbiased free retrieval approaches. Evidence for disequilibrium processes have also been found in previous transiting studies (Roudier et al. 2021;Kawashima & Min 2021).
Population studies, such as the one presented here, pave the way to future studies based on the next generation of space telescopes. In the next decade, JWST (Gardner et al. 2006), Twinkle (Edwards et al. 2019b and Ariel (Tinetti et al. 2018 will provide atmospheric data for thousands of diverse worlds, enabling the study of chemical regimes, circulation patterns and formation mechanisms well beyond the parameter space explored here.  Molecules that are only detected when Spitzer is added are marked with a star symbol ( * ). We also report if clouds were found at the planet's terminator and we indicate the retrieved temperature. The stated temperature is the retrieved atmospheric temperature, weighted by the contribution function. This temperature is not equal to the temperature obtained by the simpler Blackbody fit, which we provide in the individual planet analyses in Appendix D.

DATA AND MATERIALS AVAILABILITY
This work is based upon observations with the NASA/ESA Hubble Space Telescope, obtained at the Space Telescope Science Institute (STScI) operated by AURA, Inc. The raw data used in this work are available from the Hubble Archive which is part of the Mikulski Archive for Space Telescopes. This work is also based in part on observations made with the Spitzer Space Telescope, which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under a con-tract with NASA. We are thankful to those who operate these telescopes and their corresponding archives, the public nature of which increases scientific productivity and accessibility (Peek et al. 2019).
K e p le r-1 3 A b H D 1 8 9 7 3 3 b H D 2 0 9 4 5 8 b Tr E S -3 b W A S P -4 b W A S P -1 2 b W A S P -1 8 b W A S P -1 9 b W A S P -3 3 b W A S P -4 3 b W A S P -7 4 b W A S P -7 6 b W A S P -7 7 A b W A S P -7 9 b W A S P -1 0 3 b W A S P -1 2 1 b  Vertical Thermal gradient (K) Figure A3. Correlations at the day-side between the H2O (top) and e − (bottom) abundances and the mean retrieved temperature weighted by the contribution function. The colours indicate the retrieved thermal gradient, defined as the temperature differences between the upper and lower quantiles of the contribution function. The shaded green region is the predicted abundances from equilibrium chemistry at solar metallity and C/O ratio between 1 bar and 0.01 bar. In orange, the metallicity is increased to 100 times solar. In red, the metallicity is decreased to 0.1 times solar. In blue, the metallicity is solar and the C/O ratio is increased to 1. Note that in the e − plot (bottom), the Solar and C/O=1 cases overlap. Vertical Thermal gradient (K) Figure A4. Correlations at the day-side between the VO (top) and TiO (bottom) abundances and the mean retrieved temperature weighted by the contribution function. The colours indicate the retrieved thermal gradient, defined as the temperature differences between the upper and lower quantiles of the contribution function. The shaded green region is the predicted abundances from equilibrium chemistry at solar metallity and C/O ratio between 1 bar and 0.01 bar. In orange, the metallicity is increased to 100 times solar. In red, the metallicity is decreased to 0.1 times solar. In blue, the metallicity is solar and the C/O ratio is increased to 1. Vertical Thermal gradient (K) Figure A5. Correlations at the day-side between the FeH (top) and CO (bottom) abundances and the mean retrieved temperature weighted by the contribution function. The colours indicate the retrieved thermal gradient, defined as the temperature differences between the upper and lower quantiles of the contribution function. The shaded green region is the predicted abundances from equilibrium chemistry at solar metallity and C/O ratio between 1 bar and 0.01 bar. In orange, the metallicity is increased to 100 times solar. In red, the metallicity is decreased to 0.1 times solar. In blue, the metallicity is solar and the C/O ratio is increased to 1. Note that in the FeH plot (top), the Solar and C/O=1 cases overlap. Vertical Thermal gradient (K) Figure A6. Correlations at the day-side between the CO2 (top) and CH4 (bottom) abundances and the mean retrieved temperature weighted by the contribution function. The colours indicate the retrieved thermal gradient, defined as the temperature differences between the upper and lower quantiles of the contribution function. The shaded green region is the predicted abundances from equilibrium chemistry at solar metallity and C/O ratio between 1 bar and 0.01 bar. In orange, the metallicity is increased to 100 times solar. In red, the metallicity is decreased to 0.1 times solar. In blue, the metallicity is solar and the C/O ratio is increased to 1.  Figure A8. Metallicity (O/H) and C/O ratios recovered from our free chemistry retrievals. In colors, we highlight the error on the retrieved abundance of water σ log(H 2 O) , since only the cases where water is recovered are reliable (purple to black). C/O remains very difficult to constrain since our observations lack sensitivity to carbon bearing species. In the inset, we show the metallicity derived from free chemistry for the runs with σ log(H 2 O) < 1 versus the metallicity recovered in the equilibrium chemistry cases.

B.1. Description of the observations and their reduction
For this study, we considered the planets that have been observed in eclipse using the Hubble Space Telescope (HST) WFC3 camera with the Grism G141. This constitutes a sample of 25 planets. For all planets except HAT-P-70 b, data from the Spitzer Space Telescope is also available for at least the 3.6 µm and the 4.5 µm IRAC channels. When available, we also re-analysed the complementary HST WFC3 data from transit observations, which is the case for 17 planets of our sample. We describe below our reduction method for the HST observations and the sources for the Spitzer data.

B.1.1. Hubble Space Telescope data
For all planets, except for Kepler-13 A b and WASP-33 b, we downloaded the publicly available data from the Mikulski 3 Archive (MAST). We note that data from HST also exist for the planet CoRoT-2 b, but due to the very shallow eclipse for this planet (Wilkins et al. 2014), we were not able to reliable reduced the observations and chose to not include this planet in the study. The publicly available data consist of a series of raw detector images.
We carried out the analysis of all HST WFC3 data using Iraclis, our highly-specialised software for processing WFC3 spatially scanned spectroscopic images (Tsiaras et al. 2016b) which has been used in a number of studies (Tsiaras et al. 2016c;Damiano et al. 2017;Tsiaras et al. 2018;Tsiaras et al. 2019;Anisman et al. 2020;Changeat et al. 2020b;Skaf et al. 2020;Guilluy et al. 2021;Mugnai et al. 2021;Yip et al. 2021;Libby-Roberts et al. 2021). For this study and except for two planets (WASP-33 b and Kepler-13A b), we acquired literature spectral data only if the reduction had been performed with the Iraclis pipeline to ensure homogeneity between planets. In particular, Tsiaras et al. (2018), later referred to as T18 (Tsiaras et al. 2018), already reduced a number of the datasets for the transits we consider here with Iraclis. Else, we performed a reduction process that included the following steps: zero-read subtraction, reference-pixels correction, non-linearity correction, dark current subtraction, gain conversion, sky background subtraction, calibration, flat-field correction and bad-pixels/cosmic-rays correction. Then we extracted the white (1.088-1.68 µm) and the spectral light curves from the reduced images, taking into account the geometric distortions caused by the tilted detector of the WFC3 infrared channel.
We fitted the light curves using our light curve modelling package PyLightcurve (Tsiaras et al. 2016a) with the parameters from Table B1. During our fitting of the white light curve, the planet-to-star flux ratio and the mid-eclipse time were the only free parameters, along with a model for the systematics (Kreidberg et al. 2014c;Tsiaras et al. 2016b). It is common for WFC3 exoplanet observations to be affected by two kinds of time-dependent systematics: the long-term and short-term 'ramps'. The first affects each HST visit and is modelled by a linear function, while the second affects each HST orbit and has an exponential behaviour. The formula we used for the white light curve systematics (R w ) was the following: where t is time, n scan w is a normalization factor, T 0 is the mid-eclipse time, t o is the time when each HST orbit starts, r a is the slope of a linear systematic trend along each HST visit and (r b1 , r b2 ) are the coefficients of an exponential systematic trend along each HST orbit. The normalization factor we used (n scan w ) was changed to n f or w for upward scanning directions (forward scanning) and to n rev w ) for downward scanning directions (reverse scanning). The reason for using separate normalization factors is the slightly different effective exposure time due to the known upstream/downstream effect (McCullough & MacKenty 2012).
We fitted the white light curves using the formulae above and the uncertainties per pixel, as propagated through the data reduction process. However, it is common in HST/WFC3 data to have additional scatter that cannot be explained by the ramp model. For this reason, we scaled up the uncertainties in the individual data points, for their median to match the standard deviation of the residuals, and repeated the fitting (Tsiaras et al. 2018). The white light-curve fits obtained for the eclipse spectra are shown in Figure Table B1. List of the orbital parameters used in this paper and the associated literature references. Rs is the star radius in Sun radius. Ts is the star temperature in Kelvin. g is the stellar surface gravity in cm/s 2 . [M/H] is the stellar metallicity in dex, P is the orbital period of the planet in days. a is the semi-major axis. e is the eccentricity. i is the inclination in degrees. t mid is the mid-transit time is days. Mp is the planetary mass in masses of Jupiter. For HAT-P-70 b only a upper limit on the planetary mass is available. Next, we fitted the spectral light curves with a transit model (with the planet-to-star flux ratio being the only free parameter) along with a model for the systematics (R λ ) that included the white light curve (divide-white method (Kreidberg et al. 2014c)) and a wavelength-dependent, visit-long slope (Tsiaras et al. 2016b).
where χ λ is the slope of a wavelength-dependent linear systematic trend along each HST visit, LC w is the white light curve and M w is the best-fit model for the white light curve. Again, the normalisation factor we used (n scan λ ) was changed to (n f or λ ) for upward scanning directions (forward scanning) and to (n f or λ ) for downward scanning directions (reverse scanning). Also, in the same way as for the white light curves, we performed an initial fit using the pipeline uncertainties and then refitted while scaling these uncertainties up, for their median to match the standard deviation of the residuals.
We note that several of these datasets were acquired using the staring mode which is usually less efficient than the spatial scanning technique. Furthermore, the signal to noise ratio achieved on the eclipse varied between planets. We used the SNR on the white light curve to dictate the resolution of the spectral fitting. For those with an SNR >7, a high resolution spectrum was extracted (25 bins) while a lower resolution binning was chosen for all others (18 bins). The recovered spectral light-curve for the planets reduced with Iraclis are shown in Figure Set B2.
For planets with multiple visits, we correct for offsets by subtracting each spectrum by the corresponding white light curve depth, (F p /F ) 2 w,v , and adding the weighted average eclipse depth of all white light-curves, (F p /F ) 2 w . Finally, we compute the weighted average from all the eclipse observations which we use for all subsequent analysis. Due to possible variabilities in the star, the planet or the instrument systematics, the measured white eclipse depths could change between the different visits, so we choose not perform a joint-fit of these datasets as done in many previous studies (Cartier et al. 2017;Wakeford et al. 2018;Arcangeli et al. 2018). However, except for HD 209458 b, when multiple visits are available, the individual eclipse depths are all consistent within 1-sigma, meaning that a joint-fit would be equivalent. For HD 209458 b, we show the extracted spectra in Figure B3 along with the final spectrum obtained by averaging all five visits. We note that the spectral features are consistent in the five observations. While most of our spectra are consistent with the literature, HD 209458 b differs significantly from the spectrum derived in Line et al. (2016), as shown in Figure B3. The shape of the spectra is highly similar, but we observe a ∼160 ppm vertical offset. We identified that the difference comes from the assumption of a quadratic systematic for the long-term ramp in Line et al. (2016), while we here assume a linear behaviour. We note that Line et al. (2016) only managed to fit four of the five eclipse observations while we fitted all five. Additionally, when fitting these observations with a quadratic trend we found the white light curve depth, and associated errors, were far more unstable than the linear fits. The spectral shape recovered with both trends was very similar. HD 209458 b is the planet for which the most significant difference to the literature can be observed. Other spectra from the literature are compared to our Iraclis pipeline reductions in the Figure Set B4. Overall, we find that except for HD 209458 b and for TrES-3 b, we do not find significant differences in spectral shape and observe only minor offsets.
We also considered the transmission spectra, which are available for 17 planets. These were reduced using the same process as for the emission data.
For two systems, WASP-12 and WASP-103, a close stellar companion contaminated the data from HST WFC3. For exoplanet spectroscopy, this third light modifies the transit/eclipse depth. To account for this, we used the freely available WFC3 simulator Wayne 4 .
Wayne is capable of producing grism spectroscopic frames, both in staring and in spatial scanning modes (Varley et al. 2017). We utilised Wayne to model the contribution of each companion star to the spectral data obtained. We created simulated detector images of both the main and companion star, using these to extract the flux contribution in each spectral bin of each star. The correction to the spectra is then applied as a wavelength dependent dilution factor which is derived as a ratio of extracted flux between the stars. Such an approach has previously been used on WFC3 data, e.g.    Table B2. Summary of the observations considered in this paper. HST: Hubble Space Telescope; Spz: Spitzer Space Telescope. We used observations from the Grisms G102 and G141 aboard HST and photometric channels from 3.6 to 24 µm aboard Spitzer. For the Spitzer data, we used the data from the literature directly. For all the planets in our sample, we find that at least the 3.6 µm and the 4.5 µm channels from the Infrared Array Camera (IRAC) are available. Adding Spitzer data significantly increases the wavelength coverage of our observations, but can also lead to biases when the observations are not compatible. This well-known issue is discussed further in Complementary Text. In order to be as consistent as possible, we prioritised the inclusion of the Spitzer data from the population study of Garhart et al. (2020), later referred as G20. A few planets possess additional channels from legacy Spitzer observations with the Infrared Array Camera (IRAC: 5.8 µm and 8 µm channels), the InfraRed Spectrograph (IRS: 16 µm) and the Multiband Imaging Photometer for Spitzer (MIPS: 24 µm). In this case, we add the most complete reduction, which allows to greatly increase the information content for these planets. For HD-209458 b, we took the Spitzer data from Diamond-Lowe et al. (2014). We noted, however, that the 16 µm point used, in fact, refers to the Spitzer IRS low spectral resolution observations (Swain et al. 2008a), which span the wavelength range from 7.46 µm to 15.25 µm. We therefore chose to include the derived broad band eclipse depth from the original study. In Spitzer, photometric channels are large and can cover entire or multiple broadband molecular features. The shape of the spectral response function is therefore important and has to be accounted for during binning of the forward model. Spitzer spectral responses functions for the IRAC, IRS and MIPS instruments can be found at the NASA/IPAC Infrared Science Archive 5 .
The full list of observations that are considered in the retrievals is provided in Table B2. We detail the spectra recovered from this reduction step in Tables B3 and B4 for, respectively, eclipse and transit observations. The full tables are available in machine-readable format. In those tables, we also include for convenience the Spitzer datasets  Table B3. List of eclipse spectra used in this population study. λ refers to the central wavelength of the bin and ∆λ is the bin width. For Spitzer, the bin widths are irrelevant since we consider the spectral response of the channel.
Note- Table B3 is published in its entirety in the electronic edition of the Astrophysical Journal. A portion is shown here.  Table B4. List of transit spectra used in this population study. λ refers to the central wavelength of the bin and ∆λ is the bin width.
Note- Table B4 is published in its entirety in the electronic edition of the Astrophysical Journal. A portion is shown here.
that were obtained from the literature. These spectra constitute the inputs of our retrieval analysis, further described in the next section.

B.2. Standardised retrievals with Alfnoor
To analyse the spectra obtained by Iraclis, we perform atmospheric retrievals using the Alfnoor tool (Changeat et al. 2020a). Alfnoor extends the capabilities of the Bayesian retrieval suite TauREx3 (Al-Refaie et al. 2021b,a) and automate retrievals to large exoplanet populations. It was first built for simulations in the context of the ESA space mission Ariel (Tinetti et al. 2018(Tinetti et al. , 2021 to enable retrieval studies of the mission's entire target list (Edwards et al. 2019a), but it can also perform standardised retrievals from any real data observation. We use Alfnoor to extract the information content of our planetary atmospheres (chemical composition and temperature structure) separately for the transit and eclipse scenarios. For spectroscopic data, such as what is presented here, Bayesian retrievals are the currently most adopted analysis technique to extract un-biased information from atmospheric spectra. Another method that has been used to investigate population trends in a few studies (e.g. Crossfield & Kreidberg 2017;Mansfield et al. 2021) consists of measuring the signal amplitude of the molecules of interest, for instance water. The technique compares the observed flux in the HST water band, which is defined as between 1.35 µm and 1.48 µm for instance in Mansfield et al. (2021), with a reference flux taken outside the band to evaluate deviations from the expected signal. Such a method, however, does not rely on a full exploration of the parameter space. More importantly, the method requires the definition of signal and reference bands, both of which can contain additional molecular features as shown in our work. Here, we prefer to extract trends in our population by atmospheric retrievals with a Nested Sampling optimizer and by testing a variety of atmospheric scenarios.

B.2.1. Emission and Transmission forward models from TauREx3
For both transit and eclipse scenarios, the atmosphere is modeled assuming 1D layers (default 100 layers). In the transit case (Waldmann et al. 2015b), the total transit depth at wavelength λ is given by: where R p is the planet radius and R s is the parent star radius. a λ is the wavelength contribution from the atmosphere (transit depth), which takes the form: where z max is the altitude at the top of the atmosphere and τ λ (z) is the wavelength dependant optical depth. It is evaluated by: with τ λ,i is the optical depth of each absorber i. In eclipse (Waldmann et al. 2015a), the emission from each layer is integrated to produce the final spectrum. The wavelength dependant intensity at the top of the atmosphere from a viewing angle θ is: where µ = cos(θ), B λ (T ) is the plank function at a given temperature T, T s denotes the temperature at maximum atmospheric pressure and τ s is the total optical depth from the planetary surface to the top of the atmosphere. Then the flux is integrated for the cosine viewing angle µ using an N point Gauss-Legendre quadrature scheme: where w i and x i are our weights and abscissas respectively. The final emission spectrum is expressed as: where I s is the specific intensity from the star. In this work modelled using PHOENIX spectra (Allard et al. 2012).

B.2.2. Opacity sources
Molecular opacity sources are included in both transmission and emission using their cross sections ζ i,λ . Their contribution to the optical depth is then given by: where χ i is the column density of the species i and ρ, is the number density of the atmosphere. The contributions are integrated along the line of sight parametrised by dz . For planets with temperatures higher than 2000K, the day-side is expected to be cloud free. For planets in lower temperature regimes, we tested the presence of clouds on the day-side using two parametric models. We first tested a fully opaque grey cloud layer with a single free parameter, the cloud deck top pressure. Such model is often used when analyzing HST data due to the small wavelength coverage in near-infrared. Since we also include the infrared bands from Spitzer in the eclipse analyses, we also tested a more realistic Mie model from Lee et al. (2013). The retrieved parameters for this model are cloud top pressure, cloud particle radius, and cloud mixing ratio. Comparing the Bayesian evidence of both of the cloudy retrievals with the cloud-free case, we did not find evidence in flavor of clouds in any of the eclipse spectra. Therefore, for the rest of the paper, we did not include clouds for eclipse retrievals. We note, however, that tested cloud models only simulate extinction processes and do not account for complex scattering. In transit, since the observations are more sensitive to clouds, we included the fully opaque grey cloud model in all retrievals.

B.2.3. Retrieval model setups
In order to build our comparative study, we perform the same standardised retrievals for all the planets using Alfnoor. All the planets in our study possess eclipse observations. The planet atmospheres are modelled using plane-parallel assumption and composed of 100 layers spaced in log-pressures from 10 6 Pa to 10 −5 Pa. The stellar parameters and planetary mass are always fixed to literature values as dedicated studies provide better constraints than what can be recovered from our data (Changeat et al. 2020c).
In eclipse, due to the stronger degeneracies with the thermal profile, we fix the planetary radius to the literature values. We parametrise the temperature profile using a N-point point free profile containing three nodes. This is a heuristic profile that linearly interpolates between N freely moving points. We perform a total of eight retrievals using four different models: • Reduced: This conservative model uses free chemistry with constant with altitude mixing ratios. It includes the molecular opacities from H 2 O, CH 4 , CO and CO 2 for a total of 9 free parameters.  Figure B5. List of the free parameters and their uniform priors in the retrievals.
• Full: This model uses free chemistry with constant with altitude mixing ratios. In addition to the molecules of the reduced run, it includes the optical absorbers TiO, VO, FeH and H − for a total of 13 free parameters.
• Equilibrium (Eq): This model uses equilibrium chemistry from the GGChem model (Woitke et al. 2018). All the species available in GGChem were used to compute the chemistry and all the available opacity sources (see Opacity section) were included. To explore the thermal profile, we use the same N-Point parametrization as for the free cases, so the total number of free parameters is 7.
• Blackbody: For reference we perform a retrieval with no molecules and a simple isothermal profile (1 free parameter). This is used as a reference to compare the Bayesian evidence of the other models.
For each of those different models, we run an HST only case and an HST+Spitzer scenario. When Spitzer is included, we use the calibrated Spitzer instrument responses (Reach et al. 2005) to account for the wavelength dependant behaviour of the detector.
When available, we complement our eclipse retrievals with an analysis of the HST transit spectrum. In this case, we also retrieve the planet radius R p and a Grey opaque cloud deck as these observations are more sensitive to those properties. Since in transmission the impact of the temperature is mostly on the scale height, we restrain our model to an isothermal temperature structure (Changeat et al. 2019;Rocchetto et al. 2016). Here we only analyse the HST-WFC3 spectra and do not consider the additional Spitzer data. We run two separated retrievals: • Full: This model uses free chemistry with constant with altitude mixing ratios. We include all the absorbers: H 2 O, CH 4 , CO, CO 2 , TiO, VO, FeH and H − , for a total of 11 free parameters.
• Equilibrium (Eq): This model uses equilibrium chemistry from the GGChem model (Woitke et al. 2018) and the same isothermal temperature profile. The number of free parameters is 5.
• Featureless fit: For reference, we performed a transit fit with no molecule, fitting only for a flat line. The number of free parameters is 1.

B.2.4. Exploration of the parameter space
All the free parameters considered in this study are explored using uniform, non-informative priors. These are listed in more details in Table B5. For the chemistry, since all the planets considered are large hot Jupiters, we assume a primary envelope and therefore limit the retrieved maximum bound of the considered active molecules to 0.01. We imposed stricter priors on the optical absorbers, limiting their abundance to a maximum of 10 −4 as their abundance is expected to be much less (Woitke et al. 2018). For the lower bounds, we select 10 −12 as such abundance is low enough to not leave any spectral features at the resolution and signal-to-noise considered. The sampling is done using the nested sampling algorithm MultiNest (Feroz et al. 2009) with 750 live points and a log-likelihood tolerance of 0.5. This ensures an optimal free exploration of the parameter space.
Since we use the nested sampling algorithm MultiNest, the computation of the Bayesian Evidence for each model, here denoted E, is automatic. Then, the difference in ln(E) between two models M1 and M2 can be used for model  (1995). The Bayes Factor B is defined as the evidence ratios of the models M1 and M2. In particular ln(B) = ∆ln(E). Note that we use ln(E), the natural logarithm of the Evidence, which is the standard MultiNest outcome and is often noted log(E) in the literature.
selection and to compare the ability of the two models to explain the observed spectra (Jeffreys 1961;Kass & Raftery 1995;Feroz et al. 2009). The original Jeffreys' scale that we employ in this study is summarised in Table B5.
C. SUPPLEMENTARY TEXT

C.1. Model choices and degeneracies
Atmospheric retrieval techniques were historically developed for Earth and Solar System applications, where in-situ and high quality measurements provide a wealth of information. In this context, retrievals are used to infer probability distributions on a set of parameters for which priors knowledge is already available and tightly constraining. On the opposite, exoplanet retrievals must deal with low information content data sets and physical systems that are poorly understood. In this context, model assumptions are crucial (Min et al. 2020). These assumptions are user defined and should be selected carefully, as a retrieval will always provide a model dependent solution (Rocchetto et al. 2016;Changeat et al. 2019;Changeat et al. 2020a). In general, free assumptions attempt to heuristically describe properties of exoplanets, such as thermal or chemical profiles, using as few assumptions as possible. On the opposite, self-consistent models describe the given property using physically driven assumptions, for example chemical equilibrium or radiative equilibrium. In a low prior knowledge case, free retrievals should be favored until the understanding of the system is high enough to justify the assumptions of self-consistent models. Another issue of self-consistent modeling in retrievals is that the addition of more physical descriptions comes at the cost of computational resources, which often limits the level of details that can be included. In this work, we chose to explore both options, using free and equilibrium chemistry retrievals. Often, we find that the free approach leads to a better Bayesian evidence, thus demonstrating that the added flexibility in the free models provides a better explanation of the datasets considered here.
In this work, we focus on eclipse spectra, which as compared to transit observations, provide more constraints on the thermal profile (Burrows 2014;Rocchetto et al. 2016). Due to the narrow wavelength range of HST, only one broadband spectral feature is usually observed, which can be difficult to interpret without further constraints. In the HST only case, the retrievals can easily be misguided, confusing for instance whether a single feature should be seen in emission or absorption. This is the case for HD 209458 b (see corresponding section), where the HST only retrievals converge to two different solutions depending on the opacities that are considered. When optical absorbers are included, the retrieval interprets the spectrum with an emission feature, while it is simply fit with absorption of H 2 O in the reduced run. Adding the Spitzer data, here, provides confidence that the second scenario is more likely, but such degeneracies are often difficult to disentangle. In this work, we performed retrievals using a reduced list of opacities and a full list, as well as HST only and HST+Spitzer scenarios. In general, we find that a larger wavelength coverage is often required to extract information regarding the chemistry and thermal profile of the observed atmospheres in eclipse and we therefore focus on the HST+Spitzer runs for the main part of this study.

C.2. Biases arising from the data reduction
In previous works, data reduction pipelines have led to different answers when used on the same raw data. When considering HST, different pipelines usually recover the same spectral shape, but differences can appear as flat offsets (Changeat et al. 2020b, e.g.). See also Figure Set B5, which compares our spectra with the ones obtained in the literature.
The observed flat offsets are usually attributed to differences in the background subtraction steps or different choices to model the HST ramps. A particularly striking example of such offset can be seen in our reduction of the HD 209458 b spectrum, which differs by about 160 ppm with the one obtained by Line et al. (2016). Such large offsets are difficult to explain by differences in the background substraction alone. However, performing the reduction of the spectra using a quadratic trend for the long-term ramp, we obtained a similar flux level to Line et al. (2016). We also note that Line et al. (2016) fitted all the visits with a single transit model, which could have led to different results as their reduction did not rely on a normalization of the white light curve depths. While our pipeline does not offer this option at the moment, future investigations could assess the impact of performing joint or separate fits. In order to investigate further this particular case, we performed comparative retrievals on the Line et al. (2016) spectrum as well as our own spectrum offset by -160 ppm. In Figure C1 we show the results of this exercise.
In the same figure, we also show self-consistent estimates of the HD 209458 b emission using a two-stream model available in TauREx. While not suitable for direct data analyses studies, self-consistent forward models provide important benchmark to compare retrieval outcomes. This forward model computes the thermal profile and chemistry in radiative and chemical equilibrum from the bulk and orbital parameters of the system, which are here taken from the literature. We show two models, one with solar metallicity and C/O ratio and one with a metallicity 0.2 and a C/O ratio of 1.0, which roughly matches the derived values of our full HST+Spitzer run. Those forward models predicts a surface emission that might be compatible with the reduced observation we obtained from Iraclis. This exercise is by no mean proving that our recovered spectrum is correct, but this highlight systematic biases that can arise from different assumptions in HST reduction pipeline.
These offsets can even be more difficult to handle when combining observations from different instruments in transit (Yip et al. 2020;Changeat et al. 2020b;Yip et al. 2021;Saba et al. 2021). In eclipse or phase-curve , similar incompatibilities can occur depending on the different reduction techniques or input parameters. Similarly, if the observations are not carried out simultaneously, one encounters the risk of contamination from stellar or even planetary variabilities. Indeed, dynamic simulations predict that exoplanet atmospheres are subject to large variabilities, with time-dependent storms and modons (Cho et al. 2003;Cho et al. 2021). In this work, we do not investigate those phenomena further, and for simplicity, we assume that the obtained HST and Spitzer observations are compatible. We, however, performed a few sensitivity tests.
In order to verify the stability of our results to potential offset biases in the Spitzer data, we performed complementary retrievals on the whole population. We applied the full scenario on the HST+Spitzer data, modifying the observed Spitzer data. In our first test, we doubled the error of the Spitzer data for all planets, which is equivalent to adding  another source of unknown Gaussian noise. In the two supplementary tests, we tested potential biases in the Spitzer data, by offsetting the eclipse depths of the Spitzer data by ± 100 ppm. In those three tests, we occasionally found different solutions for individual planets but overall the trends reported using the un-modified Spitzer sets remained unchanged. This is shown in Figure C2, which displays the water-temperature map for our population in the case where the Spitzer noise has been doubled, and the case where Spitzer has been offset by +100 ppm.
C.3. On constraining metallicity and C/O ratio from free retrievals.
In our free runs, we also compute the metallicity The results of these calculations on our population are shown in Figure A8. We find that the estimates of metallicity and C/O ratios from free retrievals heavily rely on the detected molecules (see behaviour of the C/O ratio in the individual planet analyses) and come with many simplifying assumptions. Often, when a single molecule is detected in the retrieval, the derived parameter will be heavily biased towards the elemental ratio of the detected molecule. For example, detecting only CO 2 will lead to a derived C/O ratio of 0.5, which is spurious.
By comparing the derived metallicity to the one obtained by the equilibrium runs, we however believe this parameter is stable in the case where H 2 O is clearly detected. This is shown in the top left panel of Figure A8, where we show that the retrievals detecting H 2 O to high accuracies lead to similar estimates of the metallicity in free and equilibrium runs. We attribute this to the good capabilities of HST to detect this particular molecule, and our better understanding of the class of planets between 1000K and 2000K, which are believed to more closely following equilibrium chemistry.  Table D1 summarizes the recovered parameters and provides the Bayesian Evidence of each model. For molecules, non-detections are characterized by a retrieved log(VMR) below -8 for the mean, and if the mean minus the sum of the 1σ uncertainties (on both sides) is lower than -12. In this case, only the upper limit is provided.  Temperature (K)    Table D1 continued   HAT-P-41 b e log(H2O) < −6.9 < −6.9 < −4.8  Table D1 continued    Table D1 continued  log(H2O) < −6.9 < −6.1   Table D1 continued     Table D1 continued          Table D1 continued     Table D1 continued  Note- Table D1 is published in its entirety in the electronic edition of the Astrophysical Journal.

D.1. Individual analysis of CoRoT-1 b
The first transit of CoRoT-1 b was observed in 2008 ( Barge et al. 2008), unveiling a 1.5 R J inflated planet. Later observations from the ground indicated a significant emission signal corresponding to a temperature of approximately 2300 K and potentially a poor heat redistribution between the day-and night-sides (Alonso et al. 2009;Gillon et al. 2009). The planet transit from HST WFC3 (PN: 12181, PI: Deming 2010) was previously studied but no evidence for molecular absorption was found due to the high uncertainties on the data (Ranjan et al. 2014). The same proposal included eclipse observations, but no analysis of these has since been published. However, a follow-up study in 2010 analyzed Spitzer eclipse observations (Deming et al. 2010) and found potential evidence for high-altitude absorbers or an isothermal region in the atmosphere. We analyse the transmission and emission HST visits using our standardized pipeline and performed the retrieval analyses (see Materials and Methods section). For the eclipse, we considered the Spitzer observations reduced in Deming et al. (2010), which are taken without modifications.
On the day-side, the HST and HST+Spitzer are consistent with a thermal inversion. The posterior distributions indicate spectral signatures for H 2 O, VO and H − opacities, which is consistent with the findings in Deming et al. (2010). The retrieved temperatures of about 2000 K could explain the presence of the detected absorbers. For the free runs including the full list of opacities, the chemistry is consistent with a solar metallicity and a sub-solar C/O ratio. We note that, since no C-bearing species are detected, the C/O ratio from the free runs is likely biased. When equilibrium chemistry is used, a thermal inversion is also recovered, with different thermal profiles depending whether Spitzer is included or not. The equilibrium retrievals tend to converge towards an atmosphere with solar or enriched metallicity and a C/O ratio of about 1. Our analysis provides decisive evidence in favour of the free model that includes H 2 O and the optical absorbers VO and H − .
For the transit observation, the spectrum shows a downward slope towards the red end of the spectrum. In the free run, it is best fit with a cloudy atmosphere and a potential contribution from VO. This result is consistent with the detected absorption of VO on the day-side and the recovered temperature of around 1800 K. In this free scenario, the atmosphere might be consistent with a solar metallicity and C/O ratio, but we note that the uncertainties on those derived parameters are large, also allowing for sub-solar to super-solar C/O ratios. The equilibrium run leads to similar estimates of the metallicity and the C/O ratio, with better con-straints, but the recovered temperature is un-physical, with a mean higher than at the day-side.

D.2. Individual analysis of HAT-P-2 b
HAT-P-2 b was discovered in 2007 by Bakos et al. (2007). It is a massive hot Jupiter (9.1 M jup ) that orbits its host star in an highly eccentric orbit (e = 0.52) in about 5.6 days. Due to its large density, the planet is believed to require the presence of a large core. This large mass, combined with the highly eccentric orbit, raises many questions regarding the physics of this planet and its formation. For instance, along the entire orbit, the planet's equilibrium temperature varies from 1240K to 2150K Bakos et al. (2007. Studying the Rossiter-McLaughlin effect, Winn et al. (2007);Loeillet et al. (2008) found that stellar spin axis and orbital axis of the planet should be aligned, thus implying that the planet did not evolved through scattering or Kozai migration.
While being a very interesting planet, the atmosphere of HAT-P-2 b was not studied with many instruments. A phase-curve observation with Spitzer at 3.6 µm, 4.5 µm, 5.6 µm and 8 µm was presented in Lewis et al. (2013), highlighting a very complex atmosphere due to the particular orbital configuration of this planet. The study also suggested the planet might experience a temporary day side thermal inversion near periapse. In a follow-up work, Lewis et al. (2014) performed a complementary analysis with GCM models to evaluate the impact of the eccentricity on the chemistry and the thermal structure of this planet, highlighting that dis-equilibrium processes on this planet might be important.
Recently, a partial phase-curve was acquired with HST using the Grism G141 (PN: 16194, PI: Desert et al. 2020). We fitted the publicly available data using our standardised Iraclis pipeline. This data was complemented with the Spitzer observations, taken as-is from Lewis et al. (2013).
We performed our the retrieval of this planet with Alfnoor. All our retrieval models provide the same interpretation of the data and the HAT-P-2 b eclipse spectrum does not show clear signatures for any molecule. The atmosphere is best interpreted by an isothermal profile, with eventually some marginal water feature. However, comparisons of the Bayesian evidence indicate that the detection is tentative at best. The spectrum is well explained by a simple blackbody fit and we are not able to confirm the thermally inverted profile predicted by Lewis et al. (2013). Due to the fact that the HST and Spitzer observations were not carried simultaneously, we note that the atmospheric conditions might have changed between the different observations, in par-ticular considering the additional variability introduced by the highly eccentric orbit.

D.3. Individual analysis of HAT-P-7 b
HAT-P-7 b is an inflated hot Jupiter of 1.4 R J (Pál et al. 2008), which was first studied in emission during the commissioning program of Kepler when the satellite detected the eclipse as part of an optical phase curve (Borucki et al. 2009). These measurements indicated that HAT-P-7 b could have a day-side temperature of around 2650 K which confirmed predictions from (Pál et al. 2008;Fortney et al. 2008). This optical eclipse measurement was combined with Spitzer photometry over 3.5 -8 µm to infer the presence of a thermal inversion (Christiansen et al. 2010), suggested by the high flux ratio in the 4.5 µm channel of Spitzer compared to the 3.6 µm channel. In their paper, chemical equilibrium models associated these emission features with CO, H 2 O and CH 4 . A thermal inversion was also reported to provide the best fit to this data by the atmospheric models of Spiegel & Burrows (2010); Madhusudhan & Seager (2010) but all three studies noted that models without a thermal inversion could also well explain the data though only with an extremely high abundance of CH 4 . Further Kepler phase curves identified an offset in the day-side hot-spot (Esteves et al. 2013;Esteves et al. 2015) as well as changes in its location (Armstrong et al. 2016), highlighting the complex dynamics of hot Jupiter atmospheres. However, while Spitzer phase curves at 3.5 & 4.5 µm were also best fitted with a thermal inversion on the day-side and relatively inefficient day-night recirculation, Wong et al. (2016b) did not find evidence of a hot-spot offset. Two eclipses were then obtained using HST WFC3 G141 which, when combined with previous observations, were found to be best fit with a thermal inversion due to optical absorbers but at a low significance when compared to a baseline black-body fit (Mansfield et al. 2018).
We fitted the two scanning mode eclipse observations of HAT-P-7 b (PN: 14792, PI: Bean et al. 2016). We note that three staring mode eclipse observations were also taken (PN: 12181, PI: Deming 2010) but we ignore these given the higher precision of the scanning mode observations. Additionally, a staring mode transmission spectrum had been taken (PN: 12181, PI: Deming 2010) but there was no post-egress orbit and, due to significant systematics, our fit of the data presents large uncertainties. Models of HAT-P-7 b suggest the terminator region should have patchy clouds with water being well mixed throughout the atmosphere (Helling et al. 2019).
This planet, which was extensively observed with Spitzer, was not included in the population study from G20. We therefore considered the observations from Christiansen et al. (2010), which cover the wavelength range from 3.6 to 8 µm. In their paper, on top of the standard MCMC technique, they used the 'rosary bead' residual permutation method (Winn et al. 2008) to quantify potential remaining systematic errors. Their results highlighted differences in the 8 µm point, so we use the 8 µm channel from the rosary technique, as suggested in Christiansen et al. (2010).
We performed our standardised retrieval study of this planet. Analysing the full HST-only run does not clearly favour a molecule and leads to an isothermal atmosphere. When Spitzer is included, however, additional constraints can be extracted from the 4 Spitzer photometric channels. The retrievals now include some evidence for sub-solar abundance of H 2 O, high abundance of CO 2 and the presence of FeH. The associated temperature profile to best fit the HST+Spitzer data-set contains a localised thermal inversion. While the recovered temperature reaches 2500 K, we do not find evidence for H − opacities. Comparing the Bayesian evidence of the free models, it is not possible to validate the detection of molecular species from the HST only data-set, thus demonstrating that the HST spectrum is in this case, uninformative and consistent with a blackbody. It is only when Spitzer is included that we find strong evidence in favour of the models with molecular opacities. The derived metallicity for the free models is only constrained when HST+Spitzer data are considered, favouring a super-solar metallicity. Due to the detection of large abundances of CO 2 in the model that includes Spitzer data, the derived C/O ratio takes a value close to 0.5. When testing the equilibrium retrievals, the thermal profile is decreasing with altitude if Spitzer is included or possesses a localised thermal inversion if only HST is considered. In both runs, the additional constraints of this chemical model do not allow to clarify the metallicity of this atmosphere, but the associated C/O ratio is restricted to high values. Looking at the spectra, the solution obtained by the equilibrium chemistry retrievals obtain a much worse fit of the observed Spitzer data, which is also confirmed by the lower Bayesian evidence.
At the terminator, we do not recover the presence of any molecular species, the spectrum being consistent with flat. The planet could therefore have high-altitude clouds at the terminator.

D.4. Individual analysis of HAT-P-32 b
The planet HAT-P-32 b was first reported in 2011 (Hartman et al. 2011) and presents a particularly large radius: 1.98 R J (Wang et al. 2019b). The planet eclipse was then observed with HST WFC3 (Nikolov et al. 2018) and Spitzer (Zhao et al. 2014). Both analyses suggested that thermal inversions could be present in this planet. The terminator region was also observed with HST WFC3 and it is consistent with a significant water feature at 1.4 µm (Damiano et al. 2017;Tsiaras et al. 2018). This data was analysed together with HST STIS and Spitzer IRAC data, which showed a thick cloud layer and a super-solar metallicity (Alam et al. 2020).
Since the transmission spectrum (PN: 14260, PI: Deming et al. 2015) was part of the T18 study and was already reduced with Iraclis, we took the spectrum "as is", directly from this paper. For the eclipse (PN: 14767, PI: Sing et al. 2016b), we re-analysed the raw images using Iraclis following our standardised reduction technique. We include the Spitzer data from Zhao et al. (2014) in our spectral retrieval analysis.
For the eclipse, we show the fitted spectra, the temperature profiles and the posteriors in Figure Set D1. In the full runs, the recovered temperature profile is slightly inverted and is similar between the HST and HST+Spitzer runs. In addition, both runs are consistent with H − , which could be causing the observed thermal inversion. We note, however, that the recovered temperatures for this atmosphere are below the predictions for H 2 dissociation in Parmentier et al. (2018). Such results are surprising but given the small differences in Bayesian evidence between with the reduced model, we believe this detection to be marginal. The addition of the Spitzer photometric points lead to some hints for VO and CH 4 , but these detections also remain weak, with large tails in the posterior distributions. When comparing the full runs to the reduced runs, the inclusion of the Spitzer points favours the models with optical absorbers, but this is not the case when only HST is considered. For the free runs, we find that the full runs prefer a subsolar metallicity case, which contrast with the results from equilibrium chemistry retrievals. In terms of C/O ratio, solar values are also allowed with large uncertainties for the full runs, while equilibrium chemistry retrievals require a super-solar C/O ratio that is unlikely from a planetary formation perspective. The equilibrium chemistry runs also feature a decreasing temperature structure, as opposed to the full free runs. In any case, looking at the Bayesian evidence obtained by the reference black-body fits, we conclude on the likely presence of molecular absorption in this planet.
For the transmission, which is shown in Figure Set D2, our analysis recovers a large abundance of water vapour and high-altitude opaque clouds, consistent with the previous analyses of the planet. We also find marginal evidence for FeH absorption, which is surprising as it is not detected on the day-side. This might be linked to circulation processes that prevent the molecule being seen on the day-side or it could come from an unfortunate fitting of the scattered datapoints between 1.1 µm and 1.3 µm. Constraints on the other near-optical absorbers are however relatively strong, with abundances restricted to less than 10 −5 for all of them. In terms of elemental ratios, both free and equilibrium chemistry are consistent with a large range of metallicities (from sub to super-solar). However, the C/O ratio for the free run is sub-solar, which is inconsistent with the equilibrium chemistry scenario.

D.5. Individual analysis of HAT-P-41 b
In August 2012 the HATnet survey reported the identification of three new inflated, transiting hot Jupiters orbiting bright F-type hosts. One of those systems hosts HAT-P-41 b, which is in fact part of a binary system, with a K-dwarf companion (Hartman et al. 2012;Evans et al. 2016a). The HST WFC3 transmission spectrum was studied in T18, finding the presence of water vapor and high-altitude clouds. In 2020, additional transmission data was obtained using the G280 grism (200-400 nm) onboard HST/WFC3 and analysed in conjunction with the Spitzer photometry (3.6 µm and 4.5 µm). Their best fit of the transmission spectrum, obtained using a grid of self-consistent models, confirmed the presence of clouds and water and contained evidence for molecular absorption from VO/TiO and CO 2 (Wakeford et al. 2020).
As part of our work, we obtained the HST transmission spectrum from T18 and reduced the G141 observation of the eclipse from the raw data using our Iraclis pipeline (PN: 14767, PI: Sing et al. 2016b). The Spitzer data was obtained from the G20 population study of Spitzer eclipses. We then performed our retrieval analysis of the transit and eclipse spectra with TauREx3 via the Alfnoor pipeline.
For the day-side, our full scenarios agreed on a thermal inversion in both HST and HST+Spitzer cases. The solutions display evidence for optical absorbers such as VO or FeH, but there remain large tails indicating the degeneracies between those two molecules. The full runs are, however, not statistically significant as they lead to similar Bayesian evidence to the reduced cases. In the reduced cases, the retrievals do not find evidence for any molecule and seem consistent with emission from continuous CIA and Rayleigh absorption with a decreasing thermal profile. The metallicity for the free cases is solar to sub-solar and as expected from the poor features in this spectrum, the C/O ratio of the planet remains unconstrained. The equilibrium chemistry runs are consis-tent with a decreasing to isothermal thermal structure, more similar to the free reduced case. For this planet, equilibrium chemistry provides a very good fit of the HST+Spitzer data-set, with a ln(E) comparable to the analogous free runs. Assuming equilibrium chemistry, the atmosphere is found to have a solar-like metallicity but a relatively high C/O ratio.
When analysing the terminator of the planet, we confirm the results from T18, finding a strong H 2 O feature and evidence for high-altitude clouds. In the WFC3 data, VO/TiO or CO2 are not detected, suggesting that the additional constraints found in Wakeford et al. (2020) come from the addition of the G280 grism and the Spitzer data. The metallicity of the planet is consistent with solar abundances, while the C/O ratio is difficult to constrain, ranging from sub-solar in the free run to solar with large uncertainties in the equilibrium case.

D.6. Individual analysis of HAT-P-70 b
HAT-P-70 b is a very large ultra-hot Jupiter (1.87 R J ) orbiting a A-type star that was recently detected by TESS and the HATNet program. While radial velocity measurements of the star were performed in the discovery paper Zhou et al. (2019), only a upper limit on the mass of the planet (M p < 6.78 M J ) could be obtained.
The planet was observed again in 2020 at highresolution with the HARPS-N spectrograph Bello-Arufe et al. (2021). These observations detected a number of atomic and ionic species (Ca II, Cr I, Cr II, Fe I, Fe II, H I, Mg I, Na I and V I) and constrained further the orbital configuration of the system. The planet is in a highly misaligned orbit.
In July 2021, the eclipse of HAT-P-70 b was observed with HST as part of the proposal (PN: 16307, PI: Fu et al. 2020). We analyzed the raw images from this proposal using our Iraclis pipeline. Due to the retirement of Spitzer, this planet is the only one in our sample that was not observed with this telescope. We therefore include the results from the HST only full case in the population study. Additionally, as of today, the planet's mass is still not known precisely. We tested different values, up to 6.78 M J , to verify the impact of those uncertainties on our retrievals. Due to the geometry in eclipse, we did not find major differences in our test and we therefore assume a mass of 1M J for this study.
From our retrievals, the full model recovers an inverted thermal profile with H − emission. This detection not decisive as the reduced run is not consistent with this picture and obtain a ∆ln(E) = 3.2. For the full model, the metallicity is sub-solar and the C/O ratio is consistent with solar. While the equilibrium chemistry run also presents a thermal inversion, this retrieval is consistent with a C/O ratio of 1.

D.7. Individual analysis of HD 189733 b
HD 189733 b is one of the most studied exoplanets and was the first planet to exhibit evidence for molecular absorption (Tinetti et al. 2007). The water detection for the Spitzer-IRAC data was subsequently confirmed by transit observations with HST-WFC3 (McCullough et al. 2014), while HST-NICMOS data showed evidence for methane (Swain et al. 2008b).
In eclipse ( In their analysis, they combined with the Spitzer data from Charbonneau et al. (2008), which provides photometric data up to 24 µm, a rarity in the exoplanet world. Additional eclipse observations were also obtained with HST-NICMOS (Swain et al. 2008c), Spitzer-IRS (Todorov et al. 2014Grillmair et al. 2008) and other instruments . For the Spitzer-IRS dataset, the re-analysis presented in Grillmair et al. (2008) led to differences in the recovered spectrum. Overall, the observed fluxes were consistent in terms of shape but the flux ratios appeared to be 20% lower in the later study. Sources of these discrepancies could come from the additional data considered (stellar variability), or differences in the reduction techniques (ramp models). Either way, this example highlights the difficulties of analysing combined data-sets from different sources or taken at different epochs.
As the HST transit data (PN: 12881, PI: McCullough et al. 2013) was previously analysed with Iraclis in T18, we take the data directly from this paper. For the eclipse, we perform our standard reduction of the visit and obtain the emission spectrum. We note that our recovered emission spectrum of the G141 Grism is consistent with the one from Crouzet et al. (2014). For the Spitzer data, we use the latest reduction from Charbonneau et al. (2008) consisting of the 3.6 µm, 4.5 µm, 5.8 µm, 8 µm, 16 µm and 24 µm photometric channels.
When performing the retrieval analysis of the emission spectrum, we find that the HST spectrum is almost uninformative. On the HST spectrum only, the free and equilibrium runs lead to an isothermal temperature structure and the solutions found do not lead to higher evidence than the simpler black-body fits. There is some slight indication that CH 4 might be there due to the modulation at 1.4 µm and the increased absorption around 1.6 µm, where CH 4 possesses absorbing properties. When the Spitzer data is included however, the free retrievals detect H 2 O in low abundance, around 10 −5 along with CO 2 . There might also be some hints of FeH, however, the later molecule is not required as shown by the almost similar Bayesian evidence between the reduced and full runs. The derived metallicity in the full HST+Spitzer run is slightly super-solar, which matches the findings from the equivalent equilibrium run. However, the C/O ratio is found to be 0.5 due to CO 2 being the only species detected with those observations. In the equilibrium case, the thermal profile displays a localised inversion and the C/O ratio is much more constrained to super-solar values.
Our retrieval analysis of the transmission spectrum gives similar results to the one in T18. The free retrieval indicated that the terminator is cloudy with a large water absorption feature at 1.4 µm, thus matching the findings from the eclipse retrievals that include Spitzer. The rest of the molecules are not detected but strong upper limits on the optical absorbers and H − can be extracted. The equilibrium retrieval suggests a deeper cloud cover but is consistent with a similar chemistry with solar metallicity and sub-solar to solar C/O ratio.

D.8. Individual analysis of HD 209458 b
HD 209458 b is one of the historically most-studied exoplanets, together with HD 189733 b. In 1999, it was announced as the first transiting exoplanet (Charbonneau et al. 2000;Henry et al. 2000) and shortly became the prime target of many observing campaigns due to the detection of its atmosphere via absorption of Na (Charbonneau et al. 2002). Many molecules (H 2 O, CO, CH 4 , CO 2 and NH 3 ) have now been detected (Barman 2007;Swain et al. 2008a;Snellen et al. 2010;Deming et al. 2013;Line et al. 2016;Tsiaras et al. 2018;MacDonald & Madhusudhan 2017). In particular, preliminary analyses of the day-side emission obtained with five photometric Spitzer channels indicated the possible presence of thermal inversions (Burrows et al. 2007;Knutson et al. 2008;Crossfield et al. 2012b). However, later analyses of HST G141 observations and re-analyses of the Spitzer data (Line et al. 2016;Diamond-Lowe et al. 2014) were more consistent with the presence of water vapor, seen in absorption at 1.4 µm and no thermal inversions. We re-analysed the G141 data-set (PN: 13467, PI: Bean 2013) using Iraclis and were able to recover an eclipse spectrum for all five visits. We note that Line et al. (2016) also attempted to reduce the five observed eclipses, but did not manage to recover a spectrum from the 4th visit. We note that Line et al. (2016) fitted all the visits together and did not rely on a normalization of the white light curve depths. Our pipeline does not offer this option at the moment, but future works could investigate if the observed differences are due to this point. For the Spitzer data-sets, HD 209458 b was not included as part of the population study from G20, so we used the six broadband points from 3.6 µm to 25 µm used in Diamond-Lowe et al. (2014), including the datasets reduced in Crossfield et al. (2012b) and Swain et al. (2008a).
For the terminator, the HST G141 transit observation (PN: 12181, PI: Deming 2010) was previously reduced with Iraclis in T18, so we take the transmission spectrum from this study.
We present the results of our retrievals on the eclipse spectrum of HD 209458 b in Figure Set D1 but also to illustrate our methodology in the Appendix B and Appendix C. We can immediately notice that our spectrum, as compared to the one obtained in Line et al. (2016) possesses a similar shape of the water feature between 1.2 µm and 1.3 µm, but with an apparent vertical offset. In all our free runs, except the full run with HST only, we however recover a similar solution that contains water vapour in sub-solar abundances. Due to the particularly high signal-to-noise for this spectrum, the recovered abundances are very tightly-constrained. When Spitzer is included, CH 4 and CO are also detected. For the full run with HST only, we find a very different solution that does not involve H 2 O nor CH 4 . Instead, the spectrum that was interpreted with an absorption at 1.4 µm is now interpreted with an emission feature at 1.2 µm. In this case, the temperature profile is inverted and the feature is fitted with a low amount of FeH and H − opacity. We highlight that this solution is most likely an artefact of the spectrum scatter but the difference is log(E) between the HST only runs is surprisingly high in favour of the FeH and H − solution, confirming that the shape of our WFC3 data is better fit with the inverted solution. A similarly inverted temperature solution also develops in the equilibrium chemistry retrieval, further confirming the odd shape of the water feature in our data. In all the retrievals for HD 209458 b, except the HST equilibrium chemistry run, we find sub-solar metallicity solutions that are most likely attributed to the low abundance of water we retrieve. We note that the full and equilibrium Spitzer solutions are very close, both in terms of retrieved thermal structure and elemental ratios, with a C/O of about 1.
In transit, we observe water in the atmosphere of HD 209458 b, associated with a high-altitude cloudcover. No other molecules are detected. In both the free and equilibrium runs, the atmosphere is consistent with solar metallicity. A low C/O ratio is found in the free run, while the equilibrium chemistry does not seem to provide constraints on this quantity.

D.9. Individual analysis of KELT-1 b
The first low-mass object discovered by the KELT-North survey, KELT-1 b is a 27M J , 1.12R J planet with a very short period circular orbit of 29 hours. In 2012, Siverd et al. (2012) presented spectroscopy, photometry and radial velocity data and obtained an equilibrium temperature of T eq ≈ 2400 K, due to the significant amount of stellar irradiation received by this planet. Its extreme temperature and significant inflation make KELT-1 b a valuable case-study for shortperiod atmospheric characterisation. In 2014, Beatty et al. (2014) observed KELT-1 b's secondary eclipse using Spitzer, obtaining eclipse depths which are compatible with the presence of a strong sub-stellar hotspot, suggesting poor or moderate heat redistribution for this atmosphere. Subsequently their investigations favour an atmosphere without a TiO inversion layer, where a mechanism of "day-to-night TiO cold trap" is proposed. This study was followed up with groundbased spectrophotometry in 2017, with Beatty et al.
(2017) presenting a H-band emission spectrum obtained with the LUCI1 multiobject spectrograph on the Large Binocular Telescope. Modeling of the atmospheric emission using the obtained average day-side brightness temperature of 3250 K suggested a monotonically decreasing temperature-pressure profile. The team highlighted these findings as unusual, since many other hot Jupiters of similar temperatures were believed to be in possession of either an isothermal or a thermally-inverted temperature structure. The differences were attributed to the higher surface-gravity of KELT-1 b, which could contribute to the creation of TiO cold traps.
Using our standardised methodology, we re-analysed the WFC3 eclipse observations (PN: 14664, PI: Beatty et al. 2016) of KELT-1 b. To our knowledge, this dataset has not been analysed in a publication before. Since the Spitzer observations were not included in the systematic analysis of G20, we chose to include the original data-set from Beatty et al. (2014). For the transit, since the raw data (PN: 14664, PI: Beatty et al. 2016) was already analysed with Iraclis by T18, we took the spectrum obtained in this study.
For the day-side, the spectra, temperature profiles and posterior distributions of our standardised retrievals are presented in Figure Set D1. Broadly speaking, in both HST and HST+Spitzer runs, the atmosphere presents indications of a localised thermal inversion associated with VO, FeH and H − . We note that some degeneracies exist with the posteriors presenting bi-modal behaviour and solutions that include either H − , or FeH and VO. By comparing the Bayesian evidence with the reduced runs (with no metal hydrides/oxides), our retrievals demon-strate the need for these optical absorbers, especially when Spitzer is included. In the full scenarios, the atmosphere is most consistent with a sub-solar to solar metallicity but as for many other planets, the recovered C/O ratio is not well constrained and all values are allowed. When assuming equilibrium chemistry, we note that the thermal profiles decrease with altitude, contrasting with our full free runs. With this assumption, the atmosphere prefers a solar metallicity with high C/O ratio if Spitzer is not included, or a super-solar metallicity and solar C/O ratio when Spitzer is also considered. For this planet, chemical equilibrium is not a good assumption as demonstrated by the much lower Bayesian evidence obtained by these fits.
When analysing the transmission spectrum of KELT-1 b , we do not find evidence for any molecular absorption, noting that clouds are most likely masking the presence of molecules in this atmosphere.

D.10. Individual analysis of KELT-7 b
As a planet with a high equilibrium temperature and low surface gravity (Bieryla et al. 2015), KELT-7 b is an excellent candidate for atmospheric characterisation. The transmission spectrum (PN: 14767, PI: Sing et al. 2016b) of KELT-7 b was analysed by  who found a rich transmission spectrum which is consistent with a cloud-free atmosphere and suggests the presence of H 2 O and H-. The same study also analysed the WFC3 emission spectrum (PN: 14767, PI: Sing et al. 2016b), which could be explained by a varying temperature-pressure profile, collision-induced absorption (CIA) and H-.  explored the effect of including Spitzer data from G20, finding little changes when including these observations in the emission retrievals. As the raw data was reduced using the Iraclis pipeline in , we took the data as is and started our analysis from there.
Our retrieval analysis differs from the one in  in a few aspects. First, our free retrieval includes two more molecules: CH 4 and CO 2 . Here, we also considered the appropriate Spitzer instrument response functions, which were considered flat in . The star was also modelled using an interpolated Phoenix spectrum at metallicity [Fe/H] = 0.139.
From the free run, we recover the results from , which found a thermally inverted temperature profile associated with emission of H − , in both HST only and HST+Spitzer cases. In the HST+Spitzer case, we also find hints for additional absorption from TiO. When comparing the ln(E) between reduced and full runs, we find that the addition of optical absorbers is in fact justified only for the HST+Spitzer case. As discussed previously, for an atmosphere with a high level of thermal dissociation, such as KELT-7 b, the derived metallicity is likely biased. From the equilibrium chemistry runs, we find metallicities that are solar to super-solar, but we also note that the recovered thermal structure is very different to the ones found by the freechemistry runs. In terms of C/O ratio, the equilibrium chemistry retrieval favours high values.
For the transmission retrieval, we recover the same results as . The terminator region is best characterised by H 2 O and H − absorption. Again, the derived metallicity for the free run is likely to be biased by the dissociation of the tracers used for its computation. The equilibrium chemistry retrieval on the transit spectrum does not seem to well fit the observed data and does not lead to strong constraints. D.11. Individual analysis of KELT-9 b KELT-9 b is the hottest exoplanet known so far. It orbits an A0V/B9V star (T = 10140 K) and reaches day-side temperature of about ∼5000 K, being itself hotter than many stars (Gaudi et al. 2017). Given the extreme temperatures, the majority of molecules were anticipated to be dissociated, which would leave only atomic species on the day-side and a featureless broadband spectrum (Kitzmann et al. 2018). Ground-based high resolution observations have detected a number of metals including iron, titanium and calcium (Hoeijmakers et al. 2018Yan et al. 2019;Cauley et al. 2019;Turner et al. 2020;Pino et al. 2020), which are consistent with this picture. Observations of KELT-9 b phasecurve with TESS and Spitzer have revealed an asymmetric transit (Ahlers et al. 2020) induced by the fast rotation of its host star. The rotation leads to a non-uniform structure in the star, which has a larger and brighter equator than the poles, whereas KELT-9 b orbits with a 87 • spin-orbit angle. The same studies indicate that the planet possesses a low day-night temperature contrast with indications for H 2 dissociation and recombination (Wong et al. 2019a;Mansfield et al. 2020a). As KELT-9 b is subject to an intense irradiation from the star and has a large extended hydrogen envelope reaching the Roche-Lobe limit, it is experiencing extreme atmospheric escape (Yan et al. 2019). As of today, there is no transmission data for this planet. The HST eclipse observation (PN: 15820, PI: Pino et al. 2019) was previously analysed by the same authors in  and using the same methodology, so we do not reproduce the reduction of the spectrum from this planet. We note that the physical characteristics of this system leads to a particularly precise spectrum in a single observation.
The only difference with the analysis in  is the addition of CO 2 , which is required for comparing with the rest of our population. This does not change the results and we recover a strongly-inverted temperature profile, with the presence of molecules on the day-side. This was noted in  as a rather surprising finding since the associated temperatures are high enough to dissociate them in a solar composition and equilibrium chemistry. This therefore suggested the presence of dis-equilibrium processes, non-solar chemistry or biases from other sources. When analysing the equilibrium scenarios, we note that, similarly to the reduced cases, they do not perform well and lead to much worse fits of the observed data (see the Bayesian evidence). For this planet, the estimates of the metallicity and the C/O ratio from our free chemistry are highly inaccurate due to the expected dissociation of the main molecules and our lack of constraints on the elements in atomic/ionic forms. Since the planet is likely experiencing disequilibrium chemistry, the equilibrium runs, which return a a slightly sub-solar metallicity and a C/O ratio of about 1, are also most likely biased and should be considered with caution.

D.12. Individual analysis of Kepler-13A b
The exoplanet Kepler-13A b (Shporer et al. 2011) orbits a rapidly rotating A-type star in a triple star system. The planet's day-side was studied in great detail in two previous studies (Shporer et al. 2014;Beatty et al. 2017), which indicated a very high temperature but did not find thermal inversions. Due to the presence of the two stars Kepler-13 A and Kepler-13 BC, the extraction process for this planet involves many steps that are not automated in our pipeline Iraclis. We therefore choose not to re-analyse the HST emission spectrum obtained in Beatty et al. (2017), which carefully removed the contribution from the two stars to extract the spectrum. The data was taken as part of the proposal PN:13308 led by Ming Zhao (Zhao 2013). We also added the Spitzer photometric points from Shporer et al. (2014) as-is and started our analysis from these spectra. We note that this atmosphere was not studied in transit spectroscopy with HST.
The spectra, temperature profiles and posteriors from our atmospheric retrievals can be found in Figure Set D1. As in the previous studies in the literature, we recover a decreasing temperature gradient with altitude and a well-constrained abundance of water vapour in all scenarios. This is driven by the strong and well defined 1.4 µm water feature in the HST spectrum. Given the temperature of this atmosphere, one might expect the presence of metal oxide and hydrides, but none were de-tected in our analysis. The addition of the Spitzer data does not impact the results of our retrievals, except by providing further limits on the CO 2 abundance. For the free and equilibrium runs, the metallicity is similar and consistent with solar values. The C/O ratio in the free runs is difficult to constrain due to the non detection of carbon-bearing species. When considering the equilibrium chemistry retrievals, the C/O ratio is found to be roughly solar with large uncertainties.

D.13. Individual analysis of TrES-3 b
TrES-3 b is a hot Jupiter of about 2 M jup (O'Donovan et al. 2007). At the time of its discovery in 2007, it had one of the shortest orbital periods of the known planets and was deemed a good candidate for atmospheric follow-up observations. The planet was then studied using photometric measurements, leading to some of the first exoplanet detections of emission via secondary eclipse (de Mooij & Snellen 2009;Croll et al. 2010). An eclipse spectrum was then captured by the HST-WFC3 camera and studied as part of a population analysis of five planets (Ranjan et al. 2014). According to their findings, the TrES-3 b spectrum is consistent with a simple black-body emission but the authors point out that combining these observations with the Spitzer data from Fressin et al. (2010) is inconsistent with the black-body fit. In their study, they find that an atmosphere with a solar to low metallicity best represents the combined data-set. Similar results are also found in Line et al. (2014), which is also able to rule out a high C/O ratio solution, thanks to a statistical retrieval analysis.
We reduced the eclipse observation from the raw HST images (PN: 12181, PI: Deming 2010) using our Iraclis reduction pipeline. The eclipse spectrum we obtain differs strongly from the one in Ranjan et al. (2014), here displaying an emission feature at 1.4 µm. In their study, we note, however, that the fit of their white-light curve does not present data-points covering the baseline or the ingress/egress of the transit. For our analysis, we add the Spitzer photometric observations recovered in Fressin et al. (2010) and perform our standardised retrieval runs.
Due to the higher emission obtained at 1.4 µm in our observed spectrum, the run that does not include the Spitzer data is well-fit with an emission of water. The recovered thermal profile is therefore inverted and no other molecules are necessary to explain this data-set. When the Spitzer data is added, however, the water is rejected due to the low emission in the Spitzer bandpass. The spectrum is best fit with small amounts of FeH and VO. The thermal profile in this case is more complex, with a temperature decrease with altitude deep in the atmosphere and a thermal inversion at lower pressures. The metallicity and C/O ratio derived from the free retrieval with Spitzer is not well constrained as water and carbon-bearing species are not recovered. In the equilibrium chemistry runs, the solutions also differ depending on whether Spitzer is included or not. In the HST only case, the thermal profile is inverted to best fit the 1.4 µm feature and associated with a solar metallicity and C/O ratio, most similar to the results of the corresponding free chemistry. When Spitzer is included, the profile is reversed and the spectrum is best fit with a flat spectrum in the WFC3 wavelength range, associated with a super-solar metallicity and a high C/O ratio. D.14. Individual analysis of WASP-4 b WASP-4 b is an inflated hot Jupiter (Wilson et al. 2008) with an equilibrium temperature of around 1800 K. The planet day-side was observed with HST during a single visit, taken in staring mode (PN: 12181, PI: Deming 2010). The data, while having a low signal-tonoise, favoured an atmosphere with an isothermal or a weakly-inverted thermal profile and a poor water content (Ranjan et al. 2014). Those results confirmed previous findings using Spitzer (Beerer et al. 2010). In this study, we re-analysed the HST data using our Iraclis pipeline. As the Spitzer data was not part of the G20 study, we used the data from Beerer et al. (2010). We note that an additional staring mode transit spectrum was taken as part of the same HST proposal, but an adequate fit to the data could not be obtained.
In our TauREx3 retrievals, we show that the HST spectrum does not contain a lot of information. Indeed, the pure black-body fits perform similarly to the full and reduced runs in the HST-only scenarios. When optical absorbers are considered, the retrievals indicate the possible presence of TiO, but large tails can be observed. The detection is not significant since the Bayesian evidence of the reduced runs is similar to the full cases. The thermal profile for this planet, found in our free runs, is either isothermal or locally-inverted. When running the equilibrium chemistry retrievals, it is again difficult to conclude, with the HST only case favouring a weak thermal inversion and the HST+Spitzer case favouring a decrease of temperature with altitude. In those two cases the metallicity is solar with large uncertainties, while the C/O ratio is super-solar. D.15. Individual analysis of WASP-12 b WASP-12 b has previously been studied with transmission and emission spectroscopy using WFC3 G141 in staring mode (Swain et al. 2013). The study found an emission spectrum consistent with a black-body and only marginal evidence for molecular features in the transmission data-set. Phase curve observations with the Spitzer Space Telescope at 3.6 µm and 4.5 µm indicated large emission amplitudes, showing evidence for poor day-to-night heat transport (Cowan et al. 2012). Photometric eclipse observations with Spitzer over 3.6 -8 µm (Campo et al. 2011) suggested a weak thermal inversion (Madhusudhan et al. 2011) and a super-solar C/O ratio. These studies, however, did not account for the presence of a binary companion system (Bergfors et al. 2012;Bechter et al. 2014). When re-analysed to correct for this, the previous observations of WASP-12 b were inconsistent with an isothermal atmosphere due to eclipse depth variations in the Spitzer band-passes (Crossfield et al. 2012a;, and required the presence of a carbon-rich atmosphere (HCN and C 2 H 2 ) on the day-side. On the terminator, analysis detected the possible presence of H 2 O or CH 4 /HCN and metal-bearing molecules . Later studies confirmed the possibility of super-solar C/O ratio, detecting H 2 O in the terminator (Kreidberg et al. 2015) thanks to the addition of the G102 Grism, but highlighted the impact of priors and model choices (Kreidberg et al. 2015;Oreshenko et al. 2017) in analysing this particularly difficult data-set.
STIS transit observations were also taken and, when combined with the WFC3 data, showed no evidence for TiO (Sing et al. 2013). However, as the WFC3 observations were in staring mode, the error bars on the data were large (∼120-200 ppm). The G141 transmission spectrum was subsequently obtained in scanning mode and, when combined with transit observations from the G102 grism, displayed a broad water feature (log(H 2 O) = -2.7 +1.0 −1.1 ) Kreidberg et al. (2015). The spectrum showed no signs of optical absorbers with a 3σ upper limit of log(TiO) = -3.69.
As for the other planets, we reduced the HST G141 eclipse observation (PN: 12230, PI: Swain 2010) using our Iraclis pipeline. We applied the dilution factor to correct for the un-resolved companion of this system (Bergfors et al. 2012;Bechter et al. 2014). For the Spitzer data, we obtained the corrected reductions from , which we included as-is. The transmission spectrum (PN: 13467, PI: Bean 2013) with HST was already analysed in T18, so we took the spectrum from this study.
The emission of WASP-12 b at HST wavelengths is most consistent with a black-body spectrum. This is shown by the similar Bayesian evidence obtained in our reduced and full retrievals as compared with the simpler black-body fit. When Spitzer is considered, the full retrieval recovers CO 2 and possibly FeH, with a localised thermal inversion. The FeH detection is however not robust as shown by the similar Bayesian evidence obtained in the reduced run that includes the Spitzer points. In the reduced run, CO 2 is still detected, which confirms the need for this molecule to explain the observed combined spectrum. The derived metallicity for the full run is about solar with a C/O ratio of 0.5, most likely biased by the CO 2 detection. In the equilibrium chemistry runs, the HST-only case is consistent with the results from the free chemistry retrievals. When Spitzer is added, the C/O ratio is found to be closer to 1 with a decreasing-with-altitude thermal profile. We note that the free chemistry runs do not provide a significantly better fit than the chemical equilibrium.
At the terminator, the transit of WASP-12 b is clearly indicative of high altitude clouds and the presence of water vapour, which was not detected on the day-side. Our results are similar to that of T18, and we derive an atmosphere with solar metallicity but large uncertainties on the C/O ratio. The free retrievals do not require the presence of carbon-bearing species. . These have revealed a low albedo, poor redistribution of energy to the night-side and evidence for an inverted day-side temperature-pressure profile. In particular, the analysis from Sheppard et al. (2017) considered a similar data-set to us and detected a strong thermal inversion, associated with the presence of H 2 O and CO.
The HST transmission spectrum was taken, along with two eclipses, as part of a phase curve (PN: 13467, PI: Bean 2013). As part of the same proposal, three additional eclipses were also obtained. For consistency, we re-analysed the raw data using the Iraclis pipeline and our standardised methodology. We also added the Spitzer observations from Nymeyer et al. (2011) for the four photometric channels from 3.6 µm to 8 µm. The transmission spectrum (PN: 13467, PI: Bean 2013) was included in the T18 study, which also employed the Iraclis reduction pipeline, so we conserved their spectrum.
In eclipse, we find that the the atmosphere of WASP-18 b is well fit by a localised thermal inversion, with emission from H 2 O and e − . For this planet, the results are consistent between the HST and HST+Spitzer runs. As opposed to the study from Sheppard et al. (2017), we do not find evidence for CO, differences which might be due to their model not including H − opacities or differences in our reduction techniques. We find that the atmosphere has a solar to slightly super-solar metallicity, which is also confirmed by the equilibrium chemistry retrievals. For the C/O ratio, this parameter is not constrained by our free runs due to the lack of detection of carbon-bearing species, but when assuming equilibrium chemistry, we find this atmosphere consistent with a C/O ratio of about 1.
For the transit, our retrieval analysis reveals that a very flat spectrum provides a good fit to the data. This is most likely due to the presence of high-altitude clouds, opaque at those wavelengths. As noted in T18, we are not able to extract further constraints on the possible abundances of any molecular species.

D.17. Individual analysis of WASP-19 b
WASP-19 b has been the subject of a number of investigations, from both the ground and from space. Work by Anderson et al. (2013) analysed four Spitzer eclipses, taken across 3.6-8 µm, and constructed a spectral energy distribution of the planet's day-side atmosphere. They found no stratosphere, supporting the hypothesis that hot Jupiters orbiting active stars have suppressed thermal inversions (Knutson et al. 2010). Analysis of the TESS optical phase curve showed moderately efficient day-night heat transport, with a day-side temperature of 2240 K and a day to night contrast of around 1000 K (Wong et al. 2019b). This study also utilised a host of ground-based observations by Anderson et al.  et al. (2013). However, they did not utilise the HST WFC3-G141 observations of the WASP-19 b eclipse. WASP-19 b has also been studied via transmission spectroscopy. The retrievals of the STIS-G430L, G750L, WFC-G141 and Spitzer-IRAC observations suggest the presence of water at log(H 2 O) ≈ 4 but show no evidence for optical absorbers (Sing et al. 2016a;Barstow et al. 2017;Pinhas et al. 2019). Those results do not match the ground-based transits that were acquired with the European Southern Observatory's Very Large Telescope (VLT), using the low-resolution FORS2 spectrograph which covers the entire visible-wavelength domain (0.43-1.04 µm). When analysing this data, Sedaghati et al. (2017) detected the presence of TiO to a confidence level of 7.7σ.
In this study, we reduced the WFC3 data of the WASP-19 b eclipse (PN: 13431, PI: Huitson 2013) using Iraclis and obtained the eclipse spectrum. The Spitzer data is available in G20 for the 3.6 µm and the 4.5 µm channels, but since Anderson et al. (2013) also reduced the additional 5.3 µm and 8 µm channels, we choose to use the later Spitzer data in our retrieval analysis. As the transit observations of WASP-19 b were not included in the T18 population study with Iraclis, we also reduce this data-set (PN: 12181 and 13431, PI: Deming 2010; Huitson 2013) before running our standardised retrieval analysis.
At the day-side, we find decreasing temperature profiles for all scenarios, with the data being well fit by absorption of water. In all our models, we do not find evidence for optical absorbers. This is confirmed by the ln(E), which is essentially the same in both reduced and full runs. The free chemistry runs are consistent with the equilibrium chemistry scenario and we find a solar to super-solar metallicity best fits this observed spectrum. The C/O ratio is, as expected, poorly constrained, but full runs seem to favour a low C/O ratio, which is also confirmed by the HST+Spitzer equilibrium chemistry run.
At the terminator, the transmission spectrum is well fit by water vapour. In the free run, this is the only molecule that is clearly identified. Clouds could also be present with a wide range of possible altitudes, as shown by the posterior distribution. The metallicity inferred in both free and equilibrium runs is solar. For the C/O ratio, the free run prefers a sub-solar case due to the lack of detection of carbon-bearing species, while the equilibrium run does not provide strong constraints, again demonstrating the difficulty of constraining these parameters from HST only. D.18. Individual analysis of WASP-33 b WASP-33 b is the first planet discovered to orbit a δ-Scuti variable star (Cameron et al. 2010;Herrero et al. 2011). At the time of discovery, it was the hottest known exoplanet, with temperatures above 3000 K. It now belongs to the category of the extremely-hot Jupiters, with KELT-9 b and WASP-189 b. As such its hot day-side is believed to be deprived of water vapor and contain metal oxide and hydrides. Due to the complex pulsations of the star, the analysis of WASP-33 b involves complex de-convolution of the stellar signal. Early studies of its atmosphere indicated the likely presence of a thermal inversion associated with TiO emission (Deming et al. 2012;von Essen et al. 2014;Haynes et al. 2015), which greatly contributed to the debate on the importance of TiO and VO in hot Jupiters (Fortney et al. 2008;Gandhi & Madhusudhan 2019). Later follow-ups at high resolution from the ground confirmed the detection of TiO (Nugroho et al. 2017), while observations of the termina-tor during transit were also consistent with the presence of optical absorbers, this time AlO (von Essen et al. 2019). While the picture appeared clear, recent studies and re-analysis of WASP-33 b's atmosphere, however, shed some doubts on the robustness of the TiO and VO detections (Herman et al. 2020;Serindag et al. 2021). In this work, we consider the HST-WFC3 data obtained during the eclipse of WASP-33 b (PN: 12495, PI: Deming 2011). While the spectrum was extracted in Haynes et al. (2015) with another pipeline, we do not perform an Iraclis reduction from the raw data due to the complexity of accounting for the stellar pulsations. Instead, we take the reduced spectra as-is from Haynes et al. (2015) and start our retrieval analysis there. For this planet, there are no WFC3 transit observations. The results of our retrieval analysis are shown in Figure Set D1. From comparing the Bayesian evidence of the free models, it is evidence that this reduced spectrum requires a model including optical absorbers. In the full models, the solutions are the same, independently of whether Spitzer is considered or not. The models converge to a strong thermal inversion associated with emission features of water vapour, TiO and H − . In Haynes et al. (2015), their investigations also highlighted the evidence for TiO absorption in the same data-sets, but in their models they did not fit for this molecule and choose to fix its abundance to solar values. Here we find compatible results but also highlight the likely presence of H − , an opacity that was not considered in the original study. Our solutions possess a sub-solar metallicity and unconstrained C/O ratios, but we highlight that our method to determine elemental abundances is likely inaccurate for this type of atmosphere as many of the considered tracers dissociate into atomic/ionic species, similarly to KELT-9 b. In the equilibrium chemistry run, the atmosphere is found slightly sub-solar but with a well constrained C/O ratio of about 1.

D.19. Individual analysis of WASP-43 b
WASP-43 b is one of the first (Hellier et al. 2011) and most studied hot Jupiters. With a day-side temperature of about 1800 K, this 1 R jup planet of about 2 M jup has been at the center of many studies, thanks to the large number of observations in transit, eclipse and phasecurves. In particular, it is one of the rare planets that has been observed in phase-curve with HST . This study demonstrated the use of HST in phase-curve studies for the first time and detected the presence of water vapour. A follow-up study included the Spitzer phase curve to the analysis (Stevenson et al. 2017), extending the constraints to carbon-bearing species. Those studies also highlighted an offset in the hot-spot of the planet as well as a large day-night contrast. The night-side was found to be surprisingly cool, which suggested that global opaque clouds might cover the deeper layers of this side of the planet. We note however, that the use of a more complex method to analyse the same data (Feng et al. 2016), highlighted that the previous constraints on CH 4 might be biases due to the 1D assumption in the earlier studies. Using the data from the same phase curve, water vapour was also detected in the terminator region (Kreidberg et al. 2014b) but a later, more complete study (Chubb et al. 2020a), also indicated that an optical absorber, AlO, might be present in this region.
For those reasons, WASP-43 b became the go-to planet for the theoretical work on the global circulation (Mendonça et al. 2018b;Venot et al. 2020;Carone et al. 2020) in exoplanet atmospheres and for cloud/haze modelling Helling et al. 2021). Similarly, due to the availability of consistent HST and Spitzer phase-curves, the planet is often used as benchmark for the development of retrieval techniques exploiting the 3D aspects of atmospheres (Taylor et al. 2020;Irwin et al. 2020;Changeat & Al-Refaie 2020;Feng et al. 2020). However, later contradictory re-analyses of the Spitzer phase curve demonstrated the difficulties in recovering robust estimates of atmospheric emission with Spitzer (Mendonça et al. 2018a;Morello et al. 2019;May & Stevenson 2020;Bell et al. 2020). This further demonstrates the need for cautious approaches when combining HST and Spitzer data-sets. As for the thermal profile it was believed to be decreasing with altitude, up until , which displayed model-dependent behaviour and concluded that thermal inversions could not be rejected for this planet.
In this paper, we re-analyse from scratch the eclipse observations (PN: 13467, PI: Bean 2013) of WASP-43 b. For the 3.6 µm and 4.5 µm Spitzer channels, we use the values from G20. For the transmission, the data (PN: 13467, PI: Bean 2013) was first studied in Kreidberg et al. (2014b), but it was since re-analysed in T18, which used the Iraclis pipeline. We therefore proceed from the retrieval step directly, taking the spectrum from T18.
Analysing the eclipse spectrum, we find that the HST+Spitzer retrieval contains a decreasing with altitude thermal profile. The 1.4 µm feature is, in this case, associated with the presence of water vapour, as found in previous studies. The models also require the presence of CO 2 , which absorbs in the 4.5 µm Spitzer band. Now, considering HST only, the full run, including optical absorbers, converges to a different solution.
This second solution displays a thermal inversion and low abundance of FeH, seen in emission. While this is most likely due to the data-point being scatter between 1.2 µm and 1.3 µm, this result highlights why the interpretation of eclipse spectra with HST alone is difficult. This re-enforces the discussion regarding model degeneracies made in Appendix C, which explains why a spectral feature at 1.4 µm can be well interpreted by both, absorption of H 2 O or emission of refractory species. For this planet, comparisons of the Bayesian evidence in the HST only case indicate that the simpler reduced model, with absorption of H 2 O, in fact, provides an equivalent fit to the data. Regarding the derived metallicity and C/O ratio, the retrieval on the HST+Spitzer data is consistent with a super-solar metallicity and a C/O ratio about solar. The C/O ratio value is most likely due to CO 2 being the only detected species, but the the metallicity is also confirmed by the HST+Spitzer equilibrium chemistry retrieval. The HST-only run, is best represented by an isothermal thermal profile with large uncertainties on the metallicity and C/O ratio.
When considering the transit data, we find some evidence for water vapour but the features appear relatively muted, most likely due to opaque clouds. As highlighted in previous studies (Kreidberg et al. 2014b;Stevenson et al. 2017), the temperature found at the terminator is much lower than the one from the day-side. Our findings from the free chemistry are consistent with a solar to super-solar metallicity but an unconstrained C/O ratio, which matches the results found by our equilibrium chemistry model.

D.20. Individual analysis of WASP-74 b
Orbiting a slightly evolved F9 star, WASP-74 b is a 0.95 M J moderately inflated hot Jupiter (Hellier et al. 2015). While we did not find any analysis of the eclipse of this planet, the HST transmission spectrum was first analysed in T18, uncovering a moderate water feature. Later, ground-based photometry with multiple filters was obtained and, when combined with the HST spectrum, favoured models containing TiO and VO absorption (Mancini et al. 2019). This was however argued against in a follow-up study using the same data as well as further photometric observations from the ground and Spitzer IRAC 1/2. In this later study, no evidence for these optical absorbers was found, and the data favoured a strong Rayleigh slope (Luque et al. 2020).
We take the transmission spectrum from T18 (PN: 14767, PI: Sing et al. 2016b) and reduce the single G141 observation of the eclipse (PN: 14767, PI: Sing et al. 2016b) using our Iraclis pipeline. For the Spitzer obser-vations, we use the photometric emission recovered as part of the G20 population study.
Analysing the eclipse spectrum, we do not find evidence for optical absorbers. The free retrieval including HST and Spitzer observations reveals the presence of H 2 O and CH 4 as the main absorbers with a decreasing thermal profile with altitude. It is interesting to note that the HST-only case leads to a slightly inverted thermal profile with H − opacity, however, the Bayesian evidence of this retrieval as compared to the reduced case does not change, indicating that the increase in complexity is not justified. Looking at the metallicity and C/O ratio, we find a solar metallicity and an unconstrained C/O ratio. The metallicity recovered from the equilibrium chemistry runs is also solar, and the models are overall consistent with the findings from free chemistry. When equilibrium chemistry is assumed, the thermal profile is also decreasing with altitude. When only HST is considered in the equilibrium retrievals, the C/O ratio remains unconstrained, but adding Spitzer restricts the solutions to high values.
At the terminator, we find some evidence for H 2 O and FeH, but with large tails in the posterior distributions. In both free and equilibrium models, high-altitude clouds are likely to explain the flat shape of the transmission spectrum. The metallicity is found to be around solar values with large uncertainties, while the C/O ratio is unconstrained. D.21. Individual analysis of WASP-76 b WASP-76 b was first studied by HST WFC3 in transmission. These observations were analysed as part of the T18 population study (Tsiaras et al. 2018), where their retrievals suggested a water-rich atmosphere (log(H 2 O) = −2.7 ± 1.07) with a 4.4σ detection of both TiO and VO. However, as noted in the study, the abundance of TiO retrieved is likely to be nonphysical and driven by narrow spectral coverage. Retrieval analysis of this spectrum was also performed by Fisher & Heng (2018) who extracted a water abundance which was incompatible with the previous study (log(H 2 O) = −5.3 ± 0.61). Fisher & Heng (2018) did not include TiO and VO in their analysis, but instead used a non-grey cloud model to match the opacity at shorter wavelengths. In a more recent study, Edwards et al. (2020), we accounted for a faint stellar companion to WASP-76 b (Bohn et al. 2020;Southworth et al. 2020) and showed that the original transit spectrum was contaminated. A similar study was performed by Fu et al. (2021). Once the contamination was removed, the transmission spectrum no longer showed substantial evidence for optical absorbers but retained its strong water feature. WASP-76 b has also been observed in emission with WFC3-G141 and our analysis of the corrected spectrum in the same study shows strong evidence for VO. However, the best-fit abundance of log(VO) ≈ -4 is likely unfeasible, especially as the evidence for TiO in the current data-set was marginal. In this initial study, only the WFC3 data was utilised. Here, we do not re-analyse the WFC3 data (Transit PN: 14260; Eclipse PN: 14767, PI: Sing et al. 2016b) and take it directly from Edwards et al. (2020), but we also include the Spitzer points from G20. The only differences to the retrievals performed in  are that we model the star using a Phoenix spectrum, account for the proper instrument profiles for Spitzer and include CO 2 for the consistency in our approach.
The analysis of the eclipse spectrum is consistent with a thermal inversion for all models. When optical absorbers are not included, the data is well fit with CH 4 only. When optical absorbers are included, the results are similar to Edwards et al. (2020) and we confirm that H 2 O and TiO provide the best fit to the HST data, with decisive evidence in favour of the more complex model. When HST+Spitzer are included, however, a higher ln(E) is obtained for the reduced case. This solution, while simpler, can arguably be ruled out due to the high abundance of CH 4 that would be required in this case. The free retrievals obtain solar metallicity and sub-solar C/O ratios. When running our equilibrium chemistry retrievals, we also find an inverted thermal profile and the chemistry displays a solar to super-solar metallicity as well as a high C/O ratio.
Our interpretation of the transit spectrum of WASP-76 b is the same as in . The terminator region is most consistent with the presence of water vapour and high-altitude clouds. In the free run, this leads to solar metallicity and low C/O ratio but one can see in the equilibrium chemistry run that solutions with higher metallicities and a wide range of C/O ratios are also consistent with the observed spectrum.

D.22. Individual analysis of WASP-77 A b
WASP-77 A b is an inflated hot Jupiter in a wide binary system that was first discovered in 2012 (Maxted et al. 2013). The planet orbits WASP-77 A, a G8 V star. WASP-77 B, a fainter (2 mag) K-dwarf companion, is separated by 3". As such, the spatial scans from WASP-77 A and WASP-77 B overlap somewhat on the detector, as shown in Figure D3, and cannot be separated. If not corrected for, the flux from WASP-77 B would evidently adversely affect the recovered emission spectrum for the planet. Hence, to overcome this, we utilised the specialised WFC3 simulator Wayne (Varley We started by analysing the data (PN: 16168, PI: Mansfield et al. 2020b) in the normal way, extracting the contaminated eclipse spectrum of WASP-77 A b from both visits. We then preceded to account for the contamination. In each of the two eclipses of WASP-77 A b, the relative positions of WASP-77 A and WASP-77 B were different. Therefore, we calculated the positions of WASP-77A and B for each visit from the direct image. We then extracted the high-resolution (10 angstrom bins) spectra of WASP-77A and B individually from the first non-destructive read of the out-of-eclipse observations during the second eclipse observation sequence as the star separation was better than the first. Using these spectra, we simulated high resolution spatially scanned WFC3 images for each star individually. Next, utilising the same setup used for the analysis of the real data, we extracted the stellar spectra of each star for each visit. For each bin of the eclipse spectrum, we then calculated the flux ratio between the two stars, using this to apply a correction to the eclipse spectrum of WASP-77 A b for each visit. Finally, we calculated the weighted average of the planetary emission.
As for the other planets, we included the Spitzer data. We obtained the photometric measurements in Maxted et al. (2013) and interpreted the spectrum using our standardized atmospheric retrievals. We note that Line et al. (2021) analysed observations of the same planet at high resolution and detected H 2 O and CO. In their study, they concluded from constraints on the abundance of those two molecules that the planet should have a sub-solar metallicity and a solar C/O ratio. To our knowledge the HST observations of WASP-77 A b have not been presented in previous works. We obtained the Spitzer data from the population study of Garhart et al. (2020).
Our atmospheric retrievals are consistent with the presence of water vapor and a decreasing with altitude thermal profile in both HST and HST+Spitzer runs. In the HST data, we note that additional modulations of the spectrum at 1.25 µm could be associated with the molecule TiO. The detection is significant with an ln(E) > 5 between the full and reduced models. When Spitzer is included, there are additional evidence for carbonbearing species with the absorption of CH 4 and CO 2 being detected. In the full HST+Spitzer run, we recover a solar metallicity and an unconstrained C/O ratio. Chemical equilibrium runs are consistent with a sub-solar metallicity as in Line et al. (2021) but when Spitzer data is added, the C/O ratio is about 0.8. When assuming chemical equilibrium, we note that the 1.25 µm spectral modulation that was attributed to TiO in the free runs is not well fit.

D.23. Individual analysis of WASP-79 b
WASP-79 b (Smalley et al. 2012) is a very large hot Jupiter, believed to have an evaporating atmosphere (Bourrier et al. 2015). It has been shown to have a polar orbit through the Rossiter McLaughlin effect (Addison et al. 2013).
The transmission spectrum of WASP-79 b (PN: 14767, PI: Sing et al. 2016b) has been previously analysed with Iraclis as part of Skaf et al. (2020). A water rich atmosphere was found, with a > 5σ confidence level in models with H 2 O and FeH included. The latter was included as the optical absorber of choice after initial models without it struggled with nonphysical values. Sotzen et al. (2019) also analysed the HST WFC3, combining it with observations from Magella/LDSS-3C R, TESS, and Spitzer, resulting in a transmission spectrum range of 0.6-5 µm. Whilst the spectra extracted from the LDSS-3C R data indicated clouds initially, their retrieval code ATMO found that including FeH and H-as absorbers provided better data fits. Rathcke et al. (2021) studied the most complete transmission range. This study included the STIS instrument on HST, with two transits observed through the G430L grating and another through the G750L grating, the HST WFC3 instrument, Spitzer 3.6 µm and 4.5 µm channels, and LDSS-3C R data. The study used three different re-trieval tools: NEMESIS, POSEIDON, and ATMO's retrieval tool, ARC. Similar to the results found by the other two studies, their best fitting model spectrum is characterised by an absorption feature at 1.4 µm of H 2 O, as well as the relatively smooth continuum from 0.4 to 1.3 µm indicative of H-bound free absorption.
Our analysis of the transmission data (PN: 16168, PI: Sing et al. 2016b) from HST WFC3 is shown in Figure  Set D2. Whilst the water feature is relatively well constrained in the posterior, as with the other studies mentioned, there is little evidence to support the presence of any carbon-based species. The transmission spectrum also shows detection of H − opacity, which was not included in Skaf et al. (2020). In their study, they did not include H − , which explains the differences found here. However, their Bayesian Evidence (191.2 for their best model with FeH) is very similar to ours, indicating that from an observational point of view, it is difficult to conclude on whether FeH or H − (or both) is responsible for the absorption. We note that the equilibrium retrieval does not provide a good fit of the transmission spectrum.
The emission data (PN: 14767, PI: Sing et al. 2016b) was reduced with Iraclis in Bieger et al 2021 (in prep), and is shown in Figure Set D1 along with the fits from our standardised retrievals. The emission analysis from our retrievals supports what is found in the transmission-again, the most unconstrained parameters are carbon-based species. Inspection of the full retrieval indicates the detection of VO and FeH, associated with an increasing with altitude thermal profile. In the full model, the spectral modulation in the HST spectrum is fit with an emission feature from optical absorbers. On the opposite, if optical absorbers are not included, or if additional constraints from chemical equilibrium are considered, the HST spectrum is better fit with a decreasing temperature profile and the absorption of water vapour. This behaviour is independent from whether Spitzer is considered or not and is due to the degeneracies induced by the narrow spectral region of HST, as discussed in Appendix C. The full retrievals, however, provide a strong evidence in favor of optical absorbers as the difference in Bayesian Evidence is larger than 3 (∆ln(E)=4.5).
These results will require further retrievals and analysis in order to understand the role of TiO, VO, FeH or H-in the atmosphere of WASP-79 b. The JWST Early Release Science (ERS) program includes WASP-79b as a primary target; JWST will be observing it for 42 hours over four different modes, which will present more opportunity to study this planet in depth with more precise data. More detailed analysis of WASP-79 b will be presented in Bieger et al 2021 (in prep).

D.24. Individual analysis of WASP-103 b
WASP-103 b is an ultra-short period planet (P = 22.2 hr) whose orbital distance is less than 20% larger than its Roche radius, resulting in the possibility of tidal distortions and mass-loss via Roche-lobe overflow (Gillon et al. 2014b).
WASP-103 b's HST WFC3 emission spectrum was found to be featureless down to a sensitivity of 175 ppm, showing a shallow slope toward the red (Cartier et al. 2017). Work by Manjavacas et al. (2019), which performed a reanalysis of the same data-set, found that the emission spectrum of WASP-103 b was comparable to that of an M-3 dwarf. Delrez et al. (2018) obtained several ground-based high-precision photometric eclipse observations which, when added to the HST data, could be fit with an isothermal black-body or with a low water abundance atmosphere with a thermal inversion. However, their Ks band observation showed an excess of emission compared to both these models. Recently, a phase-curve analysis of the planet was taken and reported in Kreidberg et al. (2018). The study also utilized the previous HST emission spectra and confirmed a seemingly featureless day-side. The phase-resolved spectra of WASP-103 b were compared in this work to those of brown dwarfs and directly-imaged companions of similar temperature which show evidence for water absorption at 1.4 µm, whereas WASP-103 b showed no such feature. The result could be due to WASP-103 b's irradiation environment and its low surface gravity. A later study on the same data (Changeat 2022) and employing a unified phase-curve retrieval technique obtained a more complex picture of the planet. It confirmed the presence of thermal inversion and dissociation processes on the day-side of the planet and found signature of FeH emission. The study also constrained water vapor across the entire atmosphere. Ground-based transmission observations found strong evidence for sodium and potassium (Lendl et al. 2017).
Two HST G141 phase curves of WASP-103 b were obtained (PN: 14050, PI: Kreidberg et al. 2014a) which each contained a single transit and eclipse. Additionally, two further eclipses were taken with the same instrument (PN: 13660, PI: Zhao et al. 2016).
Our eclipse spectrum of WASP-103 b, is not consistent with a pure black-body fit, thus contrasting with the previous studies. In our re-analysis we find that our models prefer a thermal inversion with the presence of optical absorber such as VO or FeH. When Spitzer is included, the model also includes H − , which was predicted to be an important opacity source for this atmosphere (Kreidberg et al. 2018). We note, however, that compar-ing the Bayesian evidence with the reduced models does not provide decisive evidence in favour of the full models, meaning that an atmosphere without these active molecules could also well fit this data-set. The metallicity associated with our full runs is centred around solar values but with large uncertainties. For the C/O ratio, since CO 2 is detected in the HST+Spitzer case, the derived value converges to 0.5, but this is most likely due to the lack of detection of other species. When considering the equilibrium chemistry case, a thermal inversion is also inferred from the data. This is associated with a solar metallicity and a C/O ratio of about 1.
Analysing the transit spectrum (PN: 14050, PI: Kreidberg et al. 2014a), we observed that the downward slope is well fit by VO. Other optical absorbers might be present (TiO, H − ) but the data does not allow verification of this. The solution found possesses a wide range of metallicities; sub-solar in nature. When applying equilibrium chemistry to this data-set, the thermal profile converges to unrealistic values for this atmosphere.

D.25. Individual analysis of WASP-121 b
Significant observational time has been spent on WASP-121 b (Delrez et al. 2016). In transmission, the analyses of ground-based observations, Hubble-STIS and Hubble-WFC3 have shown the presence of H 2 O and optical absorption attributed to VO and/or FeH (Evans et al. 2016b. The authors of these studies note that chemical equilibrium models with solar abundances cannot reproduce the spectrum seen, while free chemical retrievals can only do so by converging to high abundances of VO and FeH. WASP-121 b has also been observed in eclipse and the presence of a thermal inversion provides the best fit to the data (Bourrier et al. 2019;Daylan et al. 2019;Parmentier et al. 2018;Mikal-Evans et al. 2020). Bourrier et al. (2019) attributed this to VO with a best-fit abundance of log(VO) = -6.03 while Mikal-Evans et al. (2020) performed chemical equilibrium retrievals, finding evidence for a muted water feature due to dissociation and H − opacity.
In parallel, high resolution ground-based observations of the transit have put upper limits on the abundances of TiO and VO at the terminator with log(VMR) < -7.3 and 7.9 respectively, suggesting these cannot be causing the inversion seen (Merritt et al. 2020). However, the study highlighted that these limits are largely degenerate with other atmospheric properties such as the scattering properties or the altitude of clouds on WASP-121 b. Another study found a host of atomic metals, including V, which are predicted to exist if a planet is in equilibrium and has a significant quantity of VO (Hoeijmakers et al. 2020). They too noted the absence of TiO