The Hubble PanCET Program: A Featureless Transmission Spectrum for WASP-29b and Evidence of Enhanced Atmospheric Metallicity on WASP-80b

Three transits of each system were observed with the Space Telescope Imaging Spectrograph (STIS) instrument in 2017 as part of the HST Panchromatic Comparative Exoplanetary Treasury Program (GO-14767; PIs: David Sing & Mercedes López-Morales). Two visits were conducted using the G430L grating (289–570 nm), and one visit was carried out with the G750L grating (526–1025 nm). The HST/STIS observations utilized the $52^{\prime\prime}\times 2^{\prime\prime}$ slit and included wavelength calibration and flat-field exposures that were scheduled during the last spacecraft orbit of each visit. The per-exposure integration time was $t_{\rm int}=250$ s across all visits. Most visits consisted of five HST orbits and contained 48 individual exposures, with the exception of the second STIS G430L observation of WASP-80, for which only four HST orbits were allocated, yielding 38 exposures. We discarded the first of the two pre-ingress orbits during the light-curve extraction process (Section 3.1).

We observed one transit of each system in 2016 using the G141 grism (1.0–1.7 $\mu$m) on the Wide Field Camera 3 (WFC3) instrument as part of HST program GO-14260 (PI: Drake Deming). Both visits utilized the now-standard spatial scan mode in order to maximize the total flux from the target on the 266$\times$266 pixel subarray. For WASP-29, 74 exposures were obtained across five HST orbits, each with a total exposure time of 112 s. For WASP-80, the integration time was set to 103 s, and the 74 exposures were scheduled across four HST orbits. While the WASP-29 observations were conducted using only the forward scan direction, both forward and backward scans were included for WASP-80 in order to reduce overhead. A separate direct image of the target using the F139M grism was scheduled at the beginning of each visit, which was used to compute the wavelength solution of each exposure (Section 3.2).

A pair of transit observations in the 3.6 and 4.5 $\mu$m bandpasses of the Infrared Array Camera (IRAC) were obtained for each system. All of the Spitzer/IRAC observations utilized the 32$\times$32 pixel subarray mode. For WASP-29, the transit observations were conducted in 2017 as part of Program 13044 (PI: Drake Deming) and had an effective per-exposure integration time of 1.92 s. Short $\sim$30 minute peak-up observations were scheduled immediately before the start of the science exposures, to allow for the telescope pointing to stabilize. The 3.6 and 4.5 $\mu$m WASP-80b transit observations were carried out in 2013 as part of Program 90159 (PI: Amaury Triaud). For this program, the integration times in the two bandpasses were set to 0.36 and 1.92 s, respectively, and no peak-up observations were included.

Two 3.6 $\mu$m and one 4.5 $\mu$m occultation observations of WASP-29 were obtained in 2010 and 2011 (Programs 60003 and 70084; PI: Joseph Harrington). These visits used an integration time of 0.36 s and included no peak-up exposures. An additional 4.5 $\mu$m secondary eclipse observation with peak-up and $t_{\rm int}=1.92$ s was carried out in 2014 (Program 10054: PI: Heather Knutson). For WASP-80, two secondary eclipses were observed in each of the two IRAC bandpasses in 2013 as part of the same Spitzer program that provided the transit observations for the system. An additional pair of 3.6 $\mu$m occultation observations with peak-up and $t_{\rm int}=0.36$ s were conducted as part of Program 10054.

The Level 1 data products, which include dark-subtracted, flat-fielded, linearized, and flux-calibrated Spitzer/IRAC data processed by the official IRAC pipeline, were downloaded from the Spitzer Heritage Archive (Spitzer Science Center, 2020).

An exoplanet’s transmission spectrum can be affected by chromospheric activity on the host star; in particular, the slope of the transmission spectrum in the optical and near-infrared is susceptible to biases from unocculted starspots or faculae (e.g., Pont et al., 2008; Sing et al., 2011; Bruno et al., 2018, 2020; Rackham et al., 2018, 2019). The former can mimic the negative slope caused by Rayleigh scattering, while the latter tend to have the opposite effect. Long-term photometric monitoring is an important tool for probing the level of stellar variability and assessing the extent to which any stellar surface heterogeneities may affect the interpretation of the transmission spectra.

We regularly observed WASP-80 over 16 semesters from 2013 to 2021 using the Tennessee State University Celestron 14 inch Automated Imaging Telescope (AIT) at Fairborn Observatory in Arizona, USA. The AIT uses an SBIG STL-1001E CCD camera with a 1024$\times$1024 Kodak KAF-1001E CCD detector ($1\farcs 2$ pixel scale). Images were obtained through a Cousins $R$-band filter. The original SBIG camera failed during the 2017B semester and was replaced with an identical camera and detector set; any observations from 2017B were not included in our photometric analysis. We calculated differential magnitudes of WASP-80 using the mean brightnesses of five nearby constant companion stars across 5–10 consecutive exposures. For more details about the data processing and analysis methodology, see Sing et al. (2015) and Wong et al. (2020a).

After removing outliers due to poor observing conditions and observations that coincided with a transit, we have a total of 596 nightly photometric measurements of WASP-80, with a typical single-measurement precision of 3–5 mmag on nights with good observing conditions. To search for stellar variability, we computed the Lomb–Scargle periodogram of each semester’s photometric time series. Data from six of the 16 semesters (2014B, 2016A, 2016B, 2018A, 2018B, and 2019A) showed statistically significant peaks in the periodogram. The most prominent peaks corresponded to variability periods ranging from 20.52 to 25.58 days across the six semesters, with a mean photometric period of 23.45 days. The peak-to-peak brightness modulation of these periodic signals ranged from 5.8 to 13.2 mmag. Figure 1 shows the differential photometry of WASP-80 during the six semesters for which periodicity was detected; the light curves are phase-folded on the respective best-fit variability periods.

From their detailed statistical study of stellar activity in the Kepler sample, Rackham et al. (2019) found that K6V stars such as WASP-80 have a $1\sigma$ range of peak-to-peak relative flux amplitudes spanning 0.35%–1.28%. Our measured $R$-band brightness modulation amplitudes for WASP-80 correspond to relative flux amplitudes of 0.54%–1.22%, entirely consistent with the reported range. Given the star’s measured sky-projected rotational velocity and radius ($v\sin i_{*}=1.27^{+0.14}_{-0.17}\,{\rm km}\,{\rm s}^{-1}$, $R_{*}=0.586^{+0.017}_{-0.018}\,R_{\odot}$; Triaud et al. 2015), the predicted rotation period is $P_{\rm rot}/\sin i_{*}=23.3\pm 3.2$ days, which is consistent with the average variability period we obtained from our photometric monitoring. We therefore attribute the source of the detected stellar activity to rotational brightness modulation from starspots, an interpretation that is also supported by the observed X-ray activity level (Section 2.4).

