JWST-TST DREAMS: Quartz Clouds in the Atmosphere of WASP-17b

The ExoTiC-MIRI ²³ pipeline is interoperable with the JWST Science Calibration Pipeline (jwst; Bushouse et al. 2022), enabling default processing steps to be switched out, or interleaved, with custom steps when warranted. Starting from the uncal.fits files, we processed the data through the steps dq_init, saturation, custom_drop_groups, custom_linearity, dark_current, jump, and ramp_fit, where the custom_ prefix denotes an ExoTiC-MIRI step. For the default jwst steps, standard settings from pipeline v1.8.2 were used and CRDS context 1077, except for the jump step where the rejection threshold was increased to 15 to prevent spurious flagging of cosmic rays, and the gain value was set to 3.1 electrons per data number (Bell et al. 2023).

As for the custom steps, the 175 groups are a sufficiently large number such that we could refine the group-level ramps prior to ramp fitting. Inspection of these ramps revealed evidence of the reset switch charge decay (RSCD) effect impacting the first ∼12 groups, the last-frame effect pulling down the final group (Ressler et al. 2015; Wright et al. 2023), and residual nonlinearity remaining in the groups between. This nonlinearity persisted even after the default linearity correction was applied, most likely due to the interplay between the nonlinearity effect and debiasing-induced brighter-fatter effect (BFE) changing the apparent sensitivity of each pixel (Argyriou et al. 2023). As such, we performed a self-calibration of the linearity correction. The custom_linearity step utilizes groups 12–40 as the presumed linear portion, and then extrapolates a linear fit from this region to derive a correction factor based on deviations of groups 12–174. The correction factor takes the form of a fourth-order polynomial, where the zeroth- and first-order terms are held fixed at 0 and 1, respectively. The correction was derived per detector amplifier, essentially segmenting the correction into columns of differing fluence across the core of the point-spread function owing to the dependence of the RSCD and BFE on said fluence. The first 12 and final group of each integration were dropped altogether.

After ramp fitting, the newly created rateimages.fits files were processed with the steps assign_wcs, src_type, flat_field, custom_clean_outliers, custom_background_subtract, and custom_extract_1d. The custom cleaning step aims to remove any unidentified outliers as well as replace known bad pixels from the data quality arrays. To accomplish this, a spatial profile was estimated from polynomial fits to the detector columns as per optimal extraction (Horne 1986), and pixels were iteratively replaced by this profile value if they were more than four standard deviations discrepant or had the data quality flag do_not_use. The background step subtracts a row-by-row background using the median value from columns 12 to 22 and 50 to 68. A time series of stellar spectra were then extracted using a fixed-width box aperture. This aperture was centered on column 36 and extended 3 pixels in either direction.

The Eureka! ²⁴ reduction pipeline leverages the JWST Science Calibration Pipeline (v1.8.2, jwst; Bushouse et al. 2022) for stages 1 and 2. We started our Eureka! reduction with the uncal.fits files. In the stage 1 "ramps-to-slopes" step, the default JWST MIRI time-series observations settings are applied with the first and last groups of each integration masked, default linearity correction, 15σ jump rejection threshold, and default ramp-fitting algorithm applied. In stage 2, the ramp images created in stage 1 (rateints.fits) are further minimally calibrated according to default jwst pipeline settings. We performed stage 2 reductions both with and without the photom step applied to test for any sensitivity in our solution to wavelength-dependent calibrations (e.g., from the gain; Bell et al. 2023). The output of these first two stages are calibrated images (calints.fits) from which time-series stellar spectra can be extracted.

In stage 3, background subtraction and spectral extraction are performed on the calibrated integration images. Here, we selected a half-width aperture for the spectral extraction and background exclusion regions of 4 and 10 pixels from the central pixel of the spectral trace, respectively, which served to minimize the scatter in our extracted spectra and avoided several problematic pixels. Background subtraction was performed row by row for each integration image with outlier identification performed along the time axis with a double-iteration 5σ threshold rejection scheme. Optimal spectral extraction was then performed using a spatial profile defined by a median frame constructed from the time series. A 5σ threshold was applied in flagging outliers in the median frame and then a 10σ threshold for outlier rejection was used for both constructing the spatial profile and during the optimal spectral extraction. In stage 4, spectroscopic and white light curves were then generated from the time series of the derived stellar spectra with a 5σ, five iteration clipping of outliers to a rolling temporal median spanning 10 data points.

