Evidence for a Dayside Thermal Inversion and High Metallicity for the Hot Jupiter WASP-18b

Hot Jupiters have been vital in revealing the structural and atmospheric diversity of gas-rich planets (see recent reviews by Crossfield 2015; Madhusudhan et al. 2016; Deming & Seager 2017). Since they are exposed to extreme conditions and relatively easy to observe through transit and eclipse spectroscopy, hot Jupiters provide a window into a unique part of parameter space, allowing us to better understand both atmospheric physics and planetary structure.

An outstanding question that has emerged for highly irradiated planets is the presence and origin of stratospheric thermal inversions, which have been detected in several extremely irradiated hot Jupiters (Haynes et al. 2015; Evans et al. 2017). Hubeny et al. (2003) predicted that thermal inversions in highly irradiated atmospheres would be caused by the presence of optical absorbers (e.g., TiO and VO) high in the atmosphere, but there may be other causes such as insufficient cooling (Mollière et al. 2015) or sulfur-based aerosols (Zahnle et al. 2009).

Constraints on the structure and composition of exoplanetary atmospheres allow us to test, refine, and generalize planetary formation models. Volatile ices are expected to play an important role in planet formation; thus a constraint on the composition of a hot planet's atmosphere gives us insight on how and where it was formed (Öberg et al. 2011; Madhusudhan et al. 2014). In our solar system there is an inverse mass versus atmospheric metallicity relationship, and whether or not it extends to exoplanets is informative to planetary formation and migration models. There is some evidence that the trend holds (Kreidberg et al. 2014), however that parameter space is not yet sufficiently populated to enable firm conclusions.

In this paper we use Hubble Space Telescope (HST) spectroscopy and Spitzer/IRAC photometry of secondary eclipses to explore the thermal structure and composition of the dayside atmosphere of WASP-18b, an extremely hot (${T}_{\mathrm{eq}}=2411$ K) and massive ($M=10.3\,{M}_{J}$) hot Jupiter orbiting an F-type star with an orbital period of less than one day.

We used Wide Field Camera-3 (WFC3) observations of five secondary eclipses of WASP-18b from the HST Treasury program GO-13467 (PI: J. Bean). WFC3 obtains low resolution slitless spectroscopy from 1.1 to 1.7 μm using the G141 grism (R = 130), as well as an image for wavelength calibration using the F140W filter. Grism observations were taken in the spatial scan mode (Deming et al. 2013) with a forward-reverse cadence (Kreidberg et al. 2014). The first three visits, taken between 2014 April–June, are single eclipse events. Visit 4, taken in 2014 August, contains two eclipses in an orbital phase curve, and we extract those eclipses and analyze them separately.

We also re-analyze two eclipse observations of WASP-18b taken in the 3.6 μm and 4.5 μm channels of the Spitzer Space Telescope's IRAC instrument (Program ID 60185). The 3.6 μm observation was performed on 2010 January 23, while the 4.5 μm observation was taken 2010 August 23. Both observations were taken using an exposure time of 0.36 s in the subarray mode, and were first analyzed in Maxted et al. (2013).

Our grism spectroscopy analysis utilized HST "ima" data files. We separated the data by scan direction, removed background flux, and corrected for cosmic rays and bad pixels. We removed background flux via the "difference frames" method outlined in the appendix of Deming et al. (2013), and set the aperture to maximize the amount of source photons in our analysis. The end result is two reduced light curves—one forward scan and one reverse scan—for each eclipse, which we analyze separately.

The F140W photometric image determines the location of the zero-point, which we used to assign a wavelength to each column. We confirmed the wavelengths by fitting a T = 6400 K, log g = 4.3, [Fe/H] = 0.1 ATLAS stellar spectrum (Castelli & Kurucz 2004), multiplied by the grism sensitivity curve, to an observed in-eclipse spectrum.

3.1. Light Curve Analysis

