EDP Sciences logo
Subscriber Authentication Point
Search
Menu
Home All issues Volume 707 (March 2026) A&A, 707 (2026) A303 Full HTML
Open Access
Issue
A&A
Volume 707, March 2026
Article Number A303
Number of page(s) 12
Section Planets, planetary systems, and small bodies
DOI https://doi.org/10.1051/0004-6361/202558047
Published online 24 March 2026

© The Authors 2026

Licence Creative CommonsOpen Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

This article is published in open access under the Subscribe to Open model. This email address is being protected from spambots. You need JavaScript enabled to view it. to support open access publication.

1 Introduction

Characterizations of exoplanetary atmospheres are vital for addressing some key questions in the field of exoplanets, such as the formation and evolution of planets, the properties and evolution of planet atmospheres, and the detection of possible biosignatures. Over the past two decades, atmospheric studies have predominantly focused on transiting close-in gas giants, particularly the hot Jupiters (HJs) and ultra-hot Jupiters (UHJs), due to their relatively strong atmospheric signals. Joint analysis based on data across broad wavelengths and/or multiple resolutions enables constraints on atmospheric temperature-pressure profiles, chemical composition, dynamics, and escape processes.

Ultra-hot Jupiters, which have high equilibrium temperatures Teq>2200K (Parmentier et al. 2018; Stangret et al. 2022; Tan et al. 2024) and highly inflated atmospheres, represent a distinct population of planets and are the most favorable targets for time-resolved spectroscopic observations such as transmission and emission spectroscopy (Snellen 2025). Recent abundance retrievals for some HJs and UHJs, for example WASP-121b (Smith et al. 2024; Evans-Soma et al. 2025) and WASP-18 b (Sheppard et al. 2017), seem to be inconsistent with predictions for this population, particularly the expected inverse correlation between metallicity and the carbon-to-oxygen (C/O) ratio (Espinoza et al. 2017; Cridland et al. 2019). Although mechanisms such as pebble drift (Booth et al. 2017) could potentially explain these abundances, the observed discrepancies invoke the demand for more observation input to test and refine current theories.

Moreover, several UHJs show significant spin-orbit misalignments, including MASCARA-5 b (Stangret et al. 2021), KELT-18b (Rubenzahl et al. 2024), TOI-1518b (Cabot et al. 2021), KELT-9b (Gaudi et al. 2017), WASP-121 b (Delrez et al. 2016), and WASP-33 b (Watanabe et al. 2022), indicating that UHJs experience dynamically extreme evolutionary histories, such as the eccentric Lidov-Kozai effect (Naoz et al. 2011). Atmospheric studies of these systems may therefore help constrain when and where planetary migration occurred.

UHJs experience stellar irradiation 10-100 times stronger than classic HJs, and 2-6 orders of magnitude stronger than the warm or cool planets. The intense UV flux drives high temperatures that produce thermal inversions in their upper atmospheres (Baxter et al. 2020). Within and above the temperature inversion layer on the dayside, most molecules such as H2O, TiO, and VO undergo thermal dissociation, while CO can still remain relatively abundant (Madhusudhan 2012; Moses et al. 2012; Drummond et al. 2019). Therefore, atmospheres become dominated by atomic and ionic species (e.g., Fe I/II, Mg I, Ca II, Na I, Ti I, and V I).

Observations from ground-based high-resolution spectroscopy (HRS) and space-based low to medium-resolution spectroscopy support the transition from molecular to atomic and ionic species in UHJs, reflecting extensive thermal dissociation at high temperatures, although nightside cold-trapping and rainout at the terminator that have been observed (Gandhi et al. 2023; Hoeijmakers et al. 2024) may also play a role. Some UHJs exhibit extended envelopes or tails as revealed from transit or near-transit observations, suggesting extreme mass loss (Yan & Henning 2018; Cabot et al. 2020; Yan et al. 2021). These observational atmospheric properties are highly valuable for understanding the composition of the UHJs atmosphere and refining 3D atmospheric modeling with general circulation models (GCMs).

The transiting UHJ, WASP-33 b, discovered by Christian et al. (2006) in the SuperWASP project (Pollacco et al. 2006), has an orbital period of ~1.22 days. The follow-up observation measurements yielded a mass of 2.1 MJ and a radius of 1.6 RJ (Chakrabarty & Sengupta 2019). Notably, its host star WASP-33 (HD 15082) is a bright, rapidly rotating A5 δ Scuti star (Grenier et al. 1999; Cutri et al. 2003), making it a peculiar exoplanet. Subsequent Rossiter-McLaughlin effect (RME) observations (Johnson et al. 2015; Stephan et al. 2022; Watanabe et al. 2022) revealed a near-polar orbit, showing strong resemblance to two other UHJs around A-type stars, i.e., KELT-9 b (Gaudi et al. 2017) and WASP-189b (Anderson et al. 2018; Lendl et al. 2020). This configuration suggests an intense migration history, potentially involving mechanisms such as eccentric Lidov-Kozai effects (Naoz et al. 2011). Interestingly, Mugrauer (2019) pointed out that WASP-33 is at least a binary or even a hierarchical triple star system, including a close-in M-dwarf companion (WASP-33B; ρ ~ 1.9″, ∆ Ks = 6.11 ± 0.02; Moya et al. 2011; Wöllert & Brandner 2015; Ngo et al. 2016) and a wide-separation G-dwarf companion (WASP-33C, ρ ~ 49.0″; Mugrauer 2019), invoking the possibility of dynamic perturbation on the planet in the past and/or still ongoing.

Multi-resolution studies in a broad wavelength range have been conducted extensively for the characterization of WASP-33 b. Early photometric eclipse measurements at 0.91 μm, the Ks band, and the Spitzer/IRAC broadband yielded notably high dayside brightness temperatures (≥ 3300 K), suggesting a poten-tial thermal inversion supported or a high C/O ratio (Smith et al. 2011; Deming et al. 2012; De Mooij et al. 2013). The first eclipse spectrum, obtained by Haynes et al. (2015) using HST/WFC3, showed excess flux at short wavelengths that was attributed to TiO emission. However, a subsequent reanalysis of these data by Changeat et al. (2022) and HRS has questioned the presence of TiO (Nugroho et al. 2017; Cont et al. 2021; Herman et al. 2020; Yang et al. 2024a). Recent optical and Near-infrared (NIR) HRS detected atomic species, including Fe I (Nugroho et al. 2020), Si I (Cont et al. 2022), VI and Ti I (Cont et al. 2022), molecules, including CO (van Sluijs et al. 2023; Yan et al. 2022), OH (Nugroho et al. 2021), and a weak H2O signal (Nugroho et al. 2021), which reveal thermal inversion and strong dissociation on the dayside. Transmission spectroscopy detected absorptions from Ca II (Yang et al. 2024a) and AlO (von Essen et al. 2019), and probed an extended envelope with an escape rate of ~1012 g s−1 by measuring additional absorption of Balmer lines. Notably, Yang et al. (2024b) identified a possible contribution of nightside H2O emission to transmission spectra, suggesting the presence of nightside water, and that WASP-33 b may exhibit nonisothermal temperature structures on its nightside. Phase-curve observations reveal inefficient day-night heat redistribution (Zhang et al. 2018; Herman et al. 2022; Dang et al. 2024) and CO absorption on the nightside, making WASP-33 b the second exoplanet with confirmed nightside molecules (Mraz et al. 2024). However, most constraints remain limited to species detection, with few atmospheric retrievals conducted (Haynes et al. 2015; Changeat et al. 2022; Finnerty et al. 2023). These retrievals consistently indicate a dayside thermal inversion, as predicted by Mollière et al. (2015) for atmospheres with C/O ≈1 and Teff > 1500 K, yet they differ in molecular abundances and profile details, highlighting atmospheric complexity that requires further observational constraints.

In this work, we report two ground-based high-precision eclipse observations of ultra-hot Jupiter WASP-33 b obtained with the Canada-France-Hawaii Telescope (CFHT) in 2015 using the CH4on filter and our customized CO filter. Section 2 details the observational setup, followed by a brief description of data reduction and analysis in Sect. 3. Section 3.5 describes our atmospheric retrieval methodology and results, followed by the summary in Sect. 5.

