Keck/KPIC Emission Spectroscopy of WASP-33b

WASP-33 (HD 15082) is a bright ($K$ = 7.47, Cutri et al., 2003), rapidly rotating ($v\sin i=86$ $\rm km\ s^{-1}$ ) A5 star (Collier Cameron et al., 2010; Johnson et al., 2015). WASP-33b was first discovered by Christian et al. (2006) in SuperWASP (Pollacco et al., 2006) transit photometry with a 1.22-day orbital period. Subsequent radial velocity measurements obtained a mass of 2.8 $\rm M_{J}$ (von Essen et al., 2014). The estimated planet mass and host star spectral type are similar to the HR 8799 system (Goździewski & Migaszewski, 2020), though the WASP-33b semimajor axis is more than 10³ times smaller.

As a transiting planet around a rapid rotator, WASP-33b is an ideal target for observing the Rossiter-McLaughlin Effect (RME; Ohta et al., 2005). Observations by Collier Cameron et al. (2010) found that WASP-33b has a near-polar retrograde orbit, with subsequent observations measuring orbital precession and the stellar quadrupole moment (Johnson et al., 2015; Stephan et al., 2022). The extremely small semimajor axis and near-polar orbit suggest the possibility of eccentric Kozai-Lidov dynamical effects (Naoz et al., 2011), though there are presently no additional planets known in the system to act as perturbers. This system architecture and misaligned orbit is extremely similar to the UHJ KELT-9b (Gaudi et al., 2017), which may be a result of a common formation pathway for both systems.

WASP-33b has been a frequent target for previous atmospheric characterization studies. Near-IR dayside observations have detected CO and OH in emission, but found only weak signs of H₂O, indicative of high-temperature thermochemistry (Nugroho et al., 2021; van Sluijs et al., 2022; Yan et al., 2022). Optical observations have detected a number of oxides and atomic species, including AlO (von Essen et al., 2019), Si I (Cont et al., 2022a), Fe I (Cont et al., 2021; Herman et al., 2022), Ti I, and V I (Cont et al., 2022b), while transit observations of Balmer lines suggest an extended escaping hydrogen envelope driving mass loss at a rate of $\sim 10^{12}\rm\ g/s$ (Yan et al., 2021). These observations have focused on dayside longitudes or the terminators, with phase curve observations suggesting inefficient redistribution of heat to the nightside (Zhang et al., 2018; Herman et al., 2022). These findings paint the picture of WASP-33b as an extreme ultra-hot Jupiter, with mass loss, photochemistry, and dramatic day/night differences.

WASP-33 is a rapid rotator, with rotational broadening and gravity darkening effects significantly impacting the stellar spectrum. WASP-33 is also a $\delta$ Scuti pulsator (Herrero et al., 2011). We discuss our treatment of these effects in Section 3.1.3.

We observed WASP-33 with Keck/KPIC (Delorme et al., 2021) on UT dates 2021 November 21 and 2021 November 22. KPIC (Keck Planet Imager and Characterizer) is a series of upgrades to Keck II/NIRSPEC (McLean et al., 1998; Martin et al., 2018; López et al., 2020) and Keck II AO to enable high-resolution ($R\sim 35,000$) diffraction-limited spectroscopy for Keck. Primarily intended for spectroscopic follow-up of directly imaged companions, the use of single-mode fibers offers a factor of $\sim$100 reduction in sky background, as well as long-term blaze and line-spread function (LSF) stability. Rasmussen et al. (2021) demonstrated that even minor variations in the LSF or blaze function can significantly impact detection of close-in planets through cross-correlation techniques. Previously, Finnerty et al. (2021) found that cross-correlation techniques are much more severely impacted by minor systematic errors compared with modest increases in Gaussian noise. These findings suggest that the improvements in LSF/blaze stability with KPIC may offset the increase in Poisson noise due to lower throughput in high-resolution cross-correlation applications compared with seeing-limited slit-fed spectrographs. We note that the KPIC phase I throughput was comparable to AO-fed NIRSPEC (Finnerty et al., 2022) and that KPIC phase II significantly exceeds the throughput performance of NIRSPAO (Echeverri et al., 2022).