Light-curve fitting was performed on the white light curve (5–12 μm) and at three spectroscopic resolutions with 0.125, 0.25, and 0.5 μm-wide bins for both the ExoTiC-MIRI and Eureka! reductions. These three different spectral resolutions were chosen to test the sensitivity of our derived spectra to the "odd–even row effects" that are known to exist for the JWST MIRI detector (Ressler et al. 2015). The light curves binned at 0.5 μm are shown in Figure 1. The independent ExoTiC-MIRI and Eureka! reductions both found that the light-curve model, f(t), that best represented the astrophysical and systematic trends in the observations had the form

$$\begin{eqnarray}&&f(t)=T({t}_{{\rm{s}}},{\boldsymbol{\theta }})\times S({t}_{{\rm{s}}},{t}_{{\rm{c}}}),\end{eqnarray} \tag{ 1 }$$

where T is the physical transit model, with the parameter vector θ , and the systematics model is described by

$$\begin{eqnarray}&&S({t}_{{\rm{s}}},{t}_{{\rm{c}}})={r}_{0}\exp ({r}_{1}{t}_{{\rm{s}}})\times ({c}_{0}+{c}_{1}{t}_{{\rm{c}}}).\end{eqnarray} \tag{ 2 }$$

Here, r₀, r₁, c₀, and c₁ are constants to be fit, and t_s and t_c are the observation times minus the light-curve start time and center-of-transit time, respectively.

Zoom In Zoom Out Reset image size

Both the ExoTiC-MIRI and Eureka! pipelines employed batman (Kreidberg 2015) to generate the transit model, T. The ExoTiC-MIRI light-curve fits used a fixed quadratic limb-darkening law where the coefficients for each bin are computed using ExoTiC-LD (Grant & Wakeford 2022) and Set One of the MPS-ATLAS stellar models (Kostogryz et al. 2022, 2023). The Eureka! light-curve fits also used quadratic limb-darkening but reparameterized according to Kipping (2013), and additionally these fits decorrelated against the measured position and width of the spectral trace. Both the ExoTiC-MIRI and Eureka! light-curve fits also include an error multiplier, β, to capture possible error underestimation/overestimation, in particular due to the uncertainties in the detector gain (Bell et al. 2023). The ExoTiC-MIRI and Eureka! light-curve fits cut off the first 100 (47 minutes) and 65 (30 minutes) integrations, respectively, to alleviate the model from fitting the worst of the systematic ramp. Additionally, the ExoTiC-MIRI light-curve fits manually masked the first data point from each segment, 12 in total, due to anomalous background levels, along with any remaining 4σ outliers from a running median. These same 12 integrations are flagged and removed in stage 4 of the Eureka! pipeline.

The fitting was performed on both the ExoTiC-MIRI and Eureka!-generated light curves with the Markov Chain Monte Carlo (MCMC) algorithm emcee (Foreman-Mackey et al. 2013). In both the ExoTiC-MIRI and Eureka! light-curve fitting, the white light curve was first fit to constrain the best-fit values for the center-of-transit time. Additionally, in the ExoTiC-MIRI light-curve fitting the semimajor axis and inclination were estimated from the white light curve, while Eureka! assumed the values from Alderson et al. (2022). These system parameter values were then held fixed in the spectroscopic light-curve fitting; see Table 1 for a complete list of system parameters and prior distributions used in the light-curve fitting. The transmission spectra (${\left[{R}_{p}/{R}_{* }\right]}^{2}$) derived from each reduction are displayed in Figure 2.

Table 1. Light-curve Fitting Parameter Information

	ExoTiC-MIRI		Eureka!