WASP-29 is located in the southern sky and is therefore inaccessible from the AIT. To place constraints on the level of stellar variability, we utilized the light curve provided by the Transiting Exoplanet Survey Satellite (TESS), which observed the target during Sector 2 (2018 August 22 to September 20) and Sector 29 (2020 August 26 to September 22). We considered the presearch conditioning simple aperture photometry (PDCSAP) light curve, which has been corrected by the official Science Processing and Operations Center (SPOC) pipeline (Jenkins et al., 2016) for common-mode systematics using cotrending basis vectors empirically derived from stars across the same detector array (Stumpe et al., 2014). After trimming the planetary transits and removing $3\sigma$ outliers using a 16-point-wide moving median filter, we evaluated the scatter and photometric uncertainties at various binning intervals.

At the native 2 minute cadence, the scatter and median photometric error are 1150 and 1190 ppm for Sector 2 and 1110 and 1150 ppm for Sector 29, respectively. It follows that there is no significant short-term variability on minute-long timescales. Binning at one-hour intervals, we found that the binned scatter and median of the weighted mean photometric errors are 284 and 156 ppm for Sector 2 and 210 and 151 ppm for Sector 29, indicating a low level of time-correlated noise on timescales that are comparable to the transit duration. However, it is important to acknowledge the possible presence of uncorrected instrumental systematics trends, which have been shown to persist in the PDCSAP light curves, particularly in the vicinity of scheduled interruptions due to data downlinks and reaction wheel momentum dumps (see, for example, Wong et al. 2020b, 2021). Visual inspection of the TESS WASP-29 light curve revealed flux ramps and discernible brightness discontinuities at several of the momentum dumps. Meanwhile, the Lomb–Scargle periodogram of each sector’s light curve did not show any periodicities above $3\sigma$ significance. Our variability analysis of WASP-29 is insensitive to longer-term stellar variability ($>20$ d), given that each TESS sector spans less than a month. Therefore, while we clearly do not find brightness variations above a few 100s of ppm on timescales ranging from hours to days, we cannot confidently exclude the possibility that some lower-level stellar variability may exist. These results corroborate previous ultraviolet observations of WASP-29 that revealed a magnetically quiet star at short wavelengths (dos Santos et al., 2021).

In order to further assess the activity level of the two host stars, we used archival data from the XMM-Newton observatory. WASP-80 was observed on 2014 May 13 and 2015 May 15 for a total of 49,000 s (formerly analyzed by Salz et al. 2015 and King et al. 2018). We reduced the data following standard procedures, removing intervals with high background levels. The filtered data from the two campaigns were combined, resulting in total exposure times of 39,300, 46,500, and 47,900 s in the pn, MOS1, and MOS2 detectors of the European Photon Imaging Camera (EPIC), respectively. We simultaneously fitted the spectra on the three detectors, using a two-temperature APED (Astrophysics Plasma Emission Database; Foster et al., 2012) model with $\log\,T_{1,2}=6.47\pm 0.08,6.96\pm 0.08$ (in units of Kelvin), $\log\,{\rm EM}_{1,2}=49.83^{+0.15}_{-0.21},49.70^{+0.15}_{-0.25}$ (in units of cm^-3), [Fe/H]$=-0.14$, and an interstellar medium absorption of $\log\,N_{\rm H}=18.6$ (cm^-3), at a distance of 49.716 pc (Gaia EDR3).

The spectral fit yielded an X-ray luminosity of $L_{\rm X}=3.24\times 10^{27}$ erg s^-1, with a signal-to-noise ratio of 15.6 in the spectral range 5–100 Å. We obtained a ratio of X-ray to bolometric luminosity of $\log L_{\rm X}/L_{\rm bol}=-4.9$, which indicates that WASP-80 has a moderate level of activity. A small flare was observed at X-ray wavelengths during the 2015 campaign, although no ultraviolet counterpart of the flare was detected by the optical monitor on board XMM-Newton. If we utilize the X-ray luminosity versus rotation period relation in Wright et al. (2011), a rotation period of 23.45 days (Section 2.3) would correspond to $L_{\rm X}=6\times 10^{27}$, consistent with the observed $L_{\rm X}$ value considering the level of variance in this relation.

Meanwhile, data for WASP-29 from XMM-Newton were presented and analyzed in dos Santos et al. (2021), where only an upper limit of $\log L_{\rm X}/L_{\rm bol}<-6.0$ could be derived, clearly indicating a significantly lower activity level than WASP-80.

After downloading the raw images from the Mikulski Archive for Space Telescopes, we used the latest version of CALSTIS to flat-field them. We then replaced $4\sigma$ outlier pixel values in each individual frame with the corresponding median pixel values across all frames. The column-added 1D spectra were extracted from the background-subtracted frames using the wavelength solution provided by the accompanying calibrated *sx1.fits files. The width of the extraction aperture was chosen so as to minimize the scatter in the resultant residuals from the transit light-curve fit (e.g., Deming et al., 2013). The optimal aperture widths for each visit are listed in Table 2. For data obtained with the G750L grating, we corrected for the fringing effect using the fringe flat frame that was scheduled at the end of each visit (see, for example, Nikolov et al., 2014, 2015; Sing et al., 2016).

To correct for small subpixel wavelength shifts across each visit, we set the first spectrum in the time series as the template and calculated the best-fit wavelength offsets and amplitude scaling factors using least-squares minimization to align each subsequent spectrum with the template. The normalized time series of best-fit amplitude scaling factors provided the broadband light curve for each visit. To extract the spectroscopic light curves, we used bins 100 pixels wide across the entire wavelength range spanned by the grating. To probe for the presence of alkali (Na and K) absorption lines, we included additional narrow wavelength bins at 588–591 and 766–772 nm for the G750L datasets. For these and all other light curves in our analysis, the time stamp of each spectrum was set to the mid-exposure time, converted to the BJD_TDB time standard (Eastman et al., 2010)

Following standard procedure, we discarded the first orbit prior to fitting. In addition, we removed the first data point from each orbit. Any remaining outliers in the light curves were removed by fitting each light curve separately and discarding points that were situated more than $5\sigma$ from the best-fit model. The 20th exposure was removed from the first STIS G430L visit of WASP-29. For the first STIS G430L transit observation of WASP-80, we discarded the 18th, 20th, and 35th spectra; the 20th exposure was trimmed from the second STIS G430L light curve. The raw outlier-trimmed HST/STIS broadband transit light curves of WASP-29 and WASP-80 used in our analysis are plotted in the top panels of Figures 2 and 3, respectively.

The HST/WFC3 G141 light curves were extracted from the dark- and bias-corrected *ima.fits files generated by CALWFC3. Following the usual technique for processing spatial scan mode images (e.g., Deming et al., 2013; Kreidberg et al., 2014; Evans et al., 2016), we co-added the background-subtracted difference images from each pair of consecutive nondestructive reads to construct the full frame image. The top and bottom boundaries of each subexposure were determined from the row-added flux profile and included a buffer of 15 rows in both directions beyond the locations at which the flux falls to 20% of its maximum value. We set the background flux level for each difference image as the median pixel value across two rectangular regions on either side of the spectral trace, making sure to avoid the edges of the subarray. The precise number of buffer rows and the locations of the background flux regions did not have any significant effect on the results of the light-curve fits. We corrected for cosmic rays by iteratively applying a 10 pixel wide moving median filter in the $x$ and $y$ directions and replacing $5\sigma$ outliers with the corresponding value in the row-added flux template, scaled to match the total column flux.