Our observations are summarized in Table 2. For all observations, we used an ABBA pattern nodding between KPIC science fibers 1 and 2. These fibers had the highest throughputs of the four science fibers in KPIC phase I (tested at the start of each night), and the nodding allows us to effectively subtract sky emission features and thermal emission from warm front-end optics on the Keck AO bench. We obtained 90-s exposures at each nod position, giving our observations a time resolution of approximately four minutes after read time/overheads and coadding frames. We chose a time resolution under five minutes based on simulations showing that the orbital motion of WASP-33b would begin to significantly broaden the final cross-correlation peak for longer timescales (known as Doppler smearing). Figure 1 presents a diagram of the system with the observed orbital phases shaded by observation date.

The 2021 November 21 observations began during secondary eclipse, providing a stellar spectrum with no planet contribution. Observations continued post-eclipse for just under three hours, covering roughly noon to mid-morning longitudes. Conditions were good, with throughput from the top of Earth’s atmosphere through detection peaking at consistently around 2.5%.

We attempted to obtain additional observations on 22 November covering late-evening to afternoon longitudes. However, high clouds caused significant extinction and poor AO performance, leading to top-of-atmosphere throughput consistently below 1%. We therefore omitted these observations from subsequent cross-correlation and retrieval analysis.

We began by fitting the echellogram trace location and width for each fiber and order from observations of a bright calibrator star (typically the telluric standard for the wavelength calibrator, discussed below). Because the KPIC fiber-extraction unit (FEU) does not move relative to the spectrometer backend during science operations, this only needs to be done once per night. We then AB subtracted the science frames (taken in an ABBA sequence) to remove sky emission and thermal background. Spectral extraction was performed by summing over slices in the detector-Y direction (which is well aligned with the cross-dispersion direction) centered on the best-fit trace center with a width five times the standard deviation of the best-fit Gaussian.

Each fiber was individually wavelength-calibrated by comparing to a radial velocity standard following the procedure described in Section 3.2 of Wang et al. (2021). In brief, we observed a late-type giant star (HIP 95771 in the case of the WASP-33 observations) at the start of each night and fit the observed spectrum with a combined telluric spectrum generated using PSG (Villanueva et al., 2018) and stellar model from the PHOENIX library (Husser et al., 2013). Using a late-type star greatly increases the number of lines used for calibration compared with using only the sparse $K$-band telluric absorption features, enabling a more robust calibration.

We next resampled all the extracted spectra from science fiber 1 onto the wavelength grid obtained for science fiber 2 using a cubic spline. While each fiber is stable in detector location and wavelength over a single night, manufacturing and alignment differences can lead to small variations between fibers in both fringing and LSF that we wish to average. We arrange the extracted spectra for each order into a times series, giving $n_{\mathrm{order}}$ arrays each of shape $n_{\mathrm{nods}}\times n_{\mathrm{channels}}$. These arrays, along with the Julian date/time at the midpoint of each nod pair and the corresponding barycentric radial velocity to WASP-33 from Keck, are the inputs to our retrieval code. We omit the first three orders (blueward of 2.1 $\mu\rm m$) due to strong telluric contamination, and the following two orders due to wavelength calibration inaccuracies, leaving four NIRSPEC orders spanning $\sim 2.2-2.5\ \mu$m, with significant gaps. Additional post-processing performed as part of the retrieval is discussed in Section 3.1.1.

As our retrieval framework by necessity includes a forward-modeling capability, we can easily test our retrievals on a simulated system with known inputs. This provides a check on systematic biases in the retrieval and ensures that the reported uncertainties are reasonably accurate. These checks also allow us to determine if any parameters are inherently poorly constrained by the available data. All simulated data sets used the same orbital phase sampling as the observations (2021 November 21).