Parameter	Prior	Value	Prior	Value
P (days)	Fixed (Alderson et al. 2022)	3.73548546	Fixed (Alderson et al. 2022)	3.73548546
t0 (BJDTDB − 2460016.7)	${ \mathcal U }$(−0.1, 0.1)	0.026452 ± 0.000061	${ \mathcal N }$(0.02648, 0.05)	0.026477 ± 0.000060
a/R*	${ \mathcal U }$(6, 8)	7.110 ± 0.040	Fixed (Sedaghati et al. 2016)	7.025
i (degrees)	${ \mathcal U }$(80, 90)	87.217 ± 0.135	Fixed (Alderson et al. 2022)	86.9
Rp /R*	${ \mathcal U }$(0.10, 0.15)	0.12472 ± 0.00016	${ \mathcal N }$(0.123, 0.05)	0.12506 ± 0.00019
Teff (K)	Fixed (Southworth 2012)	6550	⋯	⋯
Fe/H	Fixed (Southworth 2012)	−0.25	⋯	⋯
$\mathrm{log}g$ (cm s−2)	Fixed (Southworth 2012)	4.149	⋯	⋯
u1	Fixed	0.05	⋯	⋯
u2	Fixed	0.06	⋯	⋯
q1	⋯	⋯	${ \mathcal U }$(0, 0.05)	0.0067 ± 0.0016
q2	⋯	⋯	${ \mathcal U }$(0.15, 0.35)	0.258 ± 0.069
c0	${ \mathcal N }$(1, 0.01)	1.000091 ± 0.000035	${ \mathcal N }$(1.001, 0.01)	1.006164 ± 0.000055
c1	${ \mathcal N }$(0, 0.01)	−0.00033 ± 0.00026	${ \mathcal N }$(0, 0.01)	−0.00102 ± 0.00036
r0	${ \mathcal N }$(0, 0.1)	0.00135 ± 0.00011	${ \mathcal N }$(0, 0.01)	0.00151 ± 0.00022
r1	${ \mathcal U }$(0, 200)	40.3 ± 8.0	${ \mathcal U }$(5, 150)	27.5 ± 5.1
β	${ \mathcal U }$(0.1, 5)	1.16 ± 0.02	${ \mathcal N }$(1.5, 0.5)	1.724 ± 0.035
x0	⋯	⋯	${ \mathcal N }$(0, 0.1)	0.0439 ± 0.0069
x1	⋯	⋯	${ \mathcal N }$(0, 0.5)	0.026 ± 0.020

Notes. Values are shown for the white light-curve fits. Parameter definitions: orbital period, P; time of transit center, t₀; semimajor axis in units of stellar radius, a/R_*; inclination, i; planet radius in units of stellar radius, R_p /R_*; stellar effective temperature, T_eff; stellar metallicity, Fe/H; stellar gravity, $\mathrm{log}g;$ quadratic limb-darkening coefficients, u₁ and u₂; Kipping (2013) reparameterization of the quadratic limb-darkening coefficients, q₁ and q₂; systematic model parameters as defined in Equation (2), c₀, c₁, r₀, r₁; error multiplier, β; decorrelation coefficients for trace position and width, x₀ and x₁. Units for logarithmic parameters refer to the argument.

Download table as:  ASCIITypeset image

Our analysis with both the ExoTiC-MIRI and Eureka! pipelines produces white light curves with precision (standard deviation of the normalized residuals, or SDNR) measured at 465 and 453 ppm, respectively. The slight (12 ppm) precision increase in the Eureka! light curves is driven by two choices in the reductions. First, the ExoTiC-MIRI reduction discarded the first 12 groups from the start of each ramp, decreasing the ramp fitting's statistical power with the aim of increasing accuracy, whereas the Eureka! reduction only discarded the first group. Second, the Eureka! reduction utilized optimal extraction rather than a box aperture, and this was found to improve the precision.

The transmission spectra from the ExoTiC-MIRI and Eureka! pipelines exhibit a consistent morphological structure, albeit with one primary difference: At wavelengths <8 μm, the Eureka! spectrum is systematically deeper. Further investigation showed the different linearity corrections lead to this discrepancy. The derived custom linearity correction in ExoTiC-MIRI (see Section 2.1) is smaller than the default correction for the same increase in observed data numbers, and this correction is amplifier dependent. An analysis of the group-level ramps revealed the default linearity correction led to larger deviations from a linear ramp relative to the custom linearity correction. This is true across all detector columns but, critically, these deviations were largest (by 58%) on the dispersion axis where most of the flux is collected. These differences are most likely due to ExoTiC-MIRI's custom correction better taking into account the debiasing-induced BFE (Argyriou et al. 2023) specific to these data, which predominately impacts the bright pixels at the center of the dispersion axis. This analysis indicates that the ExoTiC-MIRI reduction is more reliable, and we choose to focus our atmospheric modeling primarily on these data.

We note that at wavelengths >8 μm the ExoTiC-MIRI and Eureka! transit depths are consistent to better than one standard deviation across all three resolutions, and the residuals between the two reductions pass the Shapiro–Wilk test for normality (Shapiro & Wilk 1965). Any small stochastic differences arise from the choices to fit or fix the limb-darkening coefficients, semimajor axis, inclination, and the type of extraction aperture. We found the transmission spectra were insensitive to decisions about the limb-darkening law within each reduction.