As described in detail in Wong et al. (2020a), we followed the method developed by Tsiaras et al. (2016) to compute the full 2D wavelength solution across the subarray. We first measured the relative horizontal offsets $\Delta x$ of each frame by comparing the centroid of the exposure with that of the first frame in the time series. Next, we utilized the calibration coefficients from the configuration file WFC3.IR.G141.V2.5.conf (Kuntschner et al., 2009) and the $x$-position of the star in the accompanying direct image to calculate the location of the spectral trace for a range of stellar $y$-positions and populate the wavelength solution at all pixels in the frame, correcting the $x$-position of the star in each frame by the previously determined offset. This wavelength solution was then used to flatten the full frame images using the coefficients provided in the WFC3.IR.G141.flat.2.fits calibration file (Kuntschner et al., 2011).

To extract the spectroscopic light curves, we defined a 30 nm wavelength grid from 1.0 to 1.7 $\mu$m and used the wavelength solution to determine the precise locations of the wavelength boundaries across each frame. The total flux in each spectroscopic bin was obtained by adding up the pixel values within the corresponding region on the subarray. To accurately account for the partial pixels near the boundaries, we calculated the local 2D cubic interpolation function around each partial pixel and integrated it over the subpixel region to determine the fractional pixel value contained within the bin. The broadband light curve for each visit was obtained by adding up all the flux within the extraction aperture for each full frame image. We discarded the first orbit, as well as the first two data points from each subsequent orbit, prior to fitting. The raw HST/WFC3 transit light curves of WASP-29 and WASP-80 are included in the top panels of Figures 2 and 3, respectively.

The methodology we utilized to extract the IRAC light curves closely mirrors the techniques described in several previous analyses of post-cryogenic Spitzer datasets (e.g., Knutson et al., 2012; Lewis et al., 2013; Todorov et al., 2013; Wong et al., 2015, 2016, 2020a). We determined the background level in each flux-calibrated subarray by masking out the star and the diffraction spikes, as well as the aberrant top row, and then calculating the median from a Gaussian fit to the pixel value distribution. After subtracting the background, we iteratively removed $3\sigma$ outliers across the time series at each pixel position using a moving median filter spanning 64 adjacent images.

The star’s position on the subarray was computed using flux-weighted centroiding (e.g., Knutson et al., 2012). We also calculated the noise pixel parameter $\beta$ (for details, see Lewis et al., 2013), which quantifies the width of the target’s point response function (PRF). We considered a range of circular apertures of radii $r_{0}$ and $r_{1}$ with which to measure the centroid position and $\beta$, respectively, yielding different versions of the photometric series. The final photometric extraction was carried out using both fixed and time-varying circular apertures of radii $r_{\rm phot}$ centered on the centroid position. The radii of the time-varying apertures were defined with respect to the square-root of the noise pixel parameter $\sqrt{\beta}$ via either a multiplicative scaling factor or a constant shift (Wong et al., 2015, 2016, 2020a). Outlier removal was performed on the photometric series, star centroid positions, and the noise pixel parameter array using a 3$\sigma$ moving median filter with a width of 64 data points.

In addition to altering $r_{0}$, $r_{1}$, and $r_{\rm phot}$, we binned the light curves prior to fitting, with bin sizes ranging from 1 to 512 points in powers of two. To remove possible residual flux ramps at the beginning of each time series, we also experimented with trimming the first 15, 30, 45, or 60 minutes of data, while ensuring that at least 30 minutes of pre-ingress photometry remained in the trimmed time series to adequately establish the out-of-eclipse baseline in the fits. For each visit, the optimal choice of the various apertures, bin size, and trimming interval was determined by fitting each version of the photometry to the full transit or secondary eclipse light-curve model (see Section 4) and choosing the dataset that produced the smallest scatter in the residuals, binned in five-minute intervals. The optimal photometric extraction parameters are listed in Table 3 for each visit. The raw Spitzer/IRAC light curves of WASP-29 and WASP-80 are shown in the top panels of Figures 2 and 3.

The ExoTEP pipeline generates the transit and secondary eclipse model light curves using batman (Kreidberg, 2015). In our global analysis of the broadband transit light curves, the transit depth for each visit (parameterized by the planet–star radius ratio $R_{p}/R_{*}$) was allowed to vary independently. We also fit for common orbital ephemeris (mid-transit time $T_{0}$ and orbital period $P$) and transit-shape parameters (impact parameter $b$ and scaled semimajor axis $a/R_{*}$), which are shared across all visits for each system. Both systems are consistent with circular orbits, and we fixed the orbital eccentricity to zero in all of our transit fits.

For the HST/STIS and HST/WFC3 datasets, the limb-darkening coefficients for both broadband and spectroscopic light curves were computed using the LDTk package (Parviainen & Aigrain, 2015), which is incorporated into ExoTEP. From the input stellar parameters of the host star, this program calculates the mean limb-darkening profile for each bandpass from interpolated $50-2600$ nm PHOENIX stellar intensity profiles (Husser et al., 2013) and fits for the limb-darkening coefficients using maximum-likelihood optimization. Following Wong et al. (2020a), we chose the standard four-parameter nonlinear limb-darkening law and fixed the limb-darkening coefficients to the values obtained from LDTk. In the case of the Spitzer transit light curves—for which LDTk is not applicable, due to the wavelength range of the PHOENIX models—we set the limb-darkening coefficients to the tabulated values¹¹1pages.jh.edu/~dsing3/David_Sing/Limb_Darkening.html provided by Sing (2010) for the nearest available combination of stellar parameters.

We experimented with fitting for the quadratic limb-darkening coefficients, but the relatively poor precision of the Spitzer data and the incomplete ingress/egress coverage of the HST light curves did not provide good constraints on the coefficient values, resulting in significantly larger uncertainties on the transit depths. Nevertheless, when comparing the results from the light-curve fits with free and fixed limb darkening coefficients, we found that the astrophysical parameter values agreed with each other to within $1\sigma$ in all cases.

Each raw light curve contains systematics trends that are characteristic of the respective instrument, which we modeled using standard detrending functions.

We employed the prescription described in Sing et al. (2008) as the basis for our HST/STIS systematics modeling. We considered three variants:

The parameters $c$ and $v$ are the normalization constant and the slope of the visit-long linear trend, respectively; $t_{v}$ is defined as the time since the start of the visit. The systematics trend within each spacecraft orbit $F(t_{\rm orb})$ is modeled as a fourth-order polynomial with coefficients $p_{1-4}$; $t_{\rm orb}$ is the time from the start of each orbit. The basic model lin includes just a visit-long linear trend. The linx version has an additional linear trend as a function of the horizontal offset $\Delta x$ of each spectrum relative to the first (see Section 3.1). An additional quadratic term in time is included in the quad variant.