As our forward modeling does not include fringing or airmass variation, we do not include the PCA step in the test retrievals. The planet spectrum is the only time-varying component in the simulated spectrum. Even with our negative-injection approach, the planet spectrum is distorted by the PCA in the absence of stronger fringe/airmass-induced variations. This limits our ability to identify biases introduced from the PCA and may result in a mild overestimation of the retrieved errors in the simulations, as the negative-injection PCA approach is expected to increase preference for the true planet model. Ongoing developments in forward modeling of the KPIC fringing should reduce future reliance on PCA for defringing.

Our first test case was a star-only (no planet) simulated data set. This provides a non-detection baseline for comparison. Figure 10 shows the full corner plot (Foreman-Mackey, 2016) in Appendix A. Despite using the same stopping criteria as other retrievals, the star-only simulation shows much flatter and noisier posteriors and reached the stopping criteria in far fewer iterations. Additionally, the scale parameter shows a weak preference towards smaller values, suggesting a preference for the absence of the planet, as expected. The retrieved pressure-temperature profile is nearly isothermal, and the maximum-likelihood planet spectrum is nearly flat, as the fitting minimizes the strength of the planet features.

Our second test case was based on the retrieved parameters for WASP-33b listed in the “Input” column of Table 4. The signal-to-noise ratio of the simulations was chosen to be similar to that of the observations. This provides both a general check for systematics or unconstrained parameters in our retrieval framework as well as a specific injection/recovery test for the retrieved dayside atmosphere. Retrieved values and uncertainties are listed in the “Retrieved” column of Table 4, and the full corner plot is presented in Figure 11.

All parameters except the OH mass fraction are retrieved to within 1$\sigma$ of the input values, with error bars similar to those retrieved from the observations. The PT profile retrieved from the simulations is slightly hotter than the input profile but is generally within the 1$\sigma$ bound and does not impact the retrieved abundances. The SiO mass fraction, H mass fraction, and H dissociation fraction are all poorly-constrained in both the observations and the simulated retrieval, suggesting our observations have limited sensitivity to these parameters.

While the OH posterior shows a peak in the observations, this is not the case in the simulated retrieval. The ¹³CO isotopologue ratio and the scale factor are also better-constrained in the observations than the simulations, though not to the same degree as the OH abundance. While this may be a result of the negative-injection PCA approach used for the observations improving sensitivity to relatively small changes in the spectrum, the poor sensitivity in simulated retrievals suggests both apparent constraints may be artifacts.

We also ran retrievals using emcee (Foreman-Mackey et al., 2013). This provides a check that any biases are not due to the particulars of our nested sampling approach, such as the stopping criteria or the number of live points. The resulting median values and uncertainties from emcee were consistent with the nested sampling results. This suggests that systematics are more likely to be related to limitations in the radiative transfer calculation, log-likelihood function, or inherent limitations in the data, rather than the statistical framework used.

We use the best-fit retrieval parameters to compute $\Delta v_{\mathrm{sys}}-\Delta K_{p}$ diagrams, both for the overall model and for each individual species. These are plotted in Figure 6, with the retrieved $\Delta v_{\mathrm{sys}}$ and $\Delta K_{p}$ values indicated in dashed green. This is the approach that previous high-resolution studies have used to make molecular detections (e.g. Brogi et al., 2012; Buzard et al., 2020, and many others). We note that attempts to quantify strengths for this type of detection have been fraught (e.g. Buzard et al., 2020; Finnerty et al., 2021; Buzard et al., 2021), but these plots can still provide useful qualitative insight.