2 Observations

We observed two secondary eclipses of WASP-33 b (Program 15AS007; PI: Wei Wang) using the Wide-field InfraRed Camera (WIRCam)1 on the CFHT. WIRCam features a 21 arcmin × 21 arcmin field of view that enables differential photometry with multiple reference stars. Observations occurred on October 25, 2015 in the CO filter (5.2 hours duration) and on 24 December, 2015 in the CH4on filter (6 hours duration). As Fig. 1 shows, both of them fully covered the secondary eclipse along with ~2.5 hours of out-of-eclipse baseline. We acquired 696 exposures (12s each) in CO and 996 exposures (3.5 s each) in CH4on. In order to achieve differential photometric precision to 10−4, the staring mode(Devost et al. 2010) was used in both observations, following previous successful attempts (e.g., Croll et al. 2010; Wang et al. 2013; Chen et al. 2014; Shi et al. 2023). The staring mode aims to maintain the telescope’s pointing stability to a few pixels, thereby minimizing the impact of intrapixel and interpixel variations of the detectors on photometric precision. However, our data show that the target star exhibits large movement with maximum shifts of ~15 pixels in the CO filter. This drift likely results from the guiding system misidentifying a cluster of hot pixels on Detector #77 as a stellar point during observations. The misidentification also leads to World Coordinate System (WCS) interpretation failures for the affected frames, rendering them unusable. In contrast, target movement remains within 5 pixels in the CH4on filter. Figure 2 shows the filter information of the CFHT/WIRCam’s CO and CH4on filters, whose λmean are 2.32 μm and 1.69 μm. The CO filter is customized to be sensitive to CO features at about 2.3 μm, while the CH4on filter is sensitive to both CH4on and H−. We use these two filters to assess the presence and abundances of CO, CH4on, and H−, and thus to provide constraints on the C/O ratio and metallicity.

Thumbnail: Fig. 1 Refer to the following caption and surrounding text. Fig. 1

Orbit phases of our two observations. Specifically, Night 1 covered phase 0.410-0.589 (October 25, 2015), Night 2 covered phase 0.403-0.573 (December 24, 2015).
Thumbnail: Fig. 2 Refer to the following caption and surrounding text. Fig. 2

Transmission curves of the CO and CH4on filters used in this work.

3 Data analysis

3.1 Raw data reduction and correction

Observations obtained with CFHT/WIRCam utilized the sample-up-the-ramp (nondestructive readout) mode (Finger et al. 2008), with each integration sequence consisting of 12 reads. Our analysis reveals that the initial read of each sequence on Detector #52 consistently exhibited elevated counts, producing photometric values 1.002-1.004 times higher than the median of other reads, which cannot be reliably corrected. Hence, we excluded the first read of each sequence and only used the remaining 11 reads in this work.

Raw data were reduced by the ‘I‘iwi pipeline (version 2.1.200)2, including flagging saturated pixels, nonlinearity correction, reference pixels subtraction, dark subtraction, flat fielding, bad pixels masking, and guide window masking. Technically, the ‘I‘iwi pipeline should assign bad pixels and saturated pixels’ count a value of Not-a-Number (NaN). However, we find that several pixels within the stellar profile are misclassified as saturated pixels by ‘I‘iwi, resulting in their assignment of NaN values. This misclassification further exacerbates the challenges associated with our data reduction and analysis. Figure 3 shows a full-frame CO filter’s reduced image illustrating our observation field of view (FOV), in which the target star is marked as the red square, and the final chosen reference stars are marked as the blue circles.

As shown in Fig. 3, the WIRCam detectors contain a large number of bad pixels. While reference stars were preferentially selected to avoid regions with clustered bad pixels, there are still some isolated defective pixels occasionally remaining within photometric apertures of some reference stars in some exposures, leading to strong systematic in the yielded lightcurves and thus preventing a high-precision differential photometry. To mitigate this effect, we used the Background2D module of the Photutils package (Bradley et al. 2024) for spatial interpolation of NaN values, following the methodology outlined by Ding et al. (2025). This technique can estimate the data values for the affected pixels, demonstrating reliable performance in handling such detector imperfections.

Thumbnail: Fig. 3 Refer to the following caption and surrounding text. Fig. 3

CO filter reduced full-frame WIRCam image of WSAP-33 taken on October 25, 2015. The target star is marked as the white square and the final chosen reference stars are marked as the grey circles.

3.2 Photometry and light curve determination

To achieve a photometric precision down to 10−4 that is required for the characterization of exoplanet atmospheres, our observations were carried out in defocused mode. With this configuration, the observation efficiency is maximized to gain more photons in the eclipse, and the uncertainties in flat-fielding, intrapixel, and interpixel sensitivity variations are minimized. The disadvantage of such a setup is that the point spread function (PSF) of a star becomes donut-shaped rather than Gaussianshaped. In our case, the resulting PSF has a radius of ~10-15 pixels. We have developed a new pipeline, written in Python, based on and improved on the IDL pipeline used in Shi et al. (2023). Our pipeline employs a parallel processing approach, reducing the photometry and light-curve modeling time for a single band from ~1 day to ≤40 minutes.

We followed the high-precision differential aperture photometry algorithm first proposed by Everett & Howell (2001) and developed by Croll et al. (2010); Wang et al. (2013); Shi et al. (2023) for the data reduction, light curve detrending and final light curve determination. We initially selected all stars with comparable brightness (∆m<2.5 mag) within the FOV as candidate reference stars, preferentially selecting those located within the same detector where the target star resides. To ensure photometric coherence, we further excluded those candidates exhibiting the out-of-eclipse light curves not similar to those of WASP-33 and other references (e.g., Reference Star 4 in Fig. 4), as such discrepancies suggest the star may be intrinsically variable and thus unsuitable as a reference. Unfortunately, WASP-33 is so bright that there are only a few (≤8) qualified candidates for each filter.

For all the selected candidate reference stars, the raw light curves for a given photometry aperture diameter (D) were derived and normalized using the median values derived from a 2-3σ clipping algorithm. Then, we performed systematic quality screening on the light curves and time-series data points. As exemplified by the CO filter data in Fig. 4, the first 1700 s data (phase 0.410-0.426) exhibited degraded quality and were therefore excluded. Ref 4 candidate displayed obviously divergent photometric trends compared to the target star and other reference stars, and was thus discarded from the reference star list. A similar approach was applied to the CH4on filter data, with slightly modified criteria: the initial and final 960s data (phase 0.403-0.412 and 0.564-0.573) were removed due to significantly worse data quality. These two excisions somewhat reduced the amount of out-of-eclipse baseline data, but improved the overall data quality and reduced the scatter.

Then, an average reference light curve for a given aperture D and reference star group (RSG), denoted Lref(D,RSG), was obtained by taking the weighted mean of the normalized light curves with weights proportional to the inverse variance of their photometric uncertainties (wi∝1/σi2Mathematical equation: $w_{i} \propto 1/\sigma_{i}^{2}$); this averaged light curve was then used to divide the target light curve so that most systematics were removed. Finally, a 3σ clipping algorithm was applied to remove obvious outliers in the normalized target light curve.

The derived target light curves were modeled using batman (Kreidberg 2015), with the modeling parameters summarized in Table 1. To provide a reasonable reference for the eclipse depths in the CO and CH4on filters, we built a forward model of WASP-33 b’s dayside atmosphere with petitRADTRANS (Mollière et al. 2019) using the retrieval results of Finnerty et al. (2023), and then compared the modeled planetary spectrum to the host star’s stellar spectrum; accounting for the filters’ bandwidths and transmission profiles with the species (Stolker et al. 2020) yielded the reference eclipse depths. Although background-related polynomial baselines of the NIR secondary eclipse light curves have been widely reported (Croll et al. 2010; Wang et al. 2013; Chen et al. 2014; Shi et al. 2023), we did not detect such features in our data. This absence can probably be attributed to the dominant influence of WASP-33’s pulsations, which may significantly overshadow instrumental or background effects.