We determined the optimal choice of systematics model for each visit by minimizing the resultant scatter in the residuals as well as the Akaike Information Criterion (Akaike, 1974): ${\rm AIC}\equiv 2\gamma-2\log L$, where $\gamma$ is the number of free parameters, and $\log L$ is the log-probability of the best-fit model. For WASP-29, we found that the two G430L light curves and the G750L light curve prefer the lin, lin, and linx detrending models, respectively. In the case of the three HST/STIS observations of WASP-80, we chose quad, linx, and quad, respectively. The particular choice of systematics models did not have any significant effect on the best-fit astrophysical parameters. However, the optimal combination of models yielded more consistent transmission spectrum shapes between the two G430L visits than other combinations.

For the HST/WFC3 broadband light curves, the systematics model was defined as follows (e.g., Deming et al., 2013; Kreidberg et al., 2014):

Here, the parameters $c$, $v$, $t_{v}$, and $t_{\rm orb}$ are analogous to the corresponding parameters in Equations (1)–(4). The additional piecewise-defined parameter $s^{*}$ was included for the WASP-80 HST/WFC3 visit, where both forward and backward scans were used, in order to account for the common-mode offset between the forward-scan and backward-scan fluxes (see Figure 3): it is set to 1 for the forward scans and $s$ for the backward scans. The second term models the exponential ramp that occurs within each spacecraft orbit: $a$ and $b$ are the rate constant and amplitude. The shape of the ramp in the first fitted orbit (after removing the initial orbit from the original time series; Section 3.2) differs from those in the subsequent orbits. This is accounted for by the piecewise-defined parameter $D^{*}$, which is set to $d$ for points in the first fitted orbit and $0$ everywhere else.

Spitzer IRAC photometry is affected by the well-known intrapixel sensitivity variations. In our analysis, we used Pixel Level Decorrelation (Deming et al., 2015) to detrend these variations. The systematics model is based on a linear combination of the nine pixel-flux arrays $\hat{P}_{k}$ extracted from the $3\times 3$ box centered on the star:

The coefficients $w_{k}$ are the weights on the respective pixel arrays. In addition, we included a visit-long linear trend with slope $v$. The pixel arrays were binned in the same manner as the photometric light curve prior to fitting (see Section 3.3).

ExoTEP simultaneously calculates the best-fit parameter values and uncertainties using the Markov Chain Monte Carlo (MCMC) package emcee (Foreman-Mackey et al., 2013). In addition to the astrophysical and systematics parameters, the pipeline includes the per-point uncertainties $\sigma_{i}$ for each light curve as additional photometric noise parameters. We allowed $\sigma_{i}$ to vary freely to ensure that the resultant reduced $\chi^{2}$ was unity and to self-consistently generate realistic error bars on the fitted parameters.

The global analysis of each planet’s transmission spectrum consisted of three stages. In the first stage, we carried out a joint fit to the HST/WFC3 and Spitzer/IRAC broadband light curves. The Spitzer photometry were prebinned by the intervals that we determined previously from optimizing the photometric extraction of each visit individually (Table 3). In this fit, we allowed all of the orbital ephemeris and transit-shape parameters to vary, along with the three transit depths in the respective bandpasses.

This methodology differs from the one described in Wong et al. (2020a) in that we did not include the HST/STIS broadband light curves in the joint fit. Our choice to exclude those datasets at this stage was motivated by the low signal-to-noise ratio (when compared to the HST/WFC3 datasets) as well as the relatively sparse coverage of the transit ingress and egress, particularly in the case of WASP-80b. When we experimented with including all broadband light curves in the joint fit, we found that the transit-shape parameters and mid-transit time shifted to better match some of the HST/STIS points during ingress and egress, at the expense of a poorer fit to the continuously sampled Spitzer transit light curves.

The second stage of our fitting analysis consisted of fitting the HST/STIS broadband light curves individually, with Gaussian priors placed on $b$, $a/R_{*}$, $T_{0}$, and $P$ using the values determined from the joint fit to the HST/WFC3 and Spitzer/IRAC light curves in the previous stage. For each planet, we analyzed the two STIS G430L visits separately in order to probe for possible stellar variability that may affect the transit depths across the different epochs. After this stage, the global broadband light-curve analysis was complete.

For the final stage, we measured the transmission spectra by fitting the spectroscopic light curves in the STIS G430L, STIS G750L, and WFC3 G141 bandpasses. Here, we allowed only $R_{p}/R_{*}$ to vary, while fixing the transit-shape parameters and orbital ephemeris to the best-fit values from the broadband light-curve fits. Following Wong et al. (2020a), we included the full HST/STIS systematics model when fitting the individual spectroscopic light curves, given that the shape of the systematics trends varies appreciably across the bandpass. Meanwhile, for the WFC3 G141 spectroscopic light-curve fits, we first divided out a common-mode correction from the spectroscopic light curves using the ratio of the raw broadband flux array and the best-fit broadband transit model (e.g., Deming et al., 2013). We then applied a simplified systematics model to account for the remaining instrumental flux variations in each spectroscopic light curve:

This two-parameter model consists of a normalization constant $c$ and a linear slope $v$ with respect to the measured subpixel shift $\Delta x$ of each spectrum in the time series relative to the first one (Section 3.2).

Using the precorrection and simplified systematics modeling for the WFC3 G141 spectroscopic light curves takes advantage of the more uniform shape of the systematics trends in the near-infrared bandpass, reducing the number of free parameters and yielding tighter constraints on the transit depths. No significant improvement to the residual scatter or AIC was found when applying the full HST/WFC3 systematics model to the spectroscopic light curves.

The results of our global broadband light-curve analyses of WASP-29 and WASP-80 are listed in Table 4. The systematics-corrected light curves and best-fit transit models are plotted in Figures 2 and 3. Notably, when comparing the transit depths from the two STIS G430L visits, we find that the values agree with one another to within $1\sigma$ for both planets, indicating the absence of significant stellar photometric variability between the two epochs.

Comparing our best-fit transit shape parameters for WASP-29b with the results from Gibson et al. (2013), we find that our values for $b$ and $a/R_{*}$ are consistent with theirs at better than the $1\sigma$ level, while having significantly smaller uncertainties. Meanwhile, we obtained an updated orbital period of $3.9227090^{+0.0000017}_{-0.0000018}$ d, which is roughly $1.5\sigma$ shorter than the recently published period measurement in Ivshina & Winn (2022): $3.92271159\pm 0.00000038$ days. The Spitzer transits of WASP-29b were previously analyzed by Baxter et al. (2021), who obtained 3.6 and 4.5 $\mu$m transit depths of $9500\pm 100$ and $9300\pm 100$ ppm, respectively, which agree with our measurements at better than the $1\sigma$ level.