Specifically, quantitative estimates of the detection strength from $K_{p}-v_{sys}$ plots typically assume values far from the planet peak are normally-distributed and use the standard deviation far from the planet as an estimate of the noise. In the case of Figure 6, this is complicated by the clear presence of systematic variations as a function of $K_{p}$ due to the detrending suppressing the planet at small $K_{p}$. We therefore estimate the noise in Figure 6 by first masking $K_{p}<50$ $\rm km\ s^{-1}$ to exclude the region where the PCA is having a significant impact, then additionally masking values above/below the 10th/90th percentile to remove the planet feature. We then calculate the standard deviation and use this to make the detection strength map. These thresholds are arbitrary, and the remaining points after masking still do not appear to be normally distributed. We emphasize that Figure 6 is intended for a qualitative assessment of the detection of different species. For a quantitative assessment, Figure 4 is better for estimating the strength of the overall detection, and Table 6 presents Bayes factors for the detection of each molecules.

The strong dependence of the detrending signal on $K_{p}$ allows it to be removed by subtracting each column of Figure 6 by its median. This results in a detection strength of 0 at $K_{p}=0$, as expected when the planet is totally removed by detrending. For the all molecule case, the planet is detected at an SNR$\sim$12, which is more consistent with expectations based on Figure 4 than the scale of Figure 6. However, this removal of the detrending systematics is not statistically rigorous. We emphasize that Figure 4 offers the most robust estimate of the overall detection strength.

The planet is clearly detected in the $K_{p}-v_{sys}$ diagram plotted in the top panel of Figure 6. In the 2D plots of cross-correlation versus time in Figure 4, planet velocity track is clearly visible, and the retrieved velocities are in good agreement with the nominal values. Additional validation of the detection is provided by the disappearance of the planet feature during secondary eclipse in Figure 4.

The cross-correlation detection is dominated by ${}^{12}\rm CO$ based on comparing the first and second panels of Figure 6. OH also produces a strong signal at the expected location in the $K_{p}-v_{sys}$ space (fourth panel of Figure 6), while the H₂O template produces a somewhat weaker signal (fifth panel of Figure 6). The ${}^{13}\rm CO$ template has a weak feature coincident with the expected planet velocity, but this feature does not dominate and we do not consider this to be an independent detection of ${}^{13}\rm CO$. The SiO template does not produce any features in the $K_{p}-v_{sys}$ space and is not shown.

Retrieved atmospheric parameters are listed in Table 5. The full corner plot (Foreman-Mackey, 2016) for the retrieval is presented in Figure 12. The retrieved dayside PT profile, emission contribution function, and median spectrum are plotted in Figure 7.

The dayside atmosphere shows a clear inversion beginning at approximately 0.1 bar, consistent with previous results (e.g. Nugroho et al., 2021; van Sluijs et al., 2022). The $K$-band emission contribution function (Figure 7, top right panel) indicates that the bulk of the observed emission arises near this pressure, though CO continues to contribute to the line cores from altitudes up to $\sim 3\ \mu$bar. Pressures greater than $\sim 1$ bar do not contribute to the observed emission, and the retrieved PT profile is therefore unlikely to be accurate at these pressures. We note that the high-resolution spectroscopy is sensitive to the line strength relative to the continuum, which is determined by the shape and contrast of the PT profile, the atmospheric metallicity, the scaling parameter, and possible impacts from H^-. This can allow “cold” PT profiles to produce similar output spectra to those of “hotter” PT profiles. While temperature will also change the relative strengths of different lines, this is a comparatively small effect over a limited bandpass. We therefore caution against relying on the exact temperatures of the retrieved PT profiles.

We obtain log mass fractions of $-4.1^{+0.7}_{-0.9}$ for H₂O, $-1.1^{+0.4}_{-0.6}$ for CO, and $-2.1^{+0.5}_{-1.1}$ for OH. This indicates that H₂O is almost entirely dissociated on the dayside of WASP-33b, consistent with previous findings from Nugroho et al. (2021). SiO is effectively unconstrained due to its relatively low $K$-band opacity, but peaks at a log mass fraction of $-3.6$, which would be consistent with the neutral Si detection reported in Cont et al. (2022a). We also obtain a weak constraint on the ¹³CO/¹²CO ratio, peaking around $10^{-1.6}$ but with a significant tail giving a median value of $10^{-1.7}$. The abundance parameters show a strong positive covariance, which enables better constraints on the C/O ratio than would be expected from the marginalized uncertainties.