Thumbnail: Fig. 4 Refer to the following caption and surrounding text. Fig. 4

Normalized light curves of target star and the candidate reference stars.
Table 1

Parameters for the WASP-33 b system.

3.3 Finding the “best” light curve

It is well recognized in our previous work and the literature that for almost every combination of D and RSG, the obtained light curve can be fitted with an eclipse-wise light curve plus various noises, with slightly different eclipse depths, uncertainties and best-fitting residuals. This highlights the need to find the “best” combination of D and RSG, which gives the “best” light curve that has the smallest fitting residuals with no noticeable hint of overfitting. To do so, we follow the methodology described in Shi et al. (2023), developed from Croll et al. (2015). The key idea of this method is to find a balance between underfitting and overfitting, by combining the two parameters, rms and β, where rms is the root-mean-square of the residual light curve after subtracting the best-fit model, and β is the ratio of the residuals to the Gaussian noise expectation as defined in Winn et al. (2008).The method to derive β is described in detail in Sect. 3.5.

When the photometric aperture D encompasses a region where stellar photons dominate the sky background, the residual rms tends to decrease with increasing D, gradually approaching the optimal aperture. However, as D continues to grow and the aperture begins to include the outer regions of the stellar profile, the red noise becomes increasingly nonnegligible. This leads to a flattening of the rms curve near its minimum or to a slow upward trend. In such cases, the β factor amplifies subtle variations in rms, helping to distinguish and exclude aperture choices with comparable rms values that are affected by correlated noise. Therefore, we use rms and β as proxies to achieve a balance between the two competing parameters. In practice, for each set of RSG and D, we derive the differential target-to-reference-assembly light curve. The next step is to model all (if possible) systematics (or correlated noise) embedding in the derived light curve, and remove them, to achieve a reliable and high-precision eclipse and the minimum rms and β. After that, we obtain an rms*β2 2D grid, which is smoothed, and thus a minimum value of rms*β2 can be located.

In the case of WASP-33 b, an additional approach has to be performed, which is to model and remove the relatively strong stellar pulsation observed in WASP-33. Fortunately, the pattern of the stellar pulsation seems to be irrelevant to the selection of RS G and D, and therefore is left to be corrected afterward.

In Fig. 5, an rms*β2 contour map for the CO filter as functions of D and NRSG is shown, where the minimum is achieved at D = 29 pixels and NRSG = 3. We use the same aperture size for the target and reference stars at this stage, but then perform subsequent refinement by individually varying D for each reference star and the target with ±1 pixel increments, yielding final apertures of D = 29 pixels for the science target and D =28,39,39 pixels for the three reference stars.

Similar approaches have been applied for the CH4on filter data. However, WASP-33’s is so bright in this band that there are only a few usable reference stars (Fig. 6), resulting in the optimal NRSG = 1 and D = 33 pixels. The final apertures chosen are D = 29 and 30 , pixel numbers for the target and reference star, respectively. It is worth noting that in both the CO and CH4on band observations, there seem to be instrumental anomalies in some nondestructive readouts, in which all the 12 sub-integrations show 5σ deviations from the nearby readouts; these outliers are flagged and excluded from model fitting.

As noted in Sect. 1, the WASP-33 system may be a triplestar system (Mugrauer 2019). Although the M dwarf companion (WASP-33 B) can be spatially resolved by a 4m class telescope, it is complicated for our WIRCam observation due to the defocus setup and thus a widespread PSF. Therefore, partial contamination may occur. We therefore carefully examined the raw images in both filters but found no hint of WASP-33 B near its expected position, possibly due to the large magnitude difference of ~6.11 mag in Ks and the very short exposure time. Following Shi et al. (2023), we quantified the M dwarf’s contaminating signal and found that it contributes only a few parts per million (ppm) to the eclipse depth ~ far below the measurement uncertainty. The other even fainter companion is quite far away with an angular separation of ~49.0″ and therefore can be ignored without any risk. We therefore find no evidence for companion-induced contamination in our photometric data and no contamination correction was applied.

Thumbnail: Fig. 5 Refer to the following caption and surrounding text. Fig. 5

Contour maps of the normalized rms*β2 distribution are superimposed on the two-dimensional parameter space defined by D and NRSG in the CO filter data. The pixel scale is 0.306 arcsec pixel−1. The minimum value of rms*β2 is reached with D = 29 and NRSG = 3.
Thumbnail: Fig. 6 Refer to the following caption and surrounding text. Fig. 6

Same as Fig. 5 for the CH4on filter. The minimum value of rms*β2 is reached with D = 33 and NRSG = 1.

3.4 Suppression of stellar pulsation

WASP-33 b is a rare exoplanet orbiting a variable star. The host star HD 15082, a δ Scuti type star, exhibits a dominant pulsation with a period of ~20 days−1 and a maximum amplitude of ~1.5 mmag (von Essen et al. 2014). This pulsation may significantly compromise the accuracy of eclipse depth measurements and therefore must be modeled and removed.

By analyzing the TESS Sector 18 data, Von Essen et al. (2020) identified 29 stellar pulsation frequencies with signal-tonoise ratios (S/N)>4 in WASP-33. However, Baluev & Sokov (2025) demonstrated that the stellar pulsations may evolve on multi-year timescales after including the TESS Sector 58 data. We tried using the pulsation models derived from Baluev & Sokov (2025) and Von Essen et al. (2020) to mitigate the notable photometric variations in the obtained light curves, but the remaining residuals are still nonnegligible and bear some patterns. Deming et al. (2012) were able to mitigate the stellar pulsations by combining long-baseline ground-based observations with nearly contemporaneous Spitzer data, enabling the identification of relatively stable pulsation periods that could be modeled and removed from the light curves. However, in our case, the available baseline is considerably shorter, the CFHT and Spitzer observations are widely separated in time, and the pulsation properties of WASP-33 are known to vary. Under these circumstances, applying a similar approach would likely require fitting a large number of additional parameters and may therefore not be practical.

Therefore, we follow Zhang’s method (Zhang et al. 2018) by using a Gaussian process (GP) to fit these pulsations non-parametrically with the public GP code celerite2, (Foreman-Mackey et al. 2017; Foreman-Mackey 2018). The core of this method is that GP can treat pulsations as a combination of various noises whose properties can be described by a parameterized covariance matrix fitted to the data. By removing the requirement to prescribe a mathematical form to the oscillations, this method enables nonstationary oscillatory behavior throughout the observations. The celerite framework constructs the covariance matrix through a kernel function exclusively parameterized by temporal separation |τ = ti - tj|. The kernel used here consists of two radial functions, inspired by a simple harmonic oscillator, augmented with a diagonal component to account for the white noise, which is k(τ)=S0ω0Qe−ω0τ2Q(cos(ηω0τ)+12ηQsin(ηω0τ)),]Mathematical equation: k (\tau) = S_{0}\omega_{0} Q e^{-\frac{\omega_{0}\tau}{2Q}}(cos(\eta\omega_{0}\tau)+\frac{1}{2 \eta Q}sin(\eta\omega_{0}\tau)) \, ,(1)

where η=|1−(4Q2)−1|12Mathematical equation: $\eta = |1 - (4Q^2)^{-1}|^{\frac{1}{2}}$.

Two sets of kernels are combined to model the stellar pulsations. In the first kernel, all three parameters (Q, S0,w0) are allowed to vary in fit, while in the second kernel, Q is fixed to 1/√2, and only S0,w0 are free parameters. These two kernels represent the oscillatory component and the nonoscillator component that decays rapidly with τ of stellar pulsations, respectively. In addition, we incorporate the eclipse depth as another free parameter in the fitting process to mitigate potential covariance with the forward model assumptions. For the CH4on band data, the broader bandwidth returns a higher S/N, which enables quite robust identification and systematic removal of pulsation-induced noise, as illustrated in Fig. 7. Conversely, the pulsation correction in the CO band is less successful, due to the narrower bandwidth and limited data volume, and thus suppressed pulsation signatures.

3.5 Light curve rebinning and fitting