For WASP-80b, the literature values of $a/R_{*}$ span a wide range from $12.0647\pm 0.0099$ (Sedaghati et al., 2017) to $12.99\pm 0.03$ (Triaud et al., 2013); our measurement of $12.451^{+0.073}_{-0.071}$ lies intermediate to those two extremes and within $1.6\sigma$ of the values presented in Kirk et al. (2018) and Triaud et al. (2015)—$12.66^{+0.12}_{-0.11}$ and $12.63^{+0.08}_{-0.13}$, respectively. The orbital inclination $i$ we derived from our fits is somewhat lower than the range of values in the literature, lying roughly $2.1\sigma$ below the inclination reported in Kirk et al. (2018). For the orbital period, our result—$3.0678569\pm 0.0000016$ days—agrees well with the most precise period value in the literature: $3.06785500\pm 0.00000036$ days (Parviainen et al., 2018). The Spitzer 3.6 and 4.5 $\mu$m transits of WASP-80b were previously analyzed in Triaud et al. (2015), who reported $R_{p}/R_{*}$ values that are 1.4 and 2.5$\sigma$ smaller than our measurements. This discrepancy is consistent with the significantly smaller inclination value we obtained from our fits.

Given the consistency between the broadband transit depths measured from the two STIS G430L visits of each target, we simply applied a weighted average to the two depths measured in each spectroscopic bin when constructing the final combined transmission spectra. At the blue end of the G430L grating and the red end of the G750L grating, the measured flux from the star drops considerably, resulting in significantly higher uncertainties in the transit depths and larger scatter in the spectrum. In both the table and the figures, we only present values for wavelength bins in which the transit depth $(R_{p}/R_{*})^{2}$ was measured to a precision of $<$250 ppm in the case of WASP-29b and $<$400 ppm in the case of WASP-80b. For the WFC3 G141 bandpass, we do not include the transit depths for the bluest and the two reddest 30 nm wavelength bins (1.10–1.13, 1.64–1.67, and 1.67–1.70 $\mu$m), because the corresponding spectroscopic light curves have significantly poorer signal-to-noise ratios relative to the others.

The measured planet–star radius ratios from the spectroscopic light-curve analyses are listed in Table 5. The full 0.4–5.0 $\mu$m transmission spectra of WASP-29b and WASP-80b are shown in Figures 4 and 5. A compilation of the spectroscopic light curves and corresponding best-fit transit models is provided in the Appendix.

In the STIS G750L bandpass, the transit depths obtained for the narrow wavelength bins centered around the main Na and K absorption regions are broadly consistent with the depths measured in the wider 100-pixel spectroscopic bins spanning those regions, suggesting no significant alkali absorption feature in the transmission spectra of either planet. These narrowband transit depths are not included in the transmission spectrum plots. However, we note that the poor signal-to-noise ratio of these transit-depth measurements makes it difficult to reach definitive conclusions about the presence or absence of alkali absorptions in the transmission spectra. In particular, the relatively deep 577–626 nm transit of WASP-29b may indicate some amount of Na absorption, a possibility that we explore in more detail in Section 6.

The reddest wavelength bin in the STIS G430L grating (538–565 nm) and the bluest wavelength bin in the STIS G750L grating (528–577 nm) overlap. For both planets, the transit depths measured from those two spectroscopic light curves are consistent to well within $1\sigma$. This indicates that there is no evidence for any significant common-mode offsets between the G430L and G750L visits.

Looking at Figure 4, we see that the transmission spectrum of WASP-29b is largely consistent with a flat line, with no detected water absorption feature in the WFC3 G141 bandpass and no evidence for a significant Rayleigh scattering slope in the optical. The shape of the optical transmission spectrum is in good agreement with the previously published ground-based 515–720 nm spectrum from Gibson et al. (2013), which is consistent with a flat line and likewise does not contain any significant alkali absorption features. Tsiaras et al. (2018) published an independent analysis of the WFC3 G141 transmission spectrum as part of their ensemble study. They also found a featureless 1.1–1.7 $\mu$m spectrum, albeit with a slightly more pronounced negative slope across the bandpass than in our spectrum. Moreover, near-infrared transit depths are systematically larger (by roughly 50–100 ppm), owing to the smaller inclination value that they used as a fixed parameter in their fits ($88\overset{\circ}{.}8$; Hellier et al. 2010).

In contrast to WASP-29b, we find a robust 1.4 $\mu$m water absorption feature for WASP-80b. As illustrated in Figure 5, the optical part of the transmission spectrum shows significant scatter, while being broadly consistent with a negative spectral slope that suggests Rayleigh scattering by fine-particle aerosols. Several ground-based transmission spectra have been obtained for WASP-80b spanning the visible wavelength range. Focusing on those with spectral resolutions greater than $\lambda/\Delta\lambda=10$, we find a range of spectral slopes. The 500–900 nm spectrum published by Parviainen et al. (2018) is consistent with a flat line, while the 500–900 nm spectrum measured by Kirk et al. (2018) shows a negative spectral slope consistent with our spectrum. The 750–1000 nm spectrum published by Sedaghati et al. (2017) has a very steep negative slope, possibly caused by uncorrected residual systematics in their VLT/FORS light curves. The WFC3 G141 transmission spectrum published by Tsiaras et al. (2018) displays a prominent water absorption feature at 1.4 $\mu$m, consistent in shape with our results. Meanwhile, our measured near-infrared transit depths are roughly 200 ppm larger on average, due to the smaller inclination we obtained from our broadband light-curve fit.

The Spitzer secondary eclipse light curves were fit using the same PLD systematics model as the transits. Due to the low signal-to-noise ratio, we fixed the transit-shape parameters and orbital ephemeris to the best-fit values obtained from our broadband transit light-curve analysis (Table 4). When fitting the individual visits to optimize the photometric extraction parameters and binning/trimming, we set the mid-eclipse time to an orbital phase of 0.5. The measured depths for all secondary eclipse light curves are listed in Table 6. The binned and systematics-corrected WASP-29b and WASP-80b eclipse light curves are plotted in Figures 6 and 7, respectively. Most of the eclipse-depth measurements are marginal, particularly at 3.6 $\mu$m. When comparing the individual eclipse depths in the same bandpass, we find that they are mutually consistent at better than the $2\sigma$ level in all cases.

For each system, we carried out a joint fit to all secondary eclipse observations. In order to probe for nonzero orbital eccentricity, we allowed the mid-eclipse orbital phase to vary freely. All eclipse light curves within a bandpass were assigned a common eclipse depth. The results from our joint fits are presented in Table 6. We have corrected the mid-eclipse orbital phases by the light travel time between the superior and inferior conjunctions. For WASP-29b, the best-fit 3.6 and 4.5 $\mu$m eclipse depths are $0.015\%\pm 0.008\%$ and $0.033\%\pm 0.006\%$, respectively, and the mid-eclipse orbital phase is $0.5034^{+0.0012}_{-0.0016}$, roughly $2.1\sigma$ after mid-orbit. For WASP-80b, the eclipse depths are $0.028\%\pm 0.005\%$ and $0.093\%^{+0.005\%}_{-0.006\%}$, respectively; the best-fit eclipse phase is consistent with 0.5 at better than the $1\sigma$ level, indicating that the orbit is consistent with circular.