Finally, we note that while the scaling parameter peaks near unity, there is a long tail towards higher values. The scaling parameter controls the strength of lines relative to the host star continuum, which can also be influenced by changes in the PT profile. Figure 12 shows weak degeneracies between the PT profile parameters and the scaling parameter, with larger values of the scaling parameter corresponding to PT profiles which produce weaker lines. This suggests the tail of the scaling parameter is a result of under-constraining the PT profile.

There are a number of systematics or physical effects that could also change the value of the scaling parameter. Our test retrievals suggest this may be a result of our data processing (see Figure 11). Physical factors include underestimates of the planet radius or temperature, H^- opacity, thermal expansion of the dayside atmosphere, extended emission from an outflow, or stellar pulsations. The retrieved posterior peaks at a scale factor $\sim 1$, suggesting that these effects are either relatively minor or largely cancel.

To assess the strength of our detections, we compute the Bayes factor comparing a flat/no planet model and the best-fit models for all molecules together and for each molecule independently. These values are presented in Table 6. The best-fit planet model is strongly preferred ($\rm\Delta BIC>100$) , with CO dominating the detection. The evidence for H₂O and OH is substantial but weaker than CO, consistent with expectations based on the weaker features in the $K_{p}-v_{sys}$ space in Figure 6.

Consistent with both previous results (Nugroho et al., 2021; van Sluijs et al., 2022) and global circulation models of UHJs (Wardenier et al., 2021; Komacek et al., 2022; Beltz et al., 2022), our retrievals confirm our assumption of a thermal inversion on the dayside. We find that constant offsets in the PT profile appear to have a limited impact on the final spectrum, suggesting our observations are not particularly sensitive to absolute temperature, with the curvature of the PT profile having a much more substantial impact. The inclusion of the scaling parameter may compound this, as it provides an alternative way to change line depths in the planet model. Figure 7 shows that our $K$-band observations are most sensitive to a relatively narrow range of pressures and therefore have limited ability to constrain the overall PT profile. Additional observations at other bands may improve the vertical constraints.

The velocity parameters prefer a slight blueshift compared to the nominal values in Table 1. While this could be a result of morning-to-evening winds, the shift is only 1-2$\sigma$ and is well within the model-dependent variations in velocity parameters reported in van Sluijs et al. (2022). Additional observations providing continuous coverage of a larger phase range would provide better constraints on dayside winds.

We adopted the $-9.2$ $\rm km\ s^{-1}$ systemic velocity from Gontcharov (2006). More recent observations have found $v_{sys}$ in the range of -3 $\rm km\ s^{-1}$ (Collier Cameron et al., 2010) to 0 $\rm km\ s^{-1}$ (Nugroho et al., 2021). van Sluijs et al. (2022) found systemic velocities from -15 $\rm km\ s^{-1}$ to +6 $\rm km\ s^{-1}$ , depending on the orbital phase of the observations and the type of model used for cross-correlation. A more positive systemic velocity would imply larger blueshifts in the retrieved planet parameters relative to the expected Keplerian motion in the host star reference frame.

The retrieved abundances yield a gas-phase $C/O=0.8^{+0.1}_{-0.2}$. The derived posteriors for C/H, O/H, and C/O are shown in Figure 8. However, these retrievals are missing several potential oxygen reservoirs in ultra-hot Jupiter atmospheres, whose presence in significant quantities would lower the atmospheric C/O ratio. Metal oxides including TiO (Cont et al., 2021; Serindag et al., 2021) and AlO (von Essen et al., 2019) have previously been detected in WASP-33b, neutral atomic oxygen has been detected in the similar UHJ KELT-9b (Borsa et al., 2021), and Si has been detected in low-pressure ($\lesssim 1$ mbar) optical observations of WASP-33b (Cont et al., 2022a), suggesting the possible presence of SiO deeper in the atmosphere. Our limited wavelength coverage precludes detection of any of these species, though a significant reduction in C/O as a result would require either significant neutral oxygen or an atmosphere substantially enriched in refractory elements compared with volatiles.

The atmospheric metallicity of WASP-33b is more difficult to infer from these observations than the C/O ratio. The retrieved abundances suggest a moderately metal-enriched atmosphere, but the lack of continuum information in high-resolution observations leads to strong correlations between species abundances, making absolute abundance estimates difficult. Previously, Cont et al. (2022a) found a solar-metallicity model best explained the observed Si signal, though they used a widely-spaced model grid. In subsequent analysis Cont et al. (2022b) estimated a metallicity of +1.5 dex based on retrievals of several atomic species from optical data, including Ti, Fe, and V, consistent with the retrieved metallicity from the KPIC observations.

In order to better understand the range of possible C/O ratios and metallicities, we compare our retrieved abundances with grids of equilibrium chemistry models from the poor_mans_nonequ_chem package of petitRADTRANS (Mollière et al., 2017, 2019). These comparisons are plotted in Figure 9. In our parameterization scheme, the CO abundance is determined primarily by the atmospheric metallicity, while the H₂O/OH abundance is significantly impacted by the C/O ratio. As the equilibrium grid does not include OH, we assume the equilibrium H₂O abundance as the prediction for the dayside OH abundance.

Figure 9 shows the retrieved abundances are most consistent with equilibrium models having $\rm[Fe/H]=1.2$ and C/O = 0.8. These values are consistent with Thorngren & Fortney (2019) modeled upper limits on metallicity for a planet the mass of WASP-33b as well as the retrieved metallicity from Cont et al. (2022b) based on optical observations. However, varying the $\rm[Fe/H]$ in the chemical model suggests that our metallicity uncertainty is approximately 1 dex, and we therefore cannot confidently exclude stellar or even sub-stellar atmospheric metallicities for WASP-33b. C/O appears to be better constrained, with a margin $\pm 0.1$ consistent with the posterior shown in Figure 8. Note that the strong covariance between species abundances suggests that this C/O ratio will hold even if the metallicity varies significantly (see Figure 8).

The chemical equilibrium models suggest significant SiO abundances at the pressures probed by our observations. At lower pressures, the SiO abundance drops due to thermal dissociation, consistent with the optical detection of Si by Cont et al. (2022a). While our $K$-band observations are not sensitive to SiO, $L$-band observations of the SiO bands from 4–4.2 $\mu$m could measure the SiO abundance. Assuming SiO is the dominant undetected oxygen reservoir at $\sim 1-0.01$ bar pressures, this would significantly improve estimates of the overall oxygen abundance.

The ${}^{13}\rm CO/^{12}CO$ posterior peaks at $\sim 10^{-1.6}$, with a tail towards lower values bringing the median to $10^{-1.7}$. While the posterior is somewhat poorly-constrained, the fact that we see even a weak preference is an argument in favor of high $\rm{}^{12}CO$ abundances. In simulated observations, mass-mixing ratios as high as $\rm\log CO=-1.5$ do not always lead to a preference in the isotopologue ratio posterior. The peak of this weak preference is consistent with estimated CO isotopologue ratios from protoplanetary disks (Woods & Willacy, 2009) and is consistent with the ${}^{13}\rm CO/^{12}CO\sim 10^{-1}-10^{-1.6}$ value measured for WASP-77Ab in Line et al. (2021). Constraints on isotopologue ratios as well as abundances may provide a complementary path to probing formation history and physical conditions in exoplanetary atmospheres (Mollière & Snellen, 2019; Zhang et al., 2021).

The high-C/O, potentially high-metallicity atmosphere suggested by our retrievals is consistent with the Pelletier et al. (2021) observations of the hot Jupiter $\tau$ Boo b. These observations conflict with predictions from pebble accretion models that atmospheric metallicity and C/O ratios should be inversely correlated, as secondary pebble accretion is expected to be dominated by oxygen-rich grains (Espinoza et al., 2017; Madhusudhan et al., 2017; Cridland et al., 2019; Khorshid et al., 2021). High C/O, high-metallicity atmospheres require formation and/or substantial accretion in an environment that is enriched in both carbon and solids. Such conditions are consistent with models that incorporate pebble drift (Booth et al., 2017), formation near CO or CO₂ snow lines (Öberg et al., 2011), or with pollution by C-rich grains interior to the H₂O snow line (Chachan et al., 2023).

As discussed in Section 5.3, one possible explanation for this discrepancy is that our $K$-band high-resolution retrievals are simply missing a significant source of oxygen. Figure 9 suggests relatively high SiO abundances for WASP-33b, which we cannot measure with our $K$-band observations. However, comparison with the equilibrium chemistry models suggests including SiO would not substantially decrease the C/O ratio under equilibrium conditions without substantial enrichment of Si relative to other metals (see Figure 9, third panel) Alternatively, additional oxygen could be present in the atomic form, as (Borsa et al., 2021) reported from transit observations of KELT-9b. In the future, Keck/HISPEC will offer simultaneous yJHK coverage at R$\sim 10^{5}$ (Mawet et al., 2019) with substantially improved throughput compared with KPIC, enabling characterization of a broader range of chemical species and reducing the likelihood of missing significant oxygen or carbon reservoirs.

Missing oxygen could also be the result of errors in linelists or molecular opacity calculations rather than chemistry. A poorly-matched H₂O or OH opacity in the retrieval setup could lead to a bias towards lower abundances. However, both this work and Pelletier et al. (2021) retrieved high-C/O atmospheres using different H₂O opacity tables and linelists, indicating that such an error would have to be systematic across multiple independent calculations, while simultaneously not significantly impacting the same calculations for CO. This suggests that the problem would have to be related to a common assumption specific to the H₂O high-temperature linelists. Computing high-resolution, high-temperature H₂O opacities for the $K$-band alone required nearly one year of CPU time using exocross (Yurchenko et al., 2018), limiting our ability to compute and compare multiple opacity tables.

Alternatively, superstellar metallicity and C/O may arise either from accretion near evaporation fronts close to snow lines (Öberg et al., 2011) or from carbon-rich grain accretion between the carbon soot line and the H₂O snow line (Chachan et al., 2023). Pebble accretion studies have found grain accretion within the snow line should boost oxygen abundances rather than carbon (e.g. Madhusudhan et al., 2017; Espinoza et al., 2017; Cridland et al., 2019), but models incorporating pebble drift can produce giant planet atmospheres with both high metallicity and high C/O. In particular, Booth et al. (2017) finds that the C/O $\sim 1$ and high carbon abundance we observe in WASP-33b is a specific indication that a planet accreted high-metallicity gas near the CO₂ snow line, which is located at $\sim 10$ AU for an A-type primary (Öberg et al., 2011). The present observed location and polar orbit of WASP-33b may then be explained through the eccentric Kozai-Lidov Effect (Naoz et al., 2011) after the protoplanetary disk dissipated. Future incorporation of refractory species abundances, such as Si and Fe, would allow the formation location to be specifically identified, as the refractory-to-volatile ratio is a diagnostic of the formation and accretion history (Lothringer et al., 2021). The presence of these species in the gas phase in the atmospheres of UHJs such as WASP-33b presents a unique opportunity to obtain a clear understanding of the evolutionary history of an exoplanet population.