Evidence of a Clear Atmosphere for WASP-62b: The Only Known Transiting Gas Giant in the JWST Continuous Viewing Zone

The Hubble/STIS transit observations consist of low-resolution (∼500) time series spectra collected on UT 2018 January 26 (visit 57) and UT 2018 May 12 (visit 58) with the G430L (2892–5700 Å) grism, and on UT 2017 November 21 (visit 59) with the G750L (5240-10270 Å) grism. Each visit consisted of five consecutive 96 minutes orbits, during which 48 stellar spectra were obtained over exposure times of 253 s. To decrease the readout times between exposures, we used a 128 pixel wide subarray. The data were taken with the 52 ×2 arcsec² slit to minimize slit light losses. We note that both of the G430L visits suffered from guide star issues at the beginning of the observations, resulting in the loss of data collection during the first orbit (visit 57) or part of the first orbit (visit 58). From the subsequent exposures taken during these visits, however, we were able to extract good quality light curves with a median precision of 311 ppm.

We reduced the STIS spectra using the methods described in Alam et al. (2018, 2020). Briefly, we bias-, dark-, and flat-field corrected the raw 2D data frames using the CALSTIS pipeline (V 3.4). We corrected for cosmic-ray events using median-combined difference images to flag and interpolate over bad pixels. We performed a 1D spectral extraction from the calibrated .flt files and extracted light curves using an aperture width of 13 pixels. From the x1d files, we obtained a wavelength solution by resampling all of the extracted spectra and cross-correlating them to a common rest frame. The cross-correlation measures the shift of each stellar spectrum with respect to the first spectrum of the time series, so we resampled the spectra to align them and remove subpixel drifts associated with the different locations of the spacecraft on its orbit (Huitson et al. 2013).

We observed two transits with Spitzer/IRAC on UT 2016 November 24 and UT 2016 December 7 in the 3.6 μm and 4.5 μm channels, respectively (Fazio et al. 2004; Werner et al. 2004). Each IRAC exposure was taken over integration times of 2 s, resulting in 20,159 images. We reduced the Spitzer photometry using the techniques of Nikolov et al. (2015), Sing et al. (2015, 2016), and Alam et al. (2018). We filtered for outliers in the data and subtracted the sky background using an iterative 3σ clipping procedure, as outlined in Knutson et al. (2012) and Todorov et al. (2013).

We extracted photometric points using two approaches: fixed and time-variable aperture photometry. In the fixed approach, we used circular apertures with radii ranging from 4 to 8 pixels in increments of 0.5. For the time-variable approach, we scaled the size of the aperture by the noise pixel parameter (the normalized effective background area of the IRAC point response function), which depends on the FWHM of the stellar PSF squared (Mighell 2005; Knutson et al. 2012; Lewis et al. 2013; Nikolov et al. 2015). We compared the light-curve residuals and the white and red noise components measured with the Carter & Winn (2009) wavelet technique to identify the best results from both methods.

We also observed one transit of WASP-62b with Hubble/WFC3 on UT 2017 April 14 (GO: 14767). This visit, however, suffered from severe guide star issues that impacted the data quality of the observations, rendering it unsuitable for reliably extracting information on the planet's atmosphere. Guide star issues during the planet's transit impacted the telescope guiding, resulting in large drifts on the order of several pixels that affect the photometric quality of the spectroscopic channels and contribute to systematic trends that drift in wavelength during the observation (see Figure 11 of Sing 2018). Further, the G141 stellar spectra exhibit large offsets in the stellar spectral structure, even though the edges of the spectrum are sharp, likely due to warping of the scan as the scan rate changes across the detector as well as positional shifts in the scan.¹⁴ Although a spectrum extracted from this observation was presented in Skaf et al. (2020), we do not use this G141 data set in our analysis for the reasons outlined above.

We were awarded XMM-Newton time (program ID 80479, P.I. J. Sanz-Forcada) to observe WASP-62 on UT 2017 April 14 for ∼8 ks. We fit European Photon Imaging Camera (EPIC) spectra using a one-temperature coronal model to calculate an X-ray luminosity (0.12−2.48 keV) of L_X =6.0 × 10²⁸ erg s⁻¹ (S/N = 5.8), for a Gaia DR2 distance of 176.53 ± 0.41 pc. The resulting $\mathrm{log}{L}_{{\rm{X}}}/{L}_{\mathrm{bol}}=-5.1$ measurement indicates that the star has a moderate activity level, similar to the young solar analog ι Hor (G0V) in which X-ray variability is present but flares are not frequently observed (Sanz-Forcada et al. 2019).