Triaud et al. (2015) presented an independent analysis of the WASP-80b eclipses obtained as part of their Spitzer program, i.e., Eclipses 1 and 2 in each bandpass, as listed in Table 3. From a global analysis of the Spitzer transit and secondary eclipse light curves, along with several other ground-based transit observations, they obtained best-fit eclipse depths of $0.0455\%\pm 0.0100\%$ and $0.0944\%^{+0.0064\%}_{-0.0065\%}$ at 3.6 and 4.5 $\mu$m, respectively. Their 4.5 $\mu$m depth agrees with our value to well within $1\sigma$, while their 3.6 $\mu$m measurement is larger than ours by roughly $1.6\sigma$. We reran the global eclipse light-curve analysis without the third and fourth 3.6 $\mu$m secondary eclipse light curves and obtained a depth of $0.35\%\pm 0.09\%$, which is consistent with the Triaud et al. (2015) value at better than the $1\sigma$ level.

From the global eclipse depths listed in Table 6, we derived the corresponding dayside brightness temperatures, following the methodology described in Wong et al. (2020a). We computed the stellar flux at 3.6 and 4.5 $\mu$m, integrated across the Spitzer bandpasses, for a grid of PHOENIX stellar models and fit for polynomial functions with respect to stellar brightness temperature, metallicity, and surface gravity. These empirical relationships were used to self-consistently propagate the uncertainties in the measured stellar parameters (Table 1) to our estimates of the planet’s dayside brightness temperature. The eclipse depth depends on the star–planet radius ratio as well, and we utilized the most precise value from our broadband light-curve analysis—the transit depth in the WFC3 G141 bandpass. Because the primary contribution to the uncertainty in the brightness temperature is the eclipse-depth uncertainty, choosing any of the other listed $R_{p}/R_{*}$ values in Table 4 did not significantly affect our temperature estimates.

The 3.6 and 4.5 $\mu$m dayside brightness temperatures for WASP-29b are $830^{+150}_{-440}$ and $940^{+48}_{-52}$ K, respectively. Fitting a single blackbody to both eclipse depths yielded a dayside brightness temperature of $937^{+45}_{-48}$ K. For WASP-80b, the dayside brightness temperatures in the two Spitzer bandpasses are $791^{+28}_{-30}$ and $871^{+15}_{-16}$ K, respectively. The two eclipse depths are consistent at better than $2\sigma$ with a single blackbody at $851^{+13}_{-14}$ K. For both WASP-29b and WASP-80b, our measured dayside brightness temperatures are in excellent agreement with the equilibrium temperatures listed in Table 1 ($970^{+32}_{-31}$ and $825\pm 19$ K, respectively). The equilibrium temperatures were calculated assuming zero Bond albedo and no heat recirculation to the nightside. Therefore, we conclude that both planets have weak day–night heat transport.

General circulation models of temperate planets with similar equilibrium temperatures to WASP-29b and WASP-80b have generally predicted more efficient global heat transport and correspondingly reduced day–night temperature contrasts when compared to their hotter counterparts (e.g., Showman et al., 2015). Therefore, at first glance, our inference of relatively weak day–night heat transport for these two planets stands in contradiction to the prevailing trend. However, metallicity has been shown to also affect the global energy budget of planetary atmospheres, with some studies showing that higher metallicities inhibit day–night heat recirculation (e.g., Zhang & Showman, 2017). The enhanced metallicities that we retrieve from the transmission spectra (Section 6) may indicate that the metal enrichment is weakening longitudinal heat transport.

More broadly, there is a marked dearth of relatively cool exoplanets that have been examined in detail with sophisticated atmospheric modeling (e.g., Zhang, 2020). This has primarily been caused by the lack of high-quality emission data for these targets. The temperature distribution of warm gas giants such as WASP-29b and WASP-80b has hitherto only been constrained with broadband Spitzer observations. The James Webb Space Telescope (JWST) will enable intensive spectroscopic study of these cooler objects and definitively probe fundamental trends in atmospheric heat transport in this temperature regime.

The mostly flat and featureless transmission spectrum of WASP-29b indicates the presence of clouds and/or very high atmospheric metallicity that attenuate the absorption features that would otherwise be detectable, particularly in the WFC3 G141 bandpass. Our retrieval runs invariably returned cloudtop pressures in the range 0.1–1 mbar. A few complications arose during our retrieval analysis. The first is the low Spitzer 4.5 $\mu$m transit depth, which pushed the limb temperature posterior against the lower edge of the temperature range allowed by PLATON (300 K). Given that WASP-29b has an equilibrium temperature of roughly 970 K (Table 1), a limb temperature at or below 300 K is physically implausible. The second complication is the low average transit depth across the STIS bandpasses. In particular, the average level of the visible transmission spectrum is lower than the average depth in the WFC3 G141 bandpass. As demonstrated by our comparisons of retrieval runs that included or excluded the STIS transmission spectrum (see below), the relatively low optical transit depths push the limb temperature down to very low values, while simultaneously forcing a mostly cloud-free atmosphere and extremely high metallicity.

To address the issue of the STIS transit depths, we included an offset (with a prior standard deviation of 100 ppm) and ran retrievals on the transmission spectrum without the Spitzer points. The resultant limb temperatures were more reasonable. However, when we reintroduced the Spitzer data points into the retrieval and allowed for offsets between the Spitzer depths and the WFC3 points, we still found a strong preference for extremely low temperatures. Therefore, we decided to exclude both Spitzer transit depths from the retrieval analysis presented here. We also removed the high outlier point at 1.52–1.55 $\mu$m.

Table 8 compiles the Bayesian evidence, maximum likelihood value, reduced $\chi^{2}$, and degrees of freedom of some of our retrievals. Due to the large uncertainties in the optical portion of the transmission spectrum, the reduced $\chi^{2}$ values are much lower than unity, indicating that the level of photometric uncertainty is preventing the retrievals from distinguishing among different atmospheric models. The extent to which a given model is disfavored relative to the preferred model (i.e., the run with the highest Bayesian evidence) is quantified by the Bayes factor $B$. The fiducial model with a STIS offset performed the best. Removing the STIS offset yielded a lower Bayesian evidence and lower likelihood and is disfavored with a Bayes factor of 3.03, which approximately corresponds to a $2.1\sigma$ significance level. We experimented with adding fine-particle hazes, but the poorly constrained shape of the transmission spectrum at optical wavelengths means that a wide range of scattering factors and slopes is allowed, with no significant improvement to the likelihood. Meanwhile, a cloud-free atmosphere is disfavored at the $1.9\sigma$ level.

When allowing for both free sodium abundance and hazes, we retrieved a broad range of $\log\,{\rm Na}$ values as high as $-$4. For comparison, the equilibrium sodium abundance at pressures below 1 mbar in the best-fit fiducial model is $\log\,{\rm Na}<-8$. While the fit quality is somewhat improved, the increase in model complexity is penalized by the Bayesian evidence. Notably, even at the upper end of the $\log\,{\rm Na}$ posterior, the model spectrum does not show a very prominent sodium absorption feature. This can be attributed to the necessity for high metallicity and/or significant cloud cover in order to match the null water vapor feature at 1.4 $\mu$m, which works to attenuate the sodium absorption as well. One would have to simultaneously invoke relatively low atmospheric metallicity, depleted water abundance, and significantly enhanced sodium abundance to yield both a strong sodium absorption feature and a flat near-infrared transmission spectrum. Given the poor precision of the measured optical transit depths and the general observation that the overall HST STIS transmission spectrum is consistent with a flat line to within the error bars, we did not pursue the issue of sodium abundance further. More intensive observations (e.g., using large ground-based facilities) are needed in order to resolve the transmission spectrum in the vicinity of the sodium absorption feature.