We used a planet-specific grid of self-consistent model atmospheres for WASP-17b generated using ATMO (Amundsen et al. 2014; Drummond et al. 2016; Tremblin et al. 2016; Goyal et al. 2018), with radiative-convective equilibrium pressure–temperature (P–T) profiles consistent with equilibrium chemistry (Goyal et al. 2020). This grid was generated for a range of heat redistribution factors (0.25, 0.5, 0.75, 1.0), metallicities (0.1×, 1×, 10×, 50×, 100×, 200× solar) and C/O ratios (0.35, 0.55, 0.70, 0.75, 1.0, 1.5). The internal temperature of the planet is assumed to be 100 K, although this value is not well known (see, e.g., Thorngren et al. 2019; Sarkis et al. 2021). We test internal temperatures of 100, 200, and 300 K and find the differences in the modeled transmission spectrum to be negligible. Using this grid of WASP-17b model atmospheres, a grid of simulated transmission spectra was generated at R ∼ 1000 for a range of Rayleigh scattering haze factors (1×, 10× nominal Rayleigh scattering) and gray cloud factors (0.0×, 0.5×, 1.0×, 5.0× H₂ Rayleigh scattering cross section at 350 nm).

We find that the best-fit ATMO model to the combined HST + Spitzer and JWST MIRI LRS data sets has a redistribution factor of 0.5, solar metallicity, supersolar C/O ratio (0.7), haze factor of 10, and gray cloud factor of 0.5. This model is shown in the top panels of Figure 3. We note that models with 0.1× to 50× solar metallicity and 0.35 to 0.7 C/O ratios lie within the 3σ range of the best-fit model. We found consistent results across all three binning regimes for both the ExoTiC-MIRI and the Eureka! reductions presented in Figure 2, but note that the best-fit offsets between the JWST MIRI LRS and HST + Spitzer data differed between ExoTiC-MIRI and Eureka! pipeline reductions. As previously discussed, the different linearity corrections applied in the ExoTiC-MIRI and Eureka! reductions resulted in offsets in the transit depth shortward of 8 μm. However, we find that allowing for an offset between the JWST MIRI LRS data and the HST + Spitzer data results in identical best-fitting atmospheric models, which indicates that our atmospheric interpretation is minimally sensitive to differences between the ExoTiC-MIRI and Eureka! reductions or the choice of spectral bin size.

We computed radiative-convective thermochemical equilibrium (RCTE) atmospheric models for WASP-17b using the well-vetted open-source model PICASO v3.1 ²⁵ (Batalha et al. 2019; Mukherjee et al. 2023), which has heritage from Fortran codes developed to study solar system giant planets (e.g., Marley & McKay 1999) and brown dwarfs (e.g., Marley & McKay 1999). We computed a grid of cloud-free models as a function of interior temperature of the planet (200 and 300 K; Thorngren et al. 2019; Sarkis et al. 2021), atmospheric metallicity (nine values between 1× and 100× solar), C/O ratio (five values between 0.25× and 2× solar), and the heat redistribution factor (0.5, 0.6, 0.7, 0.8). PICASO's RCTE module utilizes correlated-k opacities that are detailed in Marley et al. (2021) and released by Lupu et al. (2021). We included the opacity sources for 29 species, but the most important for the JWST MIRI LRS wavelength region is the line list of H₂O (Polyansky et al. 2018). The chemical equilibrium abundances are computed on a pressure–temperature-M/H-C/O grid of thermochemical equilibrium models presented in Marley et al. (2021) following the work of Gordon & McBride (1994), Fegley & Lodders (1994), Lodders (1999, 2002), Lodders & Fegley (2002), Visscher et al. (2006), and Visscher et al. (2010), using elemental abundances from Lodders (2010). From our WASP-17b climate models, we computed transmission spectra using opacities resampled to R = 60,000 (Batalha et al. 2020) from original R ∼ 10⁶ line-by-line calculations detailed in Freedman et al. (2008) and Gharib-Nezhad et al. (2021). We fit the cloud-free grid to both the ExoTiC-MIRI and Eureka! JWST MIRI LRS reductions using the "MLFriends" nested sampling algorithm (Buchner 2016, 2019) implemented in the open-source Ultranest code (Buchner 2021), and find agreement between the two, regardless of binning scheme. The best-fit model from the cloud-free grid has an internal temperature of 200 K, redistribution factor of 0.8, metallicity of 100× solar, and supersolar C/O ratio (0.7).

In order to account for the presence of clouds in WASP-17b's atmosphere, we used the cloud-free grid's pressure–temperature profiles to compute condensation cloud profiles using the Virga cloud model ²⁶ (Batalha et al. 2020; Rooney et al. 2022). The cloud methodology of Virga is detailed in Ackerman & Marley (2001). Virga models the balance between the turbulent diffusion (K_zz ) and sedimentation (f_sed) in horizontally uniform cloud decks, and therefore requires these two additional model parameters. We include two condensable species, SiO₂(s) and Al₂O₃(s), given the range of temperatures expected in WASP-17b's atmosphere and presence of a clear spectroscopic feature near 8.6 μm. For SiO₂(s), we used the α-crystal optical constants computed at 928 K by Zeidler et al. (2013). For Al₂O₃(s), we used the amorphous optical constants computed at 873 K by Koike et al. (1995) and Begemann et al. (1997). We note that we did explore a number of other cloud species, including MgSiO₃(s) and Mg₂SiO₄(s), for which we show the cross sections in the left-hand panel of Figure A1, within the Appendix. However, these cloud species could not match the observed mid-infrared spectroscopic features. We also tested the amorphous form of SiO₂(s), but this results in a worse match to the 8.6 μm feature, and the reduced χ² of the fit increases from 0.98 to 1.05 (90 degrees of freedom).

The saturation vapor pressure curve of Al₂O₃ was calculated in Wakeford et al. (2016). For SiO₂, we calculated the vapor pressure and condensation curve for the net thermochemical reaction, SiO(g) + H2O(g) = SiO₂(s) + H2(g), which is adopted because SiO and H₂O are the dominant Si- and O-bearing gases, respectively, at the temperatures pertinent to WASP-17b. This leads to condensation curves given by

$$\begin{eqnarray}&&\mathrm{log}{P}_{T}\approx 13.168\,-\,28265/T\,-\,[\mathrm{Fe}/{\rm{H}}],\end{eqnarray} \tag{ 3 }$$

or, equivalently,

$$\begin{eqnarray}&&{10}^{4}/{T}_{\mathrm{cond}}\approx 6.14\,-\,0.35\mathrm{log}{P}_{T}\,-\,0.70[\mathrm{Fe}/{\rm{H}}],\end{eqnarray} \tag{ 4 }$$

where P_T is in bars. These condensation curves are shown in the right-hand panels of Figure 3.

We fit our cloudy WASP-17b spectra to the ExoTiC-MIRI reduction combined with HST+Spitzer data from Alderson et al. (2022) using Ultranest (Buchner 2021). For the parameters confined to the grid, we obtain 1σ constraints of 30–100× solar metallicity, 0.4–0.7 C/O, a heat redistribution of 0.6–0.7, and internal temperature of 200–300 K. For the two cloud parameters, the 1σ constraints are $\mathrm{log}{K}_{{zz}}={9}_{-0.29}^{+0.75}$ cm² s⁻¹ and ${f}_{\mathrm{sed}}={0.1}_{-0.76}^{+0.22}$. We compared the cloudy model with the cloud-free model using their likelihood ratio, and we find that the SiO₂(s) cloud model is preferred by 4.2σ. The maximum-likelihood cloudy model is shown in the bottom panels of Figure 3. Additionally, in Figure 4 the opacity contributions to this model are shown, highlighting that H₂O vapor and SiO₂(s) clouds dominate WASP-17b's transmission spectrum at JWST MIRI LRS wavelengths.

We conducted free-chemistry retrievals with the inclusion of compositionally specific aerosols on HST, Spitzer, and JWST MIRI LRS observations of the transmission spectrum of WASP-17b using the open-source atmospheric retrieval code POSEIDON ²⁷ (MacDonald & Madhusudhan 2017; MacDonald 2023). Here, we build upon the existing cloud models in POSEIDON to include retrievals with Mie-scattering aerosols. We calculated a database of effective extinction cross sections of aerosols species by adapting Zhang et al.'s (2019) implementation of the Mie-scattering algorithm presented in Kitzmann & Heng (2018). The effective aerosol extinction cross section is the combined absorption and scattering cross section integrated over a log-normal radius distribution centered around a mean particle size, r_m (in microns). For computational efficiency, we precomputed effective extinction cross sections from the refractive-index databases of Wakeford & Sing (2015) and Kitzmann & Heng (2018), where the relevant species are based on work by Henning & Mutschke (1997), Palik (1998), Andersen et al. (2006), and Zeidler et al. (2013). We considered mean particle sizes ranging from 0.001 to 10 μm and wavelengths spanning 0.2–30 μm at R = 1000.

Our POSEIDON retrievals adopt a model configuration similar to Alderson et al. (2022), but with the addition of a parameterized Mie-scattering cloud model. We assume a one-dimensional H₂–He-dominated atmosphere (with He/H₂ =0.17) with the gas phase containing free abundances of H₂O, CH₄, CO₂, CO, Na, and K. The specific line lists used are detailed in the Appendix of MacDonald & Lewis (2022), with particular reference to Polyansky et al. (2018), Yurchenko et al. (2017), Chubb et al. (2021), Tashkun & Perevalov (2011), Li et al. (2015), and Barklem & Collet (2016). We assume an isothermal pressure–temperature profile, which follows from the modeling of Alderson et al. (2022). The reference pressure is set at 10 bar. We ran retrievals both with and without the JWST MIRI LRS data (0.25 μm ExoTiC-MIRI reduction), but always including the HST+Spitzer data. We consider three models: (i) a cloud-free model, (ii) the gray cloud+scattering haze prescription found in Alderson et al. (2022), and (iii) Mie-scattering SiO₂(s) clouds (both crystalline and amorphous states considered separately). Aerosol clouds are parameterized by the mean particle size, r_m, the cloud-top pressure, P_cloud, the width of the cloud in log pressure space, ${\rm{\Delta }}\mathrm{log}$ P, and the constant log mixing ratio of the aerosol in the cloud. The priors for our POSEIDON retrievals are summarized in Table 3.

We computed the marginalized likelihood for each of the models and compared them via the Bayes factor to assess the support for each model. We find that the SiO₂(s) clouds model is preferred over no clouds by 3.5σ and also preferred over the generic aerosol prescription by 2.6σ. The best-fit model finds the SiO₂(s) clouds are composed of small particles (∼0.01 μm) with a low mixing ratio ($\mathrm{log}X\sim -12$) located high in the atmosphere (P_cloud ≲ 1 mbar). Our data are highly informative of the particle size, as changes to this parameter dramatically alter both the optical scattering slope and the amplitude of the 8.6 μm vibrational-mode feature (see Figure A2 within the Appendix). Both crystalline and amorphous SiO₂(s) aerosols are able to produce good fits to the data, with the crystalline form producing the best match to the peak position of the feature near 8.6 μm. In addition to SiO₂(s), we also considered Fe₂O₃(s) aerosols based on the alignment of refractive indices with the observed feature near 8.6 μm (see Figure A1 within the Appendix). However, the best-fit model requires particle sizes of ∼0.001 μm, which is physically implausible (Sridevi et al. 2023). Thus, we rule out Fe₂O₃(s) aerosols.

For completeness, we also retrieve on our 0.25 μm Eureka! reduction. These retrievals also yield a preference for clouds, with both the SiO₂(s) and generic clouds preferred at 2.6σ and 3.0σ, respectively, over the cloud-free model. While the SiO₂(s) model produces a better fit to the Eureka! data than the generic clouds (reduced χ² = 1.14 versus 1.17), the Bayes factor does not support the SiO₂(s) model given the additional parameters in its definition. However, as detailed in Section 2.4, the Eureka! reduction may be less reliable than the ExoTiC-MIRI reduction at shorter wavelengths, and this impacts the amplitude and detectability of the 8.6 μm feature.

Retrievals that do not include JWST MIRI LRS data indicate no preference over either of our cloud models, a result consistent with predictions from Mai & Line (2019). In all cases, our retrieved abundances for the gas-phase chemical species considered are consistent with the values presented in Alderson et al. (2022). In Figure 5, we show retrieved spectra with 1σ contours and posterior distributions for gas-phase species (H₂O, CH₄, CO₂, Na, K, CO), SiO₂(s) cloud properties, and isothermal atmospheric temperature.

To complement our POSEIDON retrievals, we performed free retrievals using petitRADTRANS ²⁸ (pRT; Mollière et al. 2019), which has been used previously to explore specific aerosols and chemistry in brown dwarfs and directly imaged exoplanet atmospheres (e.g., Mollière et al. 2020). Similar to POSEIDON, we opt to use an isothermal prescription for the pressure–temperature profile and consider the case of a H₂–He-dominated atmosphere with H₂O, CO₂, SO₂, CH₄, CO, Na, and K as possible gas-phase species. The specific line lists used are detailed in Mollière et al. (2019), with particular reference to Rothman et al. (2010), Yurchenko et al. (2017), Chubb et al. (2021), Underwood et al. (2016), Tóbiás et al. (2018), and Piskunov et al. (1995). The reference pressure is set at 100 bar. We also allow for the presence of crystalline SiO₂(s) and use the same refractive indices as POSEIDON. However, we leverage the cloud model implementation built into pRT, which is based on Ackerman & Marley (2001). This cloud model includes parameters for vertical mixing, K_zz , sedimentation efficiency, f_sed, and cloud-base pressure, P_base, that give rise to a population of cloud particles with a log-normal size distribution and a specific mixing ratio across a range of pressures above the cloud base. The priors for our pRT retrievals are also summarized in Table 3.

We performed pRT retrievals on the JWST MIRI LRS data alone (0.25 μm ExoTiC-MIRI reduction) and in combination with the HST+Spitzer data presented in Alderson et al. (2022). The results are displayed in Figure 5 alongside those from POSEIDON. We find the abundances for the gas-phase species considered are consistent with estimates from POSEIDON, noting that only H₂O is constrained from the JWST MIRI LRS data alone and there is no strong evidence for SO₂, a photochemically produced species that has absorption signatures at MIRI wavelengths (Tsai et al. 2023). In our pRT retrievals, we find that the inclusion of SiO₂(s) clouds is preferred over a generic-clouds model, with χ² values of 100.0 and 129.9 (80 and 82 degrees of freedom), respectively. The SiO₂(s) clouds have a cloud base in the 1–10 mbar region of WASP-17b's limb, and an average particle size of ∼0.01 μm, which is consistent with the findings from POSEIDON. These clouds are produced within the pRT cloud model with K_zz ∼ 10⁶ cm² s⁻¹ and an f_sed on the order of 5, which is consistent with the best-fit models from PICASO+Virga considering the trade-offs between these two parameters in the context of the Ackerman & Marley (2001) cloud model.

Our suite of forward models and retrievals for WASP-17b have produced evidence for an atmosphere shaped by complex chemistry, where aerosols are required to fit the transmission spectrum of WASP-17b. Across all of our modeling efforts, we have trialed aerosols composed of SiO₂(s), Al₂O₃(s), Fe₂O₃(s), MgSiO₃(s), and Mg₂SiO₄(s). We find that SiO₂(s) is favored as the dominant cloud species, and is statistically preferred over no clouds by 3.5σ in our POSEIDON free retrievals and at 4.2σ in our PICASO equilibrium-chemistry forward models. We additionally investigated a combination of SiO₂(s) and generic aerosol parameterizations with POSEIDON, and find that SiO₂(s) is still required to fit the spectrum at the 2.6σ level. This modeling suggests that small-particle SiO₂(s) clouds are present at the limb of WASP-17b's atmosphere.

Our free-chemistry retrievals all find an overall subsolar abundance of H₂O, while our equilibrium-chemistry forward models have preferences for supersolar metallicities and C/O ratios. These differences are most likely caused by the constraints imposed by RCTE in the forward models. In particular, the PICASO+Virga cloudy model has its metallicity driven to higher values because in Virga the gas mean mass mixing ratio is fixed to a constant value. Therefore, higher-metallicity solutions may be needed to compensate for the low cloud opacity. The difference between the supersolar metallicity derived from both equilibrium-chemistry forward models and the subsolar metallicities found by the free-chemistry retrievals suggest that disequilibrium processes could be at play in WASP-17b's atmosphere.

All of our atmospheric models demonstrate how the mid-infrared observations are key to identifying the cloud species, while observations at optical wavelengths may act to further constrain the cloud's average grain size, density, and vertical location. Overall, our "four-pronged" theoretical exploration of WASP-17b's atmospheric thermal structure and composition all point to a planetary transmission spectrum strongly shaped by H₂O vapor and SiO₂(s) clouds, with suggestions of absorption from CO and CO₂ that will be more fully explored with our JWST Near Infrared Imager and Slitless Spectrograph and Near Infrared Spectrograph observations of this planet (GTO-1353).