Furthermore, the X-ray and UV light curves taken with the XMM-Newton EPIC (∼1−124 Å) and the Optical Monitor (OM)/UVM2 filter (2070−2550 Å) show some level of variability.¹⁵ The UV light curve shows some dips (up to ∼10% absorption) around the beginning of the primary transit that could indicate the presence of occulted chromospheric plages as the planet transits. (Further details will be included in J. Sanz-Forcada et al. 2021, in preparation). Considering the activity level of the host star, we therefore account for starspots and faculae in our retrieval analysis (Section 4).

We fit the light curves following the procedure detailed in Kirk et al. (2017, 2018, 2019), which we briefly describe here. We modeled the analytic transit light curves of Mandel & Agol (2002) using the Batman package (Kreidberg 2015), combined with a Gaussian process (GP) implemented with the george (Ambikasaran et al. 2015) code to model noise in the data. GPs are increasingly used in transmission spectroscopy (e.g., Gibson et al. 2012a, 2012b; Evans et al. 2015, 2017; Kirk et al. 2017, 2018, 2019; Louden et al. 2017; Parviainen et al. 2018).

For both the white light and spectroscopic light-curve fits, we used the four parameter nonlinear limb darkening law (Claret 2000). We derived the limb darkening coefficients using 3D stellar models (Magic et al. 2015), and fixed them to the theoretical values in our fits. For the white light-curve fits, we fixed the inclination i to 883, the scaled semimajor axis a/R_⋆ to 9.53, and the period P to 4.4 days (based on the literature values from Hellier et al. 2012), and fit for the time of midtransit T₀, planet-to-star radius ratio R_p/R_⋆, and the GP hyperparameters. Our GP was defined by the combination of 11 squared exponential kernels, each operating on a jitter detrending vector.¹⁶ These 11 kernels each had their own length scale but shared a common amplitude. We additionally included a white noise kernel, defined by the variance (σ²), to account for white noise unaccounted by the photometric error bars. The GP therefore added an additional 13 free parameters for each light curve.

Following similar studies (e.g., Evans et al. 2017, 2018), we standardized each GP input variable (jitter detrending variable) by subtracting the mean and dividing by the standard deviation. This method gives each standardized variable a mean of zero and a standard deviation of unity, which helps the GP to determine the inputs of importance for describing the noise characteristics. We fit for the natural log of the inverse length scale (e.g., Evans et al. 2017, 2018; Gibson et al. 2017). We placed wide, truncated uniform priors in log space on each of the hyperparameters. This choice assures uninformative priors since a uniform prior in log-space is akin to fitting for a 1/x prior in nonlogarithmic space. Since we are fitting for the natural log of the inverse length scale, this choice encourages the length scale to longer length scales so as to ensure that the GP does not over-fit the data (e.g., Evans et al. 2018; Parviainen et al. 2018; Gibson et al. 2019). The GP amplitude was bounded between 0.01 and 100× the variance of the out-of-transit data, and the length scales were bounded by the typical spacing between data points and 5× the maximum length scale. The white noise variance was bounded between 10⁻¹⁰ and 25 ppm.

We ran a Markov Chain Monte Carlo (MCMC) using the emcee Python package (Foreman-Mackey et al. 2013) to perform the fitting. For all of the light-curve fits, we began by optimizing the GP hyperparameters to the out-of-transit data to find the starting locations for the GP hyperparameters. The starting value for R_p/R_⋆ was taken from Hellier et al. (2012) and from visual inspection of the light curve for T₀. The chains were then initialized with a small scatter around these starting values. We ran the MCMC for 2000 steps with 300 walkers (20 × n_p, where n_p is the number of parameters) and calculated the 16th, 50th, and 84th percentiles for each parameter after discarding the first 1000 steps as burn in. Following the george documentation¹⁷ , we then ran a second chain with the walkers initiated with a small scatter around the 50th percentile values for another 2000 steps with 300 walkers and again discarded the first 1000 steps. We measure R_p/R_⋆ values of ${0.11176}_{-0.001394}^{+0.001374}$, ${0.111961}_{-0.000729}^{+0.000804}$, and ${0.113189}_{-0.000692}^{+0.000750}$, for visits 57, 58, and 59, respectively, from the white light curves (Figure 1).

Zoom In Zoom Out Reset image size