In Figure 4, we show the results of our preferred retrieval run. The $2\sigma$ confidence region is consistent with completely flat transmission spectra. The retrieved parameter values are listed in Table 9, while the corner plot of relevant parameters can be found in the Appendix (Figure 12). The retrieval returned a positive offset of $78\pm 37$ ppm for the STIS transit depths. The planetary radius estimate of $0.747\pm 0.036$ $R_{\rm Jup}$ is consistent with the literature value ($0.776\pm 0.043$ $R_{\rm Jup}$; Gibson et al. 2013). Due to the low precision of the transit depths at visible wavelengths and the lack of robust spectroscopic features, the limb temperature is poorly constrained, with the $1\sigma$ region spanning 390–790 K. Likewise owing to the lack of discernible absorption features amid significant cloud opacity, the full range of C/O ratios is allowed, including near-solar values and ${\rm C/O}>1$.

Overall, our retrievals indicate cloudy limbs ($P_{\rm cloud}$ from $<$0.01 mbar to $\sim$1 mbar) and unconstrained atmospheric metallicity ($\log(Z/Z_{\odot})=1.1\pm 1.4$). Inspecting the corner plot, we find a slight hint of the classic degeneracy between cloudtop pressure and atmospheric metallicity: the $3\sigma$ confidence region includes scenarios with cloud-free limbs and very high metallicity ($>$$500\times$ solar). Interior modeling of gas giants has revealed a broad range of possible atmospheric metallicities, contingent upon factors such as the level of internal mixing, with the upper limit set at the planet’s bulk metallicity value. The bulk metallicity of WASP-29b has a $2\sigma$ upper limit of $190\times$ solar (Thorngren & Fortney, 2019), therefore we consider the cloud-free scenario to be highly unlikely.

The transmission spectrum of WASP-80b stands in stark contrast to the WASP-29b spectrum and contains both a prominent scattering slope in the optical and a pronounced absorption feature around 1.4 $\mu$m. We found that the fiducial model with enhanced Rayleigh scattering was the best-performing model out of the various retrieval runs. Looking at Table 8, we can see that the model without fine-particle aerosol opacity is strongly disfavored ($B=36300$, $5.0\sigma$), while allowing the scattering slope to vary did not improve the likelihood or Bayesian evidence. Likewise, employing the more complicated Mie scattering model yielded a slightly poorer fit to the visible portion of the transmission spectrum, with Titan tholins performing slightly better than soots. Both Mie scattering retrieval runs indicated particle sizes less than 0.1 $\mu$m, in agreement with the observed Rayleigh scattering slope throughout the STIS bandpass. Regardless of the choice of scattering parameterization, other important atmospheric parameters such as atmospheric metallicity, C/O ratio, and cloudtop pressure did not change appreciably, indicating that those properties are not primarily determined by the visible portion of the transmission spectrum. Similarly, the inclusion of a common-mode offset for the STIS transit depths is marginally disfavored by the Bayesian evidence, while yielding an insignificant improvement to the maximum log-likelihood and reduced chi-squared.

The results of the preferred retrieval model are shown in Table 9 and Figure 5; the corner plot is displayed in Figure 13. The quality of the fit is excellent, and the introduction of offsets between the various instrumental datasets did not yield improved fits. The water absorption feature in the WFC3 bandpass is well-modeled, as is the relatively deep Spitzer 4.5 $\mu$m transit depth, which suggests CO and/or CO₂ absorption. To quantify the significance of the molecular detections, we carried out a series of retrievals, each time removing the opacity contribution from one major atmospheric species (selected from among H₂O, CO, CO₂, and CH₄). This process did not affect the equilibrium abundances calculated by PLATON, but instead simply set the corresponding absorption cross section to zero. The Bayesian evidences of these retrieval runs were then compared to the fiducial model with hazes. As shown in Table 8, the removal of H₂O and CO₂ leads to significantly poorer fits, with Bayes factors of 1290 and 68.7, respectively; meanwhile, we did not find a statistically significant spectroscopic contribution from CO or CH₄ in the measured transmission spectrum (but see discussion in Section 6.3).

The limb temperature of WASP-80b is well-constrained by the transmission spectrum: $930^{+120}_{-110}$ K. This value is consistent with the dayside temperature that we derived from the Spitzer secondary eclipse depths—$851^{+13}_{-14}$ K (Section 5.3)—indicating relatively uniform heat distribution across the dayside hemisphere and extending to the limbs. Due to the presence of the H₂O absorption feature, the C/O ratio is restricted to values less than 0.9 (e.g., Benneke, 2015), but otherwise unconstrained. The atmospheric metallicity is enhanced, with a mean value of roughly $100\times$ solar and a relatively narrow posterior that excludes near-solar metallicities at greater than $3\sigma$ significance. The corner plot shows a significant cloud–metallicity degeneracy. While the retrieval mostly prefers cloudy atmospheres with atmospheric metallicities of 10–$100\times$ solar, cloud-free limbs are still allowed for metallicities of $\sim$100$\times$ solar and higher. This latter observation was also supported by our cloud-free retrieval, which provided a fit quality similar to that of the fiducial model that includes both clouds and hazes and strongly constrained the atmospheric metallicity to values above roughly $80\times$ solar.

The retrieved $1\sigma$ confidence region of WASP-80b’s atmospheric metallicity (30–170$\times$ solar) is notably higher than that of any other planet larger than 0.5 $M_{\rm Jup}$. Thorngren & Fortney (2019) reported a bulk metallicity of $0.17\pm 0.04$, corresponding to a maximum atmospheric metallicity of $30\pm 7\times$ solar if the planet is fully mixed, with a $2\sigma$ upper limit of $43\times$ solar. This theoretical range is consistent with the lower end of our retrieved metallicities to within $1\sigma$ and suggests that the atmosphere of WASP-80b is relatively well-mixed. This physical argument also disfavors scenarios of cloud-free limbs with atmospheric metallicities above $100\times$ solar.

An important caveat to our retrieved atmospheric metallicity is the critical role played by the large Spitzer 4.5 $\mu$m transit depth: the high CO/CO₂ abundance needed to explain this single point strongly drives the retrieved metallicity to enhanced values. To quantify the effect of the Spitzer portion of the transmission spectrum on the retrieved metallicity, we carried out analogous retrievals on the WFC3 and WFC+STIS spectra. From the WFC3-only retrieval, the atmospheric metallicity is much less constrained, displaying a strong degeneracy with the cloudtop pressure; the $1\sigma$ confidence region of atmospheric metallicity is 3–$150\times$ solar, with a median of $20\times$ solar. Meanwhile, the WFC3+STIS spectrum yields a somewhat narrower range of atmospheric metallicities (10–$130\times$ solar), with a higher median of $40\times$ solar. It is apparent that, without the Spitzer points, the transmission spectrum of WASP-80b becomes broadly consistent with a wide range of atmospheric metallicities, from near-solar to several hundreds of times solar, while being consistent with the tighter constraints from the full dataset retrieval. The results from our Spitzer-free retrievals, along with the theoretical arguments detailed in the previous paragraph and the problematic nature of the 4.5 $\mu$m transit depth in the case of WASP-29b, suggest that the 4.5 $\mu$m transit depth of WASP-80b may be overestimated, biasing our retrieved metallicities to unphysically high values. Future JWST observations (see Section 6.3) are needed to obtain an accurate and detailed picture of the transmission spectrum in the 3–5 $\mu$m region and definitively constrain the atmospheric metallicity, cloudtop pressure, and CO/CO₂ abundances.

Our photometric monitoring of WASP-80 (Section 2.3) revealed discernible stellar variability that was synchronous with the star’s rotational period, raising the possibility of starspots on the stellar surface. The appreciable activity of WASP-80 is also corroborated by the X-ray emission observations and the detection of a flare (Section 2.4). As unocculted starspots can alter the shape of the transmission spectrum in ways that mimic a Rayleigh scattering slope, we experimented with fitting for haze-free atmospheres while including starspots. The starspot temperature was fixed to the predicted value from the empirical linear scaling law in Rackham et al. (2019): $T_{\rm spot}=3350$ K. We retrieved $f_{\rm spot}=0.033\pm 0.005$, which is consistent with the range of modeled spot-covering fractions for K6V stars in the Kepler sample ($0.014^{+0.027}_{-0.007}$; Rackham et al. 2019). As shown by the log-likelihood values in Table 8, unocculted spots provide a comparably good fit to the steep optical slope in the transmission spectrum in the absence of hazes. We note that, due to the distinct parameterization and correspondingly different prior volume of the hazeless unocculted spot model, the Bayesian evidence cannot be straightforwardly compared with that of the fiducial model with hazes. An important observation is that the spot-covering fraction must be precisely adjusted in order to simultaneously create the apparent Rayleigh scattering slope and produce transmission features that mimic those with the expected atmospheric temperature. Nevertheless, given the self-consistency between the measured photometric variability, retrieved spot-covering fraction, and the corresponding impact on the optical transmission spectrum, it is possible that contamination from unocculted starspots may be at least partially responsible for the observed spectral slope.

Looking at the results of our atmospheric retrievals, we find that the featureless transmission spectrum of WASP-29b hinders detailed quantitative comparisons between the two planets. For example, while the atmospheric metallicity of WASP-80b is largely constrained to highly enriched values, the WASP-29b transmission spectrum is consistent with the full range of allowed metallicity values, from subsolar to 1000$\times$ solar. Nevertheless, several salient differences are apparent from the overall shape of the transmission spectra.

First, the steep optical slope apparent in the WASP-80b transmission spectrum is entirely absent in the WASP-29b dataset. We reiterate that there is some ambiguity in the interpretation of the optical slope for WASP-80b, with unocculted starspots providing a plausible alternate or contributing explanation for the feature. If we assume that unocculted starspots did not significantly affect the measured transmission spectrum, then the optical slope is indicative of Rayleigh scattering from aerosols with particle sizes smaller than $\sim$0.1 $\mu$m. Such haze particles are readily produced by photochemical processes in the upper atmosphere and are expected to be the dominant source of aerosol opacity at temperatures below $\sim$950 K (e.g., Lavvas & Koskinen, 2017; Hörst et al., 2018; Adams et al., 2019; Gao et al., 2020; Ohno & Kawashima, 2020). Meanwhile, the flat spectrum of WASP-29b signals the absence of such fine-particle opacity.

Flat spectra with obscured molecular absorption features are usually attributed to condensate clouds consisting of $>$1 $\mu$m aerosol particles that extend to low pressure levels (e.g., Fortney, 2005; Helling et al., 2008; Morley et al., 2013; Charnay et al., 2015; Barstow et al., 2017; Gao & Benneke, 2018; Powell et al., 2019). Similar spectra have been measured for planets spanning a wide range of equilibrium temperatures (e.g., Deming et al., 2013; Knutson et al., 2014; Kreidberg et al., 2014; Sing et al., 2016; Wakeford et al., 2017; Tsiaras et al., 2018; Chachan et al., 2020; Wakeford et al., 2020; Alam et al., 2022). Silicate clouds (particularly Mg₂SiO₄) are expected to be the primary condensate species contributing to enhanced opacity at the equilibrium temperature of WASP-29b, and indeed across most of the temperature range from $\sim$900 to $\sim$2000 K (e.g., Gao et al., 2020). However, recent modeling has shown that photochemical hazes composed of hydrocarbon aggregates can grow to larger sizes than the spherical haze particles typically considered in earlier studies, resulting in a weaker wavelength dependence of the optical depth that can mimic the signature of condensate cloud decks (e.g., Adams et al., 2019; Ohno et al., 2020). Therefore, we cannot exclude the possibility that the gray opacity observed on WASP-29b is driven, at least partially, by photochemistry. In fact, both WASP-29b and WASP-80b lie in a temperature range where photochemical haze production is expected to contribute significantly to the atmospheric opacity on planet-wide scales; laboratory experiments and numerical models have demonstrated high levels of haze production in temperate and warm exoplanet atmospheres with $T_{\rm eq}\leq 1000$ K (e.g., Fortney et al., 2013; Morley et al., 2015; He et al., 2018; Hörst et al., 2018; Kawashima & Ikoma, 2019; Yu et al., 2021).

Another related point of discrepancy between WASP-29b and WASP-80b is the amplitude of the water-vapor absorption feature at $\sim$1.4 $\mu$m: while WASP-80b’s transmission spectrum displays a very prominent absorption feature, corresponding to a $4.2\sigma$ detection of H₂O, the WASP-29b spectrum shows no discernible feature. The so-called water amplitude $A_{H}$ quantifies the flatness of the transmission spectrum in the near-infrared and has become a commonly used metric for assessing broad trends in the impact of aerosol opacity on the growing ensemble of measured transmission spectra (e.g., Sing et al., 2016; Stevenson, 2016; Crossfield & Kreidberg, 2017; Fu et al., 2017). The most recent statistical analysis by Dymont et al. (2021) carried out a multivariate search for trends among the body of available WFC3 transmission spectra for planets with $T_{\rm eq}<1000$ K and uncovered significant covariances between $A_{H}$ and atmospheric scale height, planet gravity, and planet density.

Most of the published spectra used by Dymont et al. (2021) were taken from the large population study by Tsiaras et al. (2018), which included WASP-29b and WASP-80b. We calculated new $A_{H}$ values for these planets based on our measured transmission spectra, following the techniques described in Dymont et al. (2021). After taking the weighted average of the transit depths measured in the two wavelength bins straddling 1.25 and 1.4 $\mu$m (Table 5), we obtained $A_{H}$ values of $-0.11\pm 0.50$ and $1.00\pm 0.52$ for WASP-29b and WASP-80b, respectively. Both of these results are consistent with the values in Dymont et al. (2021): $-0.17\pm 0.47$ and $0.61\pm 0.29$.