Until now, the photometric precision of each data point of the obtained light curve is ~10−3, the same order of magnitude as the eclipse depth of a typical HJ. To obtain a robust determination of the eclipse depth, rebinning of data points before light curve fitting is essential. Although the rms*β2 method can identify an optimal bin size together with an optimal Lref (RS G, D), this initial binning is on the light curves obtained using identical D for both the science target and the reference stars, and they still remain contaminated by pulsation noise. After finalizing individual D values for every star and removing pulsation contamination, the bin size for each band is reoptimized.

As described previously, we employ the β parameter (Winn et al. 2008) to quantify the level of correlated noise. Taking the CO filter data as an example, β is calculated by averaging ratios across bin sizes ranging from 2 to 40 data points (with upper limit set to ensure time sampling size to be less than WASP-33 b’s ingress and egress; see Fig. 8). The unphysical cases where the scaled-down residuals fall below Gaussian noise expectations, i.e., β < 1, are discarded. After this evaluation, a bin size of 14 is adopted for the CO filter data. A similar procedure yields a bin size of 18 for the CH4on filter data.

The rebinned CO and CH4on band light curves, as shown in the middle panels of Figs 9 and 10, are used to determine final eclipse depths via Markov chain Monte Carlo (MCMC) analysis using the emcee package (Foreman-Mackey et al. 2013). The MCMC sampling is performed on two parameters, the eclipse depth distributed uniformly, and the central eclipse time Tmid distributed normally with sigma being the error of Tmid, propagated from the uncertainty quoted in Zhang et al. (2018). Each parameter is sampled with six walkers. An initial 1000-step run refines parameter guesses, followed by a final 6000-step sampling with 12 walkers and a step size of 10−8 > 50 times the best-fit binned data length to avoid possible autocorrelation. After discarding the first 20% of the chains as burn-in, we confirmed that the acceptance rates remained within 20%-50%, consistent with appropriately chosen step sizes. Convergence was further assessed using the Gelman-Rubin statistic (R), which was found to be close to unity for all parameters. The best-fit eclipse depths in the CO and CH4on bands are 1565.2−237.5+228.6Mathematical equation: $1565.2^{+228.6}_{-237.5}$ ppm and 914.3−57.0+56.1Mathematical equation: $914.3^{+56.1}_{-57.0}$ ppm, respectively (see Fig. 11).

Thumbnail: Fig. 7 Refer to the following caption and surrounding text. Fig. 7

Pulsation noise removal plot for CH4on filter’s data. The upper panel illustrates photometric measurements contaminated by pul-sational noise, depicted as light points (unbinned) and dark points (binned). A synthesized baseline model (red curve) combines the eclipse signal generated with batman and the pulsational signature predicted by celerite2. The middle and lower panels respectively present the extracted pulsational component and the corrected eclipse light curve after systematic noise removal.
Thumbnail: Fig. 8 Refer to the following caption and surrounding text. Fig. 8

Comparison of data-model residuals and theoretical noise levels across varying bin sizes in CO filter.
Thumbnail: Fig. 9 Refer to the following caption and surrounding text. Fig. 9

Secondary eclipse of WASP-33 b observed by CFHT/WIRCam in the CO filter on October 25, 2015. The top panel displays the unbinned (light-shaded points) and binned (dark-shaded points) light curves overlaid with the best-fit secondary eclipse model (solid red line). The middle panel isolates the binned data alongside the optimized eclipse model, while the bottom panel presents the residuals relative to the model fit. The rms of the binned residuals is 544.9 ppm.
Thumbnail: Fig. 10 Refer to the following caption and surrounding text. Fig. 10

Similar to Fig. 9 but for the secondary eclipse of WASP-33 b observed by CFHT/WIRCam with the CH4on filter on December 24, 2015. The rms of the binned residuals is 217.1 ppm.
Thumbnail: Fig. 11 Refer to the following caption and surrounding text. Fig. 11

Best-fit parameters of the binned light curve for the CO (Upper) and CH4on (Lower), with Tmid = Tc - 57319.75524 and Tmid = Tc - 57381.24756.

3.6 Modeling

In addition to our newly obtained CFHT/WIRCam CO and CH4on band measurements, for the dayside atmosphere investigation, we collected the HST/WFC3 emission spectrum obtained by Haynes et al. (2015) and the Spitzer/IRAC 3.6 and 4.5 μm photometric eclipse depths obtained by Zhang et al. (2018); Dang et al. (2024). In total, we have 13 eclipse depth measurements spanning from 1.1 μm to 4.5 μm. Note that two HST data points are excluded, as explained in Sect. 4.1.

We use the petitRADTRANS package (Mollière et al. 2019; Nasedkin et al. 2024) and construct two atmospheric models: a free-chemistry (FREE) model and an equilibrium chemistry (EQ) model. We use the python package PyMultiNest to explore the posterior distribution and calculate Bayesian evidence (Buchner et al. 2014), which implements the multimodal nested sampling algorithm based on the MultiNest library (Feroz et al. 2009).

For the FREE model, we assume a one-dimensional atmospheric structure with a temperature inversion layer spanning a pressure range of 10−6 to 103 bar, divided into 100 uniformly spaced layers in logarithmic space. Previous HR study indicates that the temperature of the upper atmosphere is lower than that of the inversion layer (Haynes et al. 2015; Nugroho et al. 2017; Herman et al. 2020; Finnerty et al. 2023), which means that the commonly adopted 2-point temperature-pressure (T-P) profile (Buchner et al. 2014) is inadequate. We therefore introduce a third point (T3,P3). Unlike the 3-point T-P profile proposed by Waldmann et al. (2015), we fix P3 at 10−6 bar rather than treating it as a free parameter. Based on previous studies on WASP-33 b’s dayside, we include CO (Rothman et al. 2010), H2O (Rothman et al. 2010), CH4 (Yurchenko et al. 2017), and TiO (McKemmish et al. 2019), which possess prominent spectral features in the wavelength coverage of interests. We have not included OH - a photodissociation product of H2O - in our model, due to the lack of OH-sensitive features in the data set we use in this work. The mass fractions of these species are treated as free parameters, while H2 and He serve as filler gases to maintain a total mass fraction, with a fixed He/H2 mass ratio of 0.305. Contributions from Rayleigh scattering by H2/He and continuum absorption from H2-H2, H2-He, and H− are considered, and the inclusion of H− necessitates treating the mass fractions of its associated species, e− and H, as free parameters. Given the high temperatures on the dayside, the atmosphere is likely to be cloud-free, and thus no cloud prescription is included in the retrieval process. For the EQ model, we use the same 3-point T-P profile as described above. The free parameters to be explored are the C/O ratio and the metallicity [Fe/H]. Both retrieval models are initially computed at a spectral resolution (λ/∆λ) of 1000 and subsequently binned to match the observational passbands. They share common fixed parameters, including planetary radius of 1.593 RJ, surface gravity of 27.114 m/s2, and orbit semi-major axis of 0.02558 au.