To derive the spectroscopic light curves, we binned the G430L and G750L spectra into 12 and 17 spectrophotometric channels, which are listed in the first column of Table 1. We then fit and detrended each spectrophotometric light curve following the same procedure as the white light curves, but kept T₀ fixed to the result from the white light-curve fit. We ran MCMCs to each light curve, following the same process as for the white light curves but with 280 walkers since there was one fewer fit parameter. The resulting R_p/R_⋆ values for each spectroscopic channel are presented in Table 1.¹⁸ We used the values of R_p/R_⋆ given in Table 1 to derive the optical transmission spectrum shown in Figure 2.

Zoom In Zoom Out Reset image size

Table 1.  Transmission Spectrum of WASP-62b Measured with STIS G430L and G750L and Spitzer IRAC

λ (Å)	Rp/R*
2900–3700	${0.11210}_{-0.00115}^{+0.00116}$
3700–4100	${0.11285}_{-0.00066}^{+0.00072}$
4100–4400	${0.11308}_{-0.00051}^{+0.00053}$
4400–4600	${0.11219}_{-0.00053}^{+0.00054}$
4600–4700	${0.11126}_{-0.00053}^{+0.00050}$
4700–4800	${0.11212}_{-0.00044}^{+0.00043}$
4800–4900	${0.11159}_{-0.00049}^{+0.00053}$
4900–5000	${0.11133}_{-0.00067}^{+0.00070}$
5000–5100	${0.11166}_{-0.00059}^{+0.00055}$
5100–5300	${0.11199}_{-0.00039}^{+0.00042}$
5300–5500	${0.11263}_{-0.00063}^{+0.00052}$
5500–5700	${0.11184}_{-0.00049}^{+0.00049}$
5700–5800	${0.11397}_{-0.00078}^{+0.00069}$
5800–5878	${0.11362}_{-0.00068}^{+0.00066}$
5878–5913	${0.11494}_{-0.00066}^{+0.00063}$
5913–6070	${0.11290}_{-0.00065}^{+0.00065}$
6070–6200	${0.11227}_{-0.00051}^{+0.00047}$
6200–6300	${0.11299}_{-0.00086}^{+0.00079}$
6300–6513	${0.11376}_{-0.00059}^{+0.00056}$
6513–6613	${0.11489}_{-0.00078}^{+0.00076}$
6613–6800	${0.11208}_{-0.00075}^{+0.00074}$
6800–7000	${0.11344}_{-0.00058}^{+0.00060}$
7000–7200	${0.11130}_{-0.00061}^{+0.00058}$
7200–7590	${0.11077}_{-0.00042}^{+0.00041}$
7590–7740	${0.10910}_{-0.00072}^{+0.00074}$
7740–8100	${0.11136}_{-0.00082}^{+0.00089}$
8100–8500	${0.11346}_{-0.00125}^{+0.00112}$
8500–8985	${0.11417}_{-0.00102}^{+0.00107}$
8985–10300	${0.11530}_{-0.00154}^{+0.00152}$
360000	0.10697 ± 0.00152
450000	0.10962 ± 0.00040

Download table as:  ASCIITypeset image

We fit the 3.6 μm and 4.5 μm IRAC light curves following the methods described in Alam et al. (2018). Briefly, we corrected for flux variations from intra-pixel sensitivity (e.g., Charbonneau et al. 2005, 2008; Reach et al. 2005; Knutson et al. 2008; Ingalls et al. 2012; Krick et al. 2016), and accounted for systematics by fitting a model of the functional form:

$$\begin{eqnarray}&&F(t)={c}_{0}\,+{c}_{1}x\,+{c}_{2}{x}^{2}\,+{c}_{3}y\,+{c}_{4}{y}^{2}\,+{c}_{5}{xy}\,+{c}_{6}t,\end{eqnarray} \tag{ 1 }$$

where F(t) is the stellar flux as a function of time, the coefficients c₀ through c₆ are free parameters, x and y are the stellar centroid positions on the detector, and t is time. We marginalized over all possible combinations of this model using the Gibson (2014) procedure. To measure R_p/R_⋆ for the transmission spectrum, we fixed P, a/R_⋆, and i to the values from Hellier et al. (2012) and fit for R_p/R_⋆ and T₀. The limb darkening coefficients were also fixed to their theoretical values based on 3D stellar atmosphere models. The measured R_p/R_⋆ values are included in Table 1.¹⁹ We note that our transit depth measurements are slightly lower than those previously published in Garhart et al. (2020), and this discrepancy may be due to differences in the data reduction procedure, such as choice of baseline ramp shapes (linear versus quadratic) or where to trim out-of-transit data.

Due to the nearly ubiquitous nature of condensation clouds and photochemical hazes in exoplanet atmospheres (e.g., Wakeford et al. 2019), few benchmark planets with cloud-free, haze-free atmospheres at the pressures probed by transmission spectroscopy are currently known. The detection of the Na i line wings at 0.59 μm in the atmosphere of WASP-62b marks the first space-based observation of the pressure-broadened wings of the Na D-lines, and suggests that WASP-62b possesses a clear terminator. Clear atmosphere exoplanets present an unmatched opportunity to obtain increasingly precise retrieved abundance constraints, since they are unhampered by cloud-composition degeneracies (e.g., Figure 10 of MacDonald & Madhusudhan 2017a).

We compare our results for WASP-62b with the ground-based VLT/FORS2 observations of WASP-96b, the clearest known exoplanet to date (Nikolov et al. 2018), in Figure 3. Although both WASP-62b and WASP-96b display the pressure-broadened Na line wings, WASP-62b shows no evidence of K i absorption at 0.76 μm whereas WASP-96b is consistent with weak evidence of K i and Li i (Nikolov et al. 2018). The missing or weak potassium absorption for these clear planets is contrary to atmospheric models, which predict the presence of both alkali features in clear exoplanets (Seager & Sasselov 2000). This observed trend may be due to a difference in the primordial abundances of the species, since Na is ∼15 times more abundant than K in a solar abundance atmosphere (Lodders 2003). The relative condensation temperatures of NaCl and KCl, as well as the photoionization energies of sodium and potassium, may further alter the relative amplitudes of the Na and K absorption features (Nikolov et al. 2018).

Zoom In Zoom Out Reset image size

In light of our tentative inference of SiH opacity in WASP-62b's transmission spectrum, here we consider theoretical predictions for gas-phase Si-bearing molecules at the equilibrium temperature (T_eq = 1440 ± 30 K; Hellier et al. 2012) of WASP-62b. Equilibrium chemistry expectations for silicates suggest that the dominant Si-bearing species are SiO and SiS at ∼1400 K (Visscher et al. 2010; Woitke et al. 2018). We investigate predictions for Si-bearing species in the atmosphere of WASP-62b using the planet-specific forward model ATMO grid (e.g., Tremblin et al. 2015; Goyal et al. 2020) of 1D radiative-convective equilibrium pressure–temperature profiles with corresponding self-consistent equilibrium chemistry. This grid of model atmospheres spans a range of recirculation factors (RCF), metallicities, and C/O ratios.²¹ The RCF parameterizes the redistribution of input stellar energy in the planetary atmosphere, where a value of 0.5 represents efficient redistribution and 1 corresponds to no redistribution.

Figure 4 (top panel) shows the equilibrium chemical abundances of H₂O, SiH, SiO, and SiS with uniform redistribution (RCF = 0.5), solar metallicity, and a solar C/O ratio. While SiH is tentatively inferred by our retrievals, equilibrium expectations predict that SiO and SiS should be orders of magnitude more abundant than SiH in WASP-62b's atmosphere. Since silicate condensation has also been found to increase the carbon-to-oxygen (C/O) ratio in gas from solar to super-solar values (Woitke et al. 2018), we also investigate how the abundances of these species change for sub-solar, solar, and super-solar C/O. Figure 4 (bottom panel) shows that SiS overtakes SiO to become the most abundant Si-bearing molecule for C/O  ≳ 0.7. This is driven by a decrease in the SiO abundance for super-solar C/O ratios, analogous to the well-known decrease in H₂O abundance for enhanced C/O ratios (Madhusudhan 2012). These results expand upon the silicon chemistry predictions from Visscher et al. (2010)—who considered rainout Si chemistry at a solar C/O ratio—demonstrating that SiS may be the most abundant gas-phase Si-bearing species for exoplanets with super-solar C/O ratios.

Zoom In Zoom Out Reset image size

These equilibrium predictions appear in tension with our inferred SiH abundance (∼10⁻⁴). Our retrievals do not infer SiO, despite its inclusion, due to its characteristic slope at UV-optical wavelengths (e.g. Sharp & Burrows 2007) differing from the feature at 0.4 μm we attribute to SiH. Evidence of SiS cannot currently be assessed at optical wavelengths, due to the lack of currently available line lists for low energy (bluer than ∼2 μm) transitions (Upadhyay et al. 2018). It is therefore possible that the spectral feature we attribute to SiH could be a misclassified SiS feature—optical line list data for SiS would resolve this ambiguity. However, we stress that the evidence for SiH from our STIS data remains low (2.1σ), and future observations will be required to assess the presence of gas-phase Si-bearing species in WASP-62b's atmosphere.

A clear atmosphere for WASP-62b opens the door to extreme precision molecular abundance measurements from early JWST observations. To quantitatively assess the potential of WASP-62b as a priority JWST target, we offer predictions from a simulated retrieval analysis of synthetic JWST ERS observations of WASP-62b.

We generated a moderate-resolution (R = 10,000) model transmission spectrum, including Na, H₂O, NH₃, FeH, and SiH, with the median retrieved atmospheric properties from our HST+Spitzer analysis. To explore the ability of JWST to constrain WASP-62's C/O ratio, we additionally injected CO, CO₂, and CH₄ abundances into our model. The representative values were taken to be the 10 mbar abundances of these species from the self-consistent model grid of Goyal et al. (2020), assumed uniform throughout the atmosphere. We then simulated JWST ERS observations with PandExo²² (Batalha et al. 2017)—in the Panchromatic Transmission configuration outlined by Bean et al. (2018)—spanning four transits across the following modes:

• 1.  
  NIRISS SOSS orders 1 and 2 (0.64–2.8 μm) [61 ppm/141 ppm].
• 2.  
  NIRSpec G235H (1.7–3.1 μm) [39 ppm].
• 3.  
  NIRCam F322W2 (2.4–4.1 μm) [67 ppm].
• 4.  
  NIRSpec G395H (2.9–5.2 μm) [57 ppm].

We binned the simulated PandExo observations to R = 100 for all modes, neglecting errors exceeding 1000 ppm (mainly the reddest NIRISS SOSS second-order data and near the NIRSpec detector gaps at 2.24 μm and 3.80 μm), for a total of 677 data points. The square brackets above denote the resulting mean data precision for each mode. We removed Gaussian scatter from the data set, centering the data on the actual model, to ensure our retrieval results are unbiased by a given noise instance (see Feng et al. 2018). Finally, we retrieved the synthetic data set using POSEIDON.

Our simulated JWST observations and best-fitting retrieved model spectrum are shown in Figure 5. The combined ERS instrument modes can detect many prominent optical and infrared spectral features, as well as obtain precise constraints on the planetary atmospheric metallicity (here taken as O/H) and C/O. Each JWST ERS mode covers at least one of the H₂O band heads at 0.95 μm, 1.15 μm, 1.4 μm, 1.9 μm, 2.7 μm, and 4.3 μm. NIRISS SOSS samples both the red wings of the Na D-lines and strong FeH features at 0.8 μm, 0.9 μm, and 1.0 μm. Both NIRISS SOSS and NIRSpec G235H are highly sensitive to NH₃ absorption, via the K-band feature at 2.2 μm (see also MacDonald & Madhusudhan 2017b). NIRCam F322W2 and NIRSpec G395H sample the CO, CO₂, and CH₄ features between 3 and 5 μm. NIRSpec G395H can also probe a broad SiH feature from ∼4.6 to 5.3 μm, testing if gas-phase silicon species are indeed present in WASP-62b's atmosphere. The combined observing configuration also provides the optical baseline necessary for obtaining precise molecular abundances and atmospheric metallicities.

Zoom In Zoom Out Reset image size

We predict that JWST observations of WASP-62b, within the scope of the ERS program, can conclusively detect Na (12.1σ), H₂O (35.6σ), FeH (22.5σ), SiH (6.3σ), NH₃ (11.1σ), CO (8.1σ), CO₂ (9.7σ), and CH₄ (3.6σ). The clear atmosphere offers remarkably precise abundance constraints: 0.12 dex for H₂O, 0.14 dex for FeH, 0.21 dex for Na, and 0.30 dex for SiH, 0.26 dex for CO, 0.18 dex for CO₂, and 0.20 dex for CH₄ (shown in Figure 5, lower panels). These predicted abundance constraints, from a single transit with each JWST mode, would immediately outclass the most precise abundances obtained for hot Jupiters with existing facilities to date (≲0.3 dex for H₂O, e.g., Welbanks et al. 2019).

One of the key goals of JWST is to measure the metallicity and C/O ratios of a population of exoplanets (Beichman et al. 2014; Bean et al. 2018). Our injected abundances of CO, CO₂, and CH₄ allow us to explore the expected constraints achievable by the JWST ERS program. We derive posterior distributions for the metallicity (O/H in solar units) and C/O from our individual abundance posteriors using the methodology of MacDonald & Madhusudhan (2019), as shown on the bottom row of Figure 5. We predict that it is possible to constrain WASP-62b's C/O to ±0.1 and the atmospheric metallicity to ±0.21 dex. Note that these constraints do not require the assumption of chemical equilibrium.