Empirical methods are necessary to correct for non-astrophysical systematic effects in WFC3 spectroscopy (Berta et al. 2012; Haynes et al. 2015). Correction methodology is especially important in emission spectroscopy, where the magnitude of systematic effects can be greater than the eclipse depth (Kreidberg et al. 2014). We thereby combined two strategies: initial removal of systematic trends using parametric marginalization (Gibson 2014; Wakeford et al. 2016), and further detrending by subtraction of scaled band-integrated residuals from wavelength bins (Mandell et al. 2013; Haynes et al. 2015). Our method accounts for uncertainty in instrument model selection, and residuals from the band-integrated analysis allow us to utilize the normally excluded first orbit of each HST data set in our spectroscopic analysis.

Fitting a band-integrated light curve provides residuals that we use to remove unidentified systematics from the spectrally resolved light curves. We calculate the HST phase (parameter for ramp and HST breathing), planetary phase (parameter for visit-long slope), and a wavelength shift derived by cross correlating each spectrum with the last spectrum for the visit (parameter for jitter) for each exposure in a time series. The grid of systematic models comprises a combination of a linear planetary phase correction and up to four powers of HST phase and wavelength shift. These models are then multiplied by a Mandel & Agol (2002) eclipse model. We simultaneously fit for the eclipse depth, all systematic coefficients, and—for two light curves with ingress and egress points—the center of the eclipse time. All of the other system parameters are fixed to literature values.

We use a Levenberg–Markwardt (L–M) least squares minimization algorithm (Markwardt 2009) to determine the parameter values, since Wakeford et al. (2016) found it to agree with the more computationally expensive Monte Carlo Markov Chain (MCMC) analysis to within 10%. An example band-integrated light curve with systematic effects removed using the best-fitting model is shown in the leftmost panel of Figure 1. The scatter (rms) of the residuals of the band-integrated curves ranges from 1.3–5.5 times the photon noise, indicating that there is excess noise beyond the photon limit present. Excess noise in the band-integrated curves is also shown by comparisons of the cumulative distributions of residuals with those of a photon-limited Gaussian (see the bottom left panel of Figure 1). However, the structure of this excess noise does not change with wavelength, allowing for its removal from the corresponding spectral light curves.

Zoom In Zoom Out Reset image size

To derive the emission spectrum, we bin the exposures in wavelength between the steep edges of the grism response and fit these spectrally resolved light curves. We remove wavelength-dependent systematics by fitting each spectral bin separately in a process that mimics the band-integrated process, with three exceptions. First, the eclipse mid-time is now fixed to the value determined by the band-integrated analysis. Second, it is possible that shifts on the detector are wavelength-dependent, so the jitter parameter is recalculated for each wavelength bin using only that portion of the spectrum in the cross-correlation procedure. Third, each systematic model now incorporates the residuals from the band-integrated fit of the same model as a decorrelation variable. The amplitude of the residuals is a free parameter, although the shape is assumed to be constant in wavelength. This removes any remaining wavelength-independent trends in the data. An example result of a reduced spectral bin light curve is shown in the central panel of Figure 1.

Finally, eclipse depths from the multiple visits are combined via an inverse-variance-weighted mean, giving the emission spectrum for WASP-18b. The spectra for all visits are shown in Figure 2.

Zoom In Zoom Out Reset image size

The average rms of the systematic-reduced spectroscopic light curves is 1.04 times the photon noise and the median rms is 0.97 times the photon noise, indicating that shot noise is typically the dominant error source. The close agreement between the cumulative distributions of residuals and those of a Gaussian with a width determined by the photon noise provides further evidence that the analysis achieved photon-limited results for the vast majority of spectral curves (see the bottom center panel of Figure 1). The remaining spectral curves have residuals with an rms greater than 1.5 times the photon limit, indicating that excess noise is present. These only constitute 6% of all spectral bins, and every one is from the single eclipse observation taken in May. We explored removing the May data set due to this increased noise, but the exclusion of these data did not affect the variance-weighted spectrum, and we chose to include this visit in subsequent analyses. Figure 2 contains the emission spectra from every visit, demonstrating the consistency of the structure of the spectrum. Our analysis routine finds that the outlier depths from the May visit have very high errors due to the presence of correlated noise, and so they contribute very little to the weighted spectrum.

To further check our methodology, we reanalyzed the published emission spectra for WASP-43b, WASP-103b, and WASP-121b. We find an agreement to the published spectra, with a mean point-by-point variation (difference/uncertainty) of 89%, 23%, and 50% for the three data sets, respectively, demonstrating the consistency of our analysis pipeline with those published by other authors.

Spitzer secondary eclipse measurements of WASP-18b were reported by Maxted et al. (2013), and we have reanalyzed key portions of those data. We confine our reanalysis to the 3.6 and 4.5 μm bands, because the instrumental systematic errors are greatest in those bands, and there are new methods to correct those systematics.

We use an updated version (Tamburo et al. 2017) of the Pixel-level Decorrelation framework (Deming et al. 2015). Our photometry uses 11 different circular aperture sizes (with radii ranging from 1.6 to 3.5 pixels). We decorrelate the instrumental systematics while simultaneously fitting for the eclipse depth, using binned data, as advocated by Deming et al. (2015) and Kammer et al. (2015). The fitting code selects the optimal aperture and bin size, and obtains an initial estimate of the eclipse depth and the pixel basis vector coefficients using linear regression. We then implement an MCMC procedure (Ford 2005) to explore parameter space, refine the best-fit values, and determine the errors. At each step, we allow the central phase, orbital inclination, and eclipse depth to vary, but lock all of the other orbital parameters to the values used in the WFC3 analysis. We also vary the multiplicative coefficients of our basis pixels (see Deming et al. 2015) and visit-long quadratic temporal baseline coefficients at every step. Our best fits use aperture radii of 2.0 and 2.5 pixels, and bin sizes of 76 and 116 points at 3.6 and 4.5 μm, respectively. The scatter in the binned data, after the removal of the best-fit eclipse, is 1.01 and 0.95 times the photon noise at 3.6 and 4.5 μm, respectively, those ratios being statistically indistinguishable from unity.

We ran three chains of 500,000 steps for both bands, confirming their convergence through the Gelman-Rubin statistic (Gelman & Rubin 1992). We combine all of the chains of eclipse depth into a unified posterior distribution for each band, and fit a Gaussian to this distribution to determine the error on eclipse depth. Our results are included in Table 1, and exhibit excellent agreement with Maxted et al. (2013), but with smaller errors.

We use our WFC3 spectrum along with the Spitzer and ground-based Ks-band photometry to constrain the composition and temperature structure of the dayside atmosphere of WASP-18b. We use the HyDRA retrieval code (Gandhi & Madhusudhan 2017b), which comprises a thermal emission model of an atmosphere coupled with a nested sampling algorithm for Bayesian inference and parameter estimation. The forward model, based on standard prescriptions for retrieval (Madhusudhan & Seager 2009; Madhusudhan et al. 2011), computes line-by-line radiative transfer in a plane parallel atmosphere under the assumptions of hydrostatic equilibrium and local thermodynamic equilibrium. The pressure–temperature (P–T) profile and chemical compositions are free parameters in the model.

The model includes 14 free parameters. For the P–T profile, we use the parametrization of Madhusudhan & Seager (2009), which involves six free parameters. The atmosphere comprises 100 layers equally spaced in log pressure between 10⁻⁶ bar and 10² bar. For the atmospheric composition we consider several species expected to be prevalent in very hot Jupiter atmospheres and with significant opacity in the observed spectral range (Madhusudhan 2012; Moses et al. 2013; Venot & Agúndez 2015). This includes H₂O, CO, CH₄, CO₂, HCN, C₂H₂, TiO, and VO. The uniform mixing ratio of each species are free parameters in the model. We assume an H₂/He-rich atmosphere with a solar He/H₂ ratio of 0.17. We consider line absorption from each of these species and collision-induced opacity from H₂–H₂ and H₂–He. The sources of opacity data are described in Gandhi & Madhusudhan (2017a); the molecular line lists are primarily from EXOMOL (Tennyson et al. 2016) and HITEMP (Rothman et al. 2010), and the CIA opacities are from Richard et al. (2012). The retrieval explores model parameter space with Bayesian nested sampling using the MultiNest code via the Python wrapper, PyMultiNest (Skilling 2004; Feroz et al. 2013; Buchner et al. 2014). We sample the multi-dimensional parameter space using 4000 live points for a total of more than one million model evaluations.

Our best-fit retrieval requires a strong thermal inversion in the dayside atmosphere. The bottom inset of Figure 3 shows the retrieved P–T profile with confidence contours, indicating an upper atmospheric temperature increase. The requirement of a thermal inversion is guided by the strong emission inferred in the 4.5 μm Spitzer IRAC band, with a brightness temperature of 3100 ± 50 K, which is significantly higher than the rest of the data. This can be explained by the presence of a thermal inversion in the atmosphere along with the presence of either CO or CO₂, which both exhibit pronounced spectral features in the 4.5 μm band (Burrows et al. 2007; Fortney et al. 2008; Madhusudhan & Seager 2010). We break this degeneracy by requiring that CO₂ be less than H₂O as expected for hot Jupiter atmospheres (Madhusudhan 2012; Heng & Lyons 2016). Another subtlety is the apparent minor trough near ∼1.6 μm, which we attribute to CO absorption below the inversion layer (∼1–10 bar), where temperature decreases outward. Emission in the 4.5 μm band is due to CO in the 0.001–0.1 bar range, which contains the thermal inversion. As part of the nested sampling analysis, we compute the Bayesian evidence value for our retrieved spectrum. By comparing this value with that obtained for a model without a thermal inversion, we conclude that a thermal inversion is favored at the 6.3σ significance level. Similarly, comparison to a model lacking CO implies that the presence of CO is favored at the 6.1σ level. Interestingly, the transition point of the inversion occurs at 0.1 bar, which is characteristic of all of the planets in the solar system with inversions as well as models of hot Jupiters (Madhusudhan & Seager 2009; Robinson & Catling 2014).

Zoom In Zoom Out Reset image size

Figure 4 shows the posterior probability distributions of all of the model parameters. The data require a CO volume mixing ratio of ${19}_{-8}^{+18}$% in the atmosphere, which is ${380}_{-160}^{+360}$ times the amount expected for a solar-abundance atmosphere at this temperature in thermochemical equilibrium. The high CO abundance is primarily constrained by the emission required to explain the 4.5 μm IRAC point as well as the absorption trough in the WFC3 band at 1.6–1.7 μm. We detect no other chemical species (see Figure 4). In particular, the non-detection of H₂O at both 1.4 and 6 μm provides a robust 3σ upper-limit of 10⁻⁶ on the volume mixing ratio. The sum-total of constraints on the chemical species lead to a super-solar metallicity in the planet (${\rm{C}}/{\rm{H}}={\rm{O}}/{\rm{H}}={283}_{-138}^{+395}\times$ solar O/H) and a C/O ratio of ∼1.

Zoom In Zoom Out Reset image size

We also conducted free-chemistry retrievals with no priors on the CO₂ abundance and find the same key results. For both cases, the data require a strong thermal inversion, a C/O ratio of ∼1, and a super-solar metallicity.

The constraints on the chemical abundances are consistent with expectations for a high C/O ratio atmosphere in the high temperature regime of WASP-18b (Madhusudhan 2012; Moses et al. 2013) where chemical equilibrium is expected to be satisfied. At high temperatures, H₂O is expected to be the most dominant oxygen-bearing molecule for a solar-abundance elemental composition (e.g., with a C/O = 0.5; Madhusudhan 2012; Moses et al. 2013). In contrast, the low-abundance of H₂O observed is possible only if the overall metallicity and O abundance were low, or if the C/O ratio were high. Given the high abundance of CO we retrieve, the only plausible solution is both a high oxygen abundance and a high C/O ratio. The constraints on all of the other species are also consistent with this scenario. While we cannot rule out a contribution from CO₂ emission in the 4.5 μm Spitzer band, the high abundance of CO₂ needed would be chemically inconsistent with the non-detection of H₂O, and we therefore believe this scenario to be unlikely.

Our inferences for this planet indicate an unusual atmosphere in several respects, calling for comment on the reliability of our conclusions. While the inference of a temperature inversion per se is no longer surprising for strongly irradiated planets (Haynes et al. 2015; Evans et al. 2017), both the very high metallicity and C/O ∼ 1 have less precedent. Those aspects are forced upon us by the lack of observed water in the WFC3 and Spitzer bandpasses, by the slight decrease at the long end of the WFC3 band, and by the Spitzer photometry point at 4.5 μm. The non-detection of WFC3 water is certainly robust—several independent eclipses show no sign of the band head that should occur at 1.35 μm (Figure 2). We reiterate that the inference of a thermal inversion hinges critically on the single Spitzer photometric point at 4.5 μm. Previously, Nymeyer et al. (2011) postulated a temperature inversion for exactly that reason. Since our eclipse depth agrees with those from previous analyses (Nymeyer et al. 2011; Maxted et al. 2013), we consider this measurement robust with regard to the analysis technique. Nevertheless, the photometry does not reveal the resolved band structure of the 4.5 μm CO band in emission that would lead to an unequivocal detection of molecular emission. However, given the data we have and the successful checks on our data analysis procedures, the unusual atmosphere of WASP-18b is a compelling conclusion. Our observations also reveal the first instance where both absorption and emission features are seen in the spectrum of an exoplanet, both due to CO. The absorption at ∼1.6 μm is caused by a weaker CO band compared to the emission in a stronger CO band in the 4.5 μm region. As shown by the contribution functions in Figure 3, the 1.6 μm region in the spectrum probes the lower atmosphere due to the lower opacity compared to the 4.5 μm band, which probes the upper atmosphere due to a higher opacity in that spectral region. Note that simultaneous absorption and emission in the same molecule is observed in the Earth's infrared spectrum, specifically in the 15 μm band of CO₂, due to the temperature structure at the tropopause and stratosphere (Hanel et al. 1972).

If confirmed, the atmospheric properties of WASP-18b open a new regime in the phase space of hot Jupiters. Classically, thermal inversions in hot Jupiters were suggested to be caused by TiO and VO in very high temperature atmospheres (Hubeny et al. 2003; Fortney et al. 2008). All studies so far have focused on the plausibility of TiO/VO as a function of various parameters and processes such as settling and cold traps (Spiegel et al. 2009), stellar chromospheric emission (Knutson et al. 2010), C/O ratio (Madhusudhan et al. 2011), dynamics (Menou 2012; Parmentier et al. 2013), etc. For TiO/VO to be abundant enough to cause thermal inversions, the C/O balance must be approximately 0.5 or lower (Madhusudhan et al. 2011). Planets with high C/O ratios were not predicted to host thermal inversions since their TiO/VO abundances would be severely depleted (Madhusudhan 2012); however, recent work suggests other processes, such as inefficient atmospheric cooling, could lead to an inverson (Mollière et al. 2015). Alternatively, oxygen-poor absorbers may play a similar role to TiO and VO (Zahnle et al. 2009). The two hot Jupiters for which thermal inversions have been detected have both showed signatures of TiO/VO in their atmospheres: WASP-33b (Haynes et al. 2015) and WASP-121 (Evans et al. 2017). WASP-18b is the first system that shows a thermal inversion along with a high C/O ratio of ∼1 with no evidence for TiO/VO, and hence provides a new test case for theories of thermal inversions in hot Jupiters.

WASP-18b's unique atmospheric composition implies an interesting constraint for planetary formation theories. Its metal enrichment is a factor of 1000 more than that predicted by the inverse mass–metallicity relationship for a 10M_J planet (Kreidberg et al. 2014). High metallicity and a C/O ratio of 1 are plausibly explained by formation from extremely CO-rich gas beyond the water condensation line (Madhusudhan et al. 2014) or upper atmospheric enrichment in carbon and oxygen due to ablation of icy planetesimals during late-stage accretion (Pinhas et al. 2016). Future eclipse observations with the James Webb Space Telescope and improved modeling of giant planet accretion processes will help clarify the details of WASP-18b's formation history.

These observations are associated with HST program GO-13467 (PI: J. Bean) and Spitzer program GO-60185 (PI: P. Maxted). The authors thank Hannah Wakeford for assistance on data reduction algorithms and Eric Lopez for helpful discussions on interpretation. K.S., A.M., P.T., and D.D. acknowledge funding from HST Program GO-14260 and a grant from the NASA Astrophysics Data Analysis Program. S.G. acknowledges financial support from the Science and Technology Facilities Council (STFC), UK. A.P. is grateful for research support from the Gates Cambridge Trust.