To identify potential species from the emission spectrum obtained, we first computed the Bayesian evidence (Kass & Raftery 1995) for the full model (including all species). We then iteratively excluded one species at a time and recalculated the Bayes factors for each reduced model. Following the criteria proposed by Kass & Raftery (1995), the strength of evidence is categorized as strong (|∆lnZ| ≥ 5), moderate (3 ≤ |∆lnZ| < 5), weak (1 ≤ |∆lnZ| < 3), or inconclusive (|∆lnZ| ≤ 1.

4 Retrieval analysis

4.1 FREE model retrieval

The FREE retrieval analysis yields a best-fit spectrum that is largely featureless, with residuals exceeding 3σ at the blue-end datapoints (1.155 and 1.199 μm). These two points dominate the χ2 value, leading to a generally poor fit.

Haynes et al. (2015) attributed the unusually large depths at 1.155 and 1.199 μm to a high abundance of TiO, but subsequent joint retrievals (Changeat et al. 2022) combining HST and Spitzer data and later HRS studies have contested the presence of dayside TiO (Herman et al. 2020; Cont et al. 2021). To test this hypothesis, we ran several forward models with temperature inversion layers that included only H, He, and TiO, assuming different TiO abundances, and found that the TiO abundance must be ≥10−4 to account for the depths observed 1.155 and 1.199 μm. However, such a high TiO abundance should lead the eclipse depths of the CH4on and CO bands (especially the latter) to be significantly higher than our measured depths, and TiO should be detected in HR studies (Nugroho et al. 2017; Herman et al. 2020; Cont et al. 2021; Serindag et al. 2021).

We therefore argue that the large depths of 1.155 and 1.199 μm cannot be due to a high TiO abundance; rather, they come from their imperfect correction of the pulsation signal. Haynes et al. (2015) noticed the pulsation signals in their light curves and modeled the pulsation signals using sine functions. However, such a procedure was later found to only be capable of removing ~10% of pulsation noise (Zhang et al. 2018), thus leaving substantial pulsational residuals in at least some of the HST/WFC3 light curves. Given that WASP-33 is an A-type star, pulsational noise affects the precision of eclipse depth measurements much more significantly at shorter wavelengths, and thus the 1.155 and 1.199 μm depths should still be affected by pulsation noise. To mitigate this and given that the two bluest datapoints are inconsistent with the other datapoints, we excluded the 1.155 and 1.199 μm datapoints in the subsequent retrieval analysis.

In addition to HST measurements, two ground-based secondary-eclipse measurements of WASP-33 b have been reported in the literature: a 0.9 μm datapoint from Smith et al. (2011) and a Ks-band measurement from Deming et al. (2012). The 0.9 μm measurement is not included in our joint retrieval because the transmission curve of the filter used in that observation is not available, preventing self-consistent modeling within petitRADTRANS. In contrast, Ks-band eclipse measurement can be readily incorporated. We therefore performed additional retrievals including this datapoint and found that the inferred atmospheric parameters remained fully consistent with our nominal results, with all differences well within the 1σ uncertainties. This demonstrates that our conclusions are insensitive to the inclusion of the Ks-band measurement; consequently, it is not included in the final retrieval presented here.

The best-fit result for the FREE model is shown in Fig. 12, with Fig. 13 & Table 2 illustrating the corresponding mass mixing ratios (MMRs) and T-P profile. Table 3 reports the detection significance represented by Bayesian evidence for each chemical species.

In the FREE retrieval, WASP-33 b’s dayside atmosphere exhibits a strong temperature inversion layer and a featureless emission spectrum, resulting in no or very poor constraints on the MMRs for all species considered except H− and e−, which have quite high concentrations. The lack of concrete detection of molecules is in line with the low resolution (LR) spacebased results from Changeat et al. (2022), but seems to conflict with the HR results (Finnerty et al. 2023). The results may be caused by: (1) high temperatures in the emission-contribution layers that thermally dissociate most molecules and thus reduce their detectability; (2) residual pulsation noise in the HST spectra, which could still obscure molecular features, even after correction; and (3) the featureless spectrum forcing the retrieval to fit continuum opacity by the inclusion of a large amount of H−.

EQ retrieval analysis yields consistent results with those from the HR spectroscopic study (Changeat et al. 2022). The best-fit model spectrum, the corresponding T-P profile, and the pressure-dependent chemical abundances are presented in Figs. 14, 15 and Table 4, respectively. The retrieved T-P profile suggests that the temperature inversion layer initiates near 1 bar, consistent with previous results (e.g. Nugroho et al. 2021 ; van Sluijs et al. 2023; Finnerty et al. 2023), although the layer above the inversion is relatively cooler in our results. As shown in Table 4, the dayside atmosphere of WASP-33 b should be highly metal-rich and exhibit a super-solar C/O ratio.

Table 5 summarizes the Bayesian-evidence validation for each species, showing strong support for H−, H2O, and CO, weak support for TiO, and no support for CH4. Figure 16 displays the MMRs of these species as a function of the pressure derived from the retrieved T-P profile. Given that CH4 is disfavored by Bayesian evidence and is thermochemically implausible in the dayside atmosphere of UHJs, we therefore do not further examine the CH4 MMR. Although H− abundance in the EQ model is appreciably lower than in the FREE model, the strong Bayesian evidence for H− in both frameworks underscores its necessity in the dayside atmosphere. The very high and nearly constant MMR of CO with altitude, together with the comparatively lower abundance of H2O, is consistent with the Keck Planet Imager and Characterizer (KPIC) detection reported by Finnerty et al. (2023). The low concentration of TiO is likely a natural consequence of thermal dissociation in the hot upper atmosphere. Note that the retrieved emission spectrum is dominated by atmospheric layers between 1 μbar and 1 bar; therefore, the species abundances within this pressure range are expected to be the most reliable.

Table 2

Best-fit parameters from the FREE model retrieval on WASP-33 b.

Thumbnail: Fig. 12 Refer to the following caption and surrounding text. Fig. 12

Best-fit FREE model spectrum together with the data points from HST, CFHT, and Spitzer shown in filled dots with error bars and distinct colors. The orange curve and surrounding shaded region represent the best-fit model and its 3σ regime. The shaded areas in different colors below the orange curve indicate the reference models containing only H2O, CO, TiO, or H−, respectively.
Thumbnail: Fig. 13 Refer to the following caption and surrounding text. Fig. 13

Retrieved temperature-pressure (T-P) profile for the FREE model. The solid blue line shows the best-fit solution, and the shaded area denotes the 1σ confidence interval derived from the posterior distribution.
Table 3

Bayesian model comparison of the chemical species considered for the FREE model.

4.2 EQ model retrieval

In conventional planet formation theories, such as pebble accretion (Ormel & Klahr 2010; Lambrechts & Johansen 2012), metallicity and the C/O ratio generally show an inverse correlation, making it unlikely for planets to simultaneously possess high C/O ratios and high metallicity. However, the pebble drift theory proposed by Booth et al. (2017) indicates that the metal-licity of planets can be enriched not only by accreting solids but also by accreting metal-rich gas generated through the pebble drift mechanism. Synthesizing the dynamical and chemical properties of WASP-33 b, our conclusions align with those outlined in Finnerty et al. (2023): the planet likely formed in a carbon- and solid-rich zone with high C/O ratio, such as near the CO2 snow line (~ 10 au for an A-type primary star; Öberg et al. 2011; Mollière et al. 2022), accreted a high-metallicity atmosphere from its vicinity and migrated to its present location via a process like eccentric Lidov-Kozai effect (Naoz et al. 2011). It should be noted, as suggested by Finnerty et al. (2023), that the observed high C/O ratio in the dayside atmosphere could also be biased due to data interpretation, namely by weak constraints on oxygen-bearing species. For example, our retrieval lacks constraints on OH, a photodissociation product of H2O. Given that a significant amount of H2O should be photodissociated into OH in upper atmospheres, this lack may lead to an overestimated C/O ratio.

On the other hand, using only metallicity and the C/O ratio as diagnostics is insufficient to accurately or unbiasedly constrain planet formation and evolution history (Feinstein et al. 2025). Incorporating measurements of additional tracers (e.g., volatile elements and refractory element abundances) would better enrich the understanding of WASP-33 b’s formation history. We recommend that future JWST observations of WASP-33 b achieve more precise atmospheric constraints.

Thumbnail: Fig. 14 Refer to the following caption and surrounding text. Fig. 14

Same as Fig. 12, but for the EQ retrieval. Note that H2O, H , and CO all exhibit prominent features in the model spectrum.
Thumbnail: Fig. 15 Refer to the following caption and surrounding text. Fig. 15

Same as Fig. 13 but for the retrieved T-P profile from the EQ model.
Table 4

Best-fit parameters from the EQ model retrieval on WASP-33 b.

Table 5

Bayesian model comparison of the chemical species retrieved with EQ model.

Thumbnail: Fig. 16 Refer to the following caption and surrounding text. Fig. 16

Mass mixing ratios of chemical species under EQ model as functions of atmospheric height. Apart from CO maintaining high concentrations and H exhibiting a marked increase with rising temperature, all other molecules undergo thermal dissociation within the thermal inversion layer.

5 Summary and conclusion

We conducted two secondary-eclipse observations of the ultrahot Jupiter WASP-33 b with CFHT/WIRCam to constrain the planet’s dayside atmosphere. To handle the thousands of iterations required for systematic correction and atmospheric retrievals, we developed a parallelized data reduction and analysis pipeline that substantially reduced computational time. After carefully removing the host star’s pulsation noise, we measured eclipse depths of 1565.2−237.5+228.6Mathematical equation: $1565.2^{+228.6}_{-237.5}$ ppm in the CO filter and 914.3−57.0+56.1Mathematical equation: $914.3^{+56.1}_{-57.0}$ ppm in the CH4on filter.

In addition to our CFHT/WIRCam observations, we incorporated HST/WFC3 spectroscopic data (Haynes et al. 2015) and Spitzer/IRAC photometric measurements - 3.6 μm (Zhang et al. 2018) and 4.5 μm (Dang et al. 2024) - to jointly perform both free-chemistry and equilibrium-chemistry atmospheric retrievals using petitRADTRANS. Given the extreme dayside temperatures typical of ultra-hot Jupiters, we consider the equilibriumchemistry model to be a more physically representative description of WASP-33 b’s dayside atmosphere, compared to the free chemistry model. Accordingly, all results presented in the main findings in the following are based on the equilibrium model. Our main findings are summarized below.

  1. Thermal inversion and opacity sources. WASP-33 b exhibits a clear thermal inversion from both free and equilibrium retrievals. Under the free-chemistry model, the dayside emission spectrum appears to be largely featureless, with Bayesian evidence strongly favoring the presence of H− opacity. However, the retrieved abundances of H− and e− are considered unreliable due to model limitations and spectral coverage. Incorporating blue or optical spectra in future joint retrievals will better constrain these species.

  2. High C/O ratio and metallicity in equilibrium chemistry. The equilibrium-chemistry retrieval yields a high C/O ratio of 0.78−0.04+0.03Mathematical equation: $0.78^{+0.03}_{-0.04}$ and a metallicity of ~26× stellar, both more precisely constrained than in previous studies. The emission spectrum is primarily shaped by H2O, CO, and H−, while the existence of TiO remains uncertain due to weak Bayesian evidence. In the upper atmosphere, most molecules undergo thermal dissociation except CO.

  3. Implications for formation and migration. The combination of high metallicity and high C/O ratio is consistent with planetary enrichment through the accretion of metal-rich gas, as proposed by the pebble drift theory (Booth et al. 2017). A possible formation and migration scenario suggests that WASP-33 b formed and/or accreted material in a carbon- and solid-rich region - possibly near the CO2 snow line - before migrating inward via a mechanism such as the eccentric Lidov-Kozai effect (Naoz et al. 2011). To test this hypothesis and rule out apparent carbon enrichment caused by weak constraints on oxygenbearing species, we recommend follow-up JWST observations to obtain tighter constraints on additional elemental abundances. WASP-33 b remains an unusual exoplanet, orbiting a δ Scuti variable star. Its high metallicity, elevated C/O ratio, and misaligned orbit all indicate a distinct formation and evolutionary history. Although stellar pulsations previously posed major challenges for measuring transit and eclipse depths, improved methodologies now allow for better mitigation of these signals and stronger isolation of the planetary component. Future observations-either at higher spectral resolution in targeted bands or with broader wavelength coverage at lower resolution - will enable more detailed characterization of WASP-33 b’s atmosphere and provide a deeper insight into its formation and evolution.

Acknowledgements

We thank the referee and editor for their concise review and instructive suggestion. This research is supported by the National Key R&D Program of China (2025YFE0102100, 2024YFA1611802, 2025YFE0213204), the National Natural Science Foundation of China under grant 62127901, 12588202, National Astronomical Observatories Chinese Academy of Sciences No. E4TQ2101, and the Pre-research project on Civil Aerospace Technologies No. D010301 funded by the China National Space Administration (CNSA), by the China Manned Space Program with grant no. CMS-CSST-2025-A16. We acknowledge financial support from the Agencia Estatal de Investigación of the Ministerio de Ciencia e Innovación MCIN/AEI/10.13039/501100011033 and the ERDF ’A way of making Europe’ through projects PID2021-125627OB-C3 2 and PID2024-158486OB-C32. Q.Y. Zou is supported by the Program of China Scholarship Council (Grant CSC202510740003). MZ was supported by the Chinese Academy of Sciences (CAS), through a grant to the CAS South America Center for Astronomy (CASSACA) in Santiago, Chile. This research uses data obtained through the Telescope Acess Program(TAP), which has been funded by the National Astronomical Observatories of China, the Chinese Academy of Sciences, and the Special Fund for Astronomy from the Ministry of Finance. Two nights (Oct. 25 and Dec. 24, 2015) with WIRCam on the 3.6 m CFHT telescope were distributed to us for scientific studies of exoplanetary atmosphere via the TAP. We also thank Dr. Jie Zheng, Jingyuan Zhao, Hengkai Ding and Yingjie Cai for their constructive comments on this work.

References

  1. Anderson, D. R., Temple, L. Y., Nielsen, L. D., et al. 2018, arXiv e-prints [arXiv:1809.04897] [Google Scholar]
  2. Baluev, R. V., & Sokov, E. N. 2025, MNRAS, 537, 171 [Google Scholar]
  3. Baxter, C., Désert, J.-M., Parmentier, V., et al. 2020, A&A, 639, A36 [EDP Sciences] [Google Scholar]
  4. Booth, R. A., Clarke, C. J., Madhusudhan, N., & Ilee, J. D. 2017, MNRAS, 469, 3994 [Google Scholar]
  5. Bradley, L., Sipocz, B., Robitaille, T., et al. 2024, astropy/photutils: 1.11.0 [Google Scholar]
  6. Buchner, J., Georgakakis, A., Nandra, K., et al. 2014, A&A, 564, A125 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  7. Cabot, S. H., Madhusudhan, N., Welbanks, L., Piette, A., & Gandhi, S. 2020, MNRAS, 494, 363 [NASA ADS] [CrossRef] [Google Scholar]
  8. Cabot, S. H., Bello-Arufe, A., Mendonça, J. M., et al. 2021, AJ, 162, 218 [NASA ADS] [CrossRef] [Google Scholar]
  9. Chakrabarty, A., & Sengupta, S. 2019, AJ, 158, 39 [NASA ADS] [CrossRef] [Google Scholar]
  10. Changeat, Q., Edwards, B., Al-Refaie, A. F., et al. 2022, ApJS, 260, 3 [NASA ADS] [CrossRef] [Google Scholar]
  11. Chen, G., van Boekel, R., Wang, H., et al. 2014, A&A, 563, A40 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  12. Christian, D. J., Pollacco, D., Skillen, I., et al. 2006, MNRAS, 372, 1117 [NASA ADS] [CrossRef] [Google Scholar]
  13. Cont, D., Yan, F., Reiners, A., et al. 2021, A&A, 651, A33 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  14. Cont, D., Yan, F., Reiners, A., et al. 2022, A&A, 668, A53 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  15. Cont, D., Yan, F., Reiners, A., et al. 2022, A&A, 657, L2 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  16. Cridland, A. J., Eistrup, C., & van Dishoeck, E. F. 2019, A&A, 627, A127 [Google Scholar]
  17. Croll, B., Albert, L., Lafreniere, D., Jayawardhana, R., & Fortney, J. J. 2010, ApJ, 717, 1084 [Google Scholar]
  18. Croll, B., Albert, L., Jayawardhana, R., et al. 2015, ApJ, 802, 28 [Google Scholar]
  19. Cutri, R., Skrutskie, M., Van Dyk, S., et al. 2003, VizieR Online Data Catalog: 2246, II [Google Scholar]
  20. Dang, L., Bell, T. J., Shu, Y. Z., et al. 2024, AJ, 169, 32 [Google Scholar]
  21. De Mooij, E., Brogi, M., De Kok, R., et al. 2013, A&A, 550, A54 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  22. Delrez, L., Santerne, A., Almenara, J. M., et al. 2016, MNRAS, 458, 4025 [NASA ADS] [CrossRef] [Google Scholar]
  23. Deming, D., Fraine, J. D., Sada, P. V., et al. 2012, ApJ, 754, 106 [NASA ADS] [CrossRef] [Google Scholar]
  24. Devost, D., Albert, L., Teeple, D., & Croll, B. 2010, in SPIECS, 7737, Observatory Operations: Strategies, Processes, and Systems III, eds. D. R. Silva, A. B. Peck, & B. T. Soifer, 77372D [Google Scholar]
  25. Ding, H., Shu, Y., Chen, Y., et al. 2025, RAA, 25, 065013 [Google Scholar]
  26. Drummond, B., Carter, A. L., Hébrard, E., et al. 2019, MNRAS, 486, 1123 [NASA ADS] [CrossRef] [Google Scholar]
  27. Espinoza, N., Fortney, J. J., Miguel, Y., Thorngren, D., & Murray-Clay, R. 2017, ApJ, 838, L9 [Google Scholar]
  28. Evans-Soma, T. M., Sing, D. K., Barstow, J. K., et al. 2025, Nat. Astron., 1 [Google Scholar]
  29. Everett, M. E., & Howell, S. B. 2001, PASP, 113, 1428 [Google Scholar]
  30. Feinstein, A. D., Booth, R. A., Bergner, J. B., et al. 2025, arXiv preprint [arXiv:2506.00669] [Google Scholar]
  31. Feroz, F., Hobson, M., & Bridges, M. 2009, MNRAS, 398, 1601 [NASA ADS] [CrossRef] [Google Scholar]
  32. Finger, G., Dorn, R. J., Eschbaumer, S., et al. 2008, in High Energy, Optical, and Infrared Detectors for Astronomy III, 7021, SPIE, 236 [Google Scholar]
  33. Finnerty, L., Schofield, T., Sappey, B., et al. 2023, AJ, 166, 31 [NASA ADS] [CrossRef] [Google Scholar]
  34. Foreman-Mackey, D. 2018, RNAAS, 2, 31 [NASA ADS] [Google Scholar]
  35. Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306 [Google Scholar]
  36. Foreman-Mackey, D., Agol, E., Ambikasaran, S., & Angus, R. 2017, AJ, 154, 220 [Google Scholar]
  37. Gaia Collaboration 2020, VizieR Online Data Catalog: 1350, I [Google Scholar]
  38. Gandhi, S., Kesseli, A., Zhang, Y., et al. 2023, AJ, 165, 242 [NASA ADS] [CrossRef] [Google Scholar]
  39. Gaudi, B. S., Stassun, K. G., Collins, K. A., et al. 2017, Nature, 546, 514 [NASA ADS] [Google Scholar]
  40. Grenier, S., Baylac, M.-O., Rolland, L., et al. 1999, A&AS, 137, 451 [Google Scholar]
  41. Haynes, K., Mandell, A. M., Madhusudhan, N., Deming, D., & Knutson, H. 2015, ApJ, 806, 146 [NASA ADS] [CrossRef] [Google Scholar]
  42. Herman, M. K., de Mooij, E. J. W., Jayawardhana, R., & Brogi, M. 2020, AJ, 160, 93 [Google Scholar]
  43. Herman, M. K., de Mooij, E. J., Nugroho, S. K., Gibson, N. P., & Jayawardhana, R. 2022, AJ, 163, 248 [NASA ADS] [CrossRef] [Google Scholar]
  44. Hoeijmakers, H., Kitzmann, D., Morris, B., et al. 2024, A&A, 685, A139 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  45. Johnson, M. C., Cochran, W. D., Collier Cameron, A., & Bayliss, D. 2015, ApJ, 810, L23 [Google Scholar]
  46. Kass, R. E., & Raftery, A. E. 1995, JASA, 90, 773 [CrossRef] [Google Scholar]
  47. Kreidberg, L. 2015, PASP, 127, 1161 [Google Scholar]
  48. Lambrechts, M., & Johansen, A. 2012, A&A, 544, A32 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  49. Lendl, M., Csizmadia, S., Deline, A., et al. 2020, A&A, 643, A94 [EDP Sciences] [Google Scholar]
  50. Madhusudhan, N. 2012, ApJ, 758, 36 [NASA ADS] [CrossRef] [Google Scholar]
  51. McKemmish, L. K., Masseron, T., Hoeijmakers, H. J., et al. 2019, MNRAS, 488, 2836 [Google Scholar]
  52. Mollière, P., van Boekel, R., Dullemond, C., Henning, T., & Mordasini, C. 2015, ApJ, 813, 47 [Google Scholar]
  53. Mollière, P., Wardenier, J., Van Boekel, R., et al. 2019, A&A, 627, A67 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  54. Mollière, P., Molyarova, T., Bitsch, B., et al. 2022, ApJ, 934, 74 [CrossRef] [Google Scholar]
  55. Moses, J., Madhusudhan, N., Visscher, C., & Freedman, R. 2012, ApJ, 763, 25 [Google Scholar]
  56. Moya, A., Bouy, H., Marchis, F., Vicente, B., & Barrado, D. 2011, A&A, 535, A110 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  57. Mraz, G., Darveau-Bernier, A., Boucher, A., et al. 2024, ApJ, 975, L42 [Google Scholar]
  58. Mugrauer, M. 2019, MNRAS, 490, 5088 [Google Scholar]
  59. Naoz, S., Farr, W. M., Lithwick, Y., Rasio, F. A., & Teyssandier, J. 2011, Nature, 473, 187 [NASA ADS] [CrossRef] [Google Scholar]
  60. Nasedkin, E., Mollière, P., & Blain, D. 2024, J. Open Source Softw., 9, 5875 [CrossRef] [Google Scholar]
  61. Ngo, H., Knutson, H. A., Hinkley, S., et al. 2016, ApJ, 827, 8 [Google Scholar]
  62. Nugroho, S. K., Kawahara, H., Masuda, K., et al. 2017, AJ, 154, 221 [Google Scholar]
  63. Nugroho, S. K., Gibson, N. P., de Mooij, E. J. W., et al. 2020, ApJ, 898, L31 [Google Scholar]
  64. Nugroho, S. K., Kawahara, H., Gibson, N. P., et al. 2021, ApJ, 910, L9 [CrossRef] [Google Scholar]
  65. Öberg, K. I., Murray-Clay, R., & Bergin, E. A. 2011, ApJ, 743, L16 [Google Scholar]
  66. Ormel, C., & Klahr, H. 2010, A&A, 520, A43 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  67. Parmentier, V., Line, M. R., Bean, J. L., et al. 2018, A&A, 617, A110 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  68. Pollacco, D. L., Skillen, I., Cameron, A. C., et al. 2006, PASP, 118, 1407 [Google Scholar]
  69. Rothman, L. S., Gordon, I., Barber, R., et al. 2010, JQSRT, 111, 2139 [NASA ADS] [CrossRef] [Google Scholar]
  70. Rubenzahl, R. A., Dai, F., Halverson, S., et al. 2024, AJ, 168, 188 [Google Scholar]
  71. Serindag, D. B., Nugroho, S. K., Mollière, P., et al. 2021, A&A, 645, A90 [EDP Sciences] [Google Scholar]
  72. Sheppard, K. B., Mandell, A. M., Tamburo, P., et al. 2017, ApJ, 850, L32 [Google Scholar]
  73. Shi, Y., Wang, W., Zhao, G., et al. 2023, MNRAS, 522, 1491 [NASA ADS] [CrossRef] [Google Scholar]
  74. Smith, A., Anderson, D., Skillen, I., Cameron, A. C., & Smalley, B. 2011, MNRAS, 416, 2096 [NASA ADS] [CrossRef] [Google Scholar]
  75. Smith, P. C. B., Sanchez, J. A., Line, M. R., et al. 2024, AJ, 168, 293 [NASA ADS] [CrossRef] [Google Scholar]
  76. Snellen, I. A. 2025, ARA&A, 63 [Google Scholar]
  77. Stangret, M., Pallé, E., Casasayas-Barris, N., et al. 2021, A&A, 654, A73 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  78. Stangret, M., Casasayas-Barris, N., Pallé, E., et al. 2022, A&A, 662, A101 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  79. Stassun, K. G., Oelkers, R. J., Paegert, M., et al. 2019, AJ, 158, 138 [Google Scholar]
  80. Stephan, A. P., Wang, J., Cauley, P. W., et al. 2022, ApJ, 931, 111 [NASA ADS] [CrossRef] [Google Scholar]
  81. Stolker, T., Quanz, S. P., Todorov, K. O., et al. 2020, A&A, 635, A182 [EDP Sciences] [Google Scholar]
  82. Tan, X., Komacek, T. D., Batalha, N. E., et al. 2024, MNRAS, 528, 1016 [NASA ADS] [CrossRef] [Google Scholar]
  83. van Sluijs, L., Birkby, J. L., Lothringer, J., et al. 2023, MNRAS, 522, 2145 [NASA ADS] [CrossRef] [Google Scholar]
  84. von Essen, C., Czesla, S., Wolter, U., et al. 2014, A&A, 561, A48 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  85. von Essen, C., Mallonn, M., Welbanks, L., et al. 2019, A&A, 622, A71 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  86. Von Essen, C., Mallonn, M., Borre, C., et al. 2020, A&A, 639, A34 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  87. Waldmann, I. P., Rocchetto, M., Tinetti, G., et al. 2015, ApJ, 813, 13 [Google Scholar]
  88. Wang, W., van Boekel, R., Madhusudhan, N., et al. 2013, ApJ, 770, 70 [NASA ADS] [CrossRef] [Google Scholar]
  89. Watanabe, N., Narita, N., Palle, E., et al. 2022, MNRAS, 512, 4404 [NASA ADS] [CrossRef] [Google Scholar]
  90. Winn, J. N., Holman, M. J., Torres, G., et al. 2008, ApJ, 683, 1076 [Google Scholar]
  91. Wöllert, M., & Brandner, W. 2015, A&A, 579, A129 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  92. Yan, F., & Henning, T. 2018, Nat. Astron., 2, 714 [Google Scholar]
  93. Yan, F., Wyttenbach, A., Casasayas-Barris, N., et al. 2021, A&A, 645, A22 [EDP Sciences] [Google Scholar]
  94. Yan, F., Pallé, E., Reiners, A., et al. 2022, A&A, 661, L6 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  95. Yang, Y., Chen, G., Wang, S., & Yan, F. 2024a, AJ, 167, 36 [CrossRef] [Google Scholar]
  96. Yang, Y., Chen, G., Yan, F., Tan, X., & Ji, J. 2024b, ApJ, 971, L8 [Google Scholar]
  97. Yurchenko, S. N., Amundsen, D. S., Tennyson, J., & Waldmann, I. P. 2017, A&A, 605, A95 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  98. Zhang, M., Knutson, H. A., Kataria, T., et al. 2018, AJ, 155, 83 [Google Scholar]

All Tables

Table 1

Parameters for the WASP-33 b system.

In the text
Table 2

Best-fit parameters from the FREE model retrieval on WASP-33 b.

In the text
Table 3

Bayesian model comparison of the chemical species considered for the FREE model.

In the text
Table 4

Best-fit parameters from the EQ model retrieval on WASP-33 b.

In the text
Table 5

Bayesian model comparison of the chemical species retrieved with EQ model.

In the text

All Figures

Thumbnail: Fig. 1 Refer to the following caption and surrounding text. Fig. 1

Orbit phases of our two observations. Specifically, Night 1 covered phase 0.410-0.589 (October 25, 2015), Night 2 covered phase 0.403-0.573 (December 24, 2015).
In the text
Thumbnail: Fig. 2 Refer to the following caption and surrounding text. Fig. 2

Transmission curves of the CO and CH4on filters used in this work.
In the text
Thumbnail: Fig. 3 Refer to the following caption and surrounding text. Fig. 3

CO filter reduced full-frame WIRCam image of WSAP-33 taken on October 25, 2015. The target star is marked as the white square and the final chosen reference stars are marked as the grey circles.
In the text
Thumbnail: Fig. 4 Refer to the following caption and surrounding text. Fig. 4

Normalized light curves of target star and the candidate reference stars.
In the text
Thumbnail: Fig. 5 Refer to the following caption and surrounding text. Fig. 5

Contour maps of the normalized rms*β2 distribution are superimposed on the two-dimensional parameter space defined by D and NRSG in the CO filter data. The pixel scale is 0.306 arcsec pixel−1. The minimum value of rms*β2 is reached with D = 29 and NRSG = 3.
In the text
Thumbnail: Fig. 6 Refer to the following caption and surrounding text. Fig. 6

Same as Fig. 5 for the CH4on filter. The minimum value of rms*β2 is reached with D = 33 and NRSG = 1.
In the text
Thumbnail: Fig. 7 Refer to the following caption and surrounding text. Fig. 7

Pulsation noise removal plot for CH4on filter’s data. The upper panel illustrates photometric measurements contaminated by pul-sational noise, depicted as light points (unbinned) and dark points (binned). A synthesized baseline model (red curve) combines the eclipse signal generated with batman and the pulsational signature predicted by celerite2. The middle and lower panels respectively present the extracted pulsational component and the corrected eclipse light curve after systematic noise removal.
In the text
Thumbnail: Fig. 8 Refer to the following caption and surrounding text. Fig. 8

Comparison of data-model residuals and theoretical noise levels across varying bin sizes in CO filter.
In the text
Thumbnail: Fig. 9 Refer to the following caption and surrounding text. Fig. 9

Secondary eclipse of WASP-33 b observed by CFHT/WIRCam in the CO filter on October 25, 2015. The top panel displays the unbinned (light-shaded points) and binned (dark-shaded points) light curves overlaid with the best-fit secondary eclipse model (solid red line). The middle panel isolates the binned data alongside the optimized eclipse model, while the bottom panel presents the residuals relative to the model fit. The rms of the binned residuals is 544.9 ppm.
In the text
Thumbnail: Fig. 10 Refer to the following caption and surrounding text. Fig. 10

Similar to Fig. 9 but for the secondary eclipse of WASP-33 b observed by CFHT/WIRCam with the CH4on filter on December 24, 2015. The rms of the binned residuals is 217.1 ppm.
In the text
Thumbnail: Fig. 11 Refer to the following caption and surrounding text. Fig. 11

Best-fit parameters of the binned light curve for the CO (Upper) and CH4on (Lower), with Tmid = Tc - 57319.75524 and Tmid = Tc - 57381.24756.
In the text
Thumbnail: Fig. 12 Refer to the following caption and surrounding text. Fig. 12

Best-fit FREE model spectrum together with the data points from HST, CFHT, and Spitzer shown in filled dots with error bars and distinct colors. The orange curve and surrounding shaded region represent the best-fit model and its 3σ regime. The shaded areas in different colors below the orange curve indicate the reference models containing only H2O, CO, TiO, or H−, respectively.
In the text
Thumbnail: Fig. 13 Refer to the following caption and surrounding text. Fig. 13

Retrieved temperature-pressure (T-P) profile for the FREE model. The solid blue line shows the best-fit solution, and the shaded area denotes the 1σ confidence interval derived from the posterior distribution.
In the text
Thumbnail: Fig. 14 Refer to the following caption and surrounding text. Fig. 14

Same as Fig. 12, but for the EQ retrieval. Note that H2O, H , and CO all exhibit prominent features in the model spectrum.
In the text
Thumbnail: Fig. 15 Refer to the following caption and surrounding text. Fig. 15

Same as Fig. 13 but for the retrieved T-P profile from the EQ model.
In the text
Thumbnail: Fig. 16 Refer to the following caption and surrounding text. Fig. 16

Mass mixing ratios of chemical species under EQ model as functions of atmospheric height. Apart from CO maintaining high concentrations and H exhibiting a marked increase with rising temperature, all other molecules undergo thermal dissociation within the thermal inversion layer.
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.