Atmospheric Metallicity and C/O of HD 189733 b from High-resolution Spectroscopy

We present high-resolution K-band emission spectra of the quintessential hot Jupiter HD 189733 b from the Keck Planet Imager and Characterizer. Using a Bayesian retrieval framework, we fit the dayside pressure–temperature profile, orbital kinematics, mass-mixing ratios of H2O, CO, CH4, NH3, HCN, and H2S, and the 13CO/12CO ratio. We measure mass fractions of logH2O=−2.0−0.4+0.4 and logCO=−2.2−0.5+0.5 , and place upper limits on the remaining species. Notably, we find logCH4 < −4.5 at 99% confidence, despite its anticipated presence at the equilibrium temperature of HD 189733 b assuming local thermal equilibrium. We make a tentative (∼3σ) detection of 13CO, and the retrieved posteriors suggest a 12C/13C ratio similar to or substantially less than the local interstellar value. The possible 13C enrichment would be consistent with accretion of fractionated material in ices or in the protoplanetary disk midplane. The retrieved abundances correspond to a substantially substellar atmospheric C/O = 0.3 ± 0.1, while the carbon and oxygen abundances are stellar to slightly superstellar, consistent with core-accretion models which predict an inverse correlation between C/O and metallicity. The specific combination of low C/O and high metallicity suggests significant accretion of solid material may have occurred late in the formation process of HD 189733 b.

First discovered in 2005 through both transit and radial velocity observations (Bouchy et al. 2005), the hot Jupiter HD 189733 b is a frequent target for characterization studies due to its bright K-type host star (K mag = 5.5, Cutri et al. 2003) and large transit depth (∼ 2%).Winn et al. (2006) measured the Rossiter-McLaughlin effect for the system, finding the planetary orbit to be well aligned with the stellar rotation axis (λ = −1.4• ± 1.1 • ).Deming et al. (2006) detected the secondary eclipse in Spitzer observations, eventually leading to the first full-disk temperature map for an exoplanet (Knutson et al. 2007) which found the hottest part of the planet is offset by 16 • ± 6 • east from the substellar point due to supersonic winds.
Early efforts to characterize the chemical composition of HD 189733 b were met with mixed results, a selection of which we summarize here.Some studies reported detection of H 2 O (Tinetti et al. 2007;Grillmair et al. 2008) or CO (Désert et al. 2009), while others reported a flat transmission spectrum (Grillmair et al. 2007;Pont et al. 2008;Sing et al. 2009).The flat spectra were found to be consistent with the presence of hazes (Pont et al. 2008;Sing et al. 2009), preventing the clear detection of individual molecular species from low-resolution transmission spectra.Subsequent observations suggest this haze extends from the optical into the near-infrared H band (Gibson et al. 2012).
High-resolution emission spectroscopy provides a potential pathway to perform atmospheric characterization in the presence of clouds or hazes.Since high-resolution spectroscopy probes line cores at lower pressures compared to transmission spectroscopy, emission originating above the cloud layer can be detected and characterized (Gandhi et al. 2020).de Kok et al. (2013) used highresolution transmission spectroscopy to detect CO emission and place an upper limit on the H 2 O line contrast.Brogi et al. (2016) used transmission spectroscopy to measure a rotation rate consistent with tidal locking and a blueshift suggesting a ∼2 km s −1 day-to-night wind.Water was eventually detected in high-resolution transmission spectroscopy from GIANO (Brogi et al. 2018) and has been repeatedly confirmed in both emission and transmission spectroscopy (Flowers et al. 2019;Cabot et al. 2019;Brogi & Line 2019;Boucher et al. 2021;Klein et al. 2023).CO has similarly been confirmed in multiple studies (Brogi & Line 2019;Cabot et al. 2019;Flowers et al. 2019), while Cabot et al. (2019) also reported tentative evidence for HCN.
Despite the repeated confirmation of the presence of both H 2 O and CO, none of these studies were able to obtain simultaneously bounded constraints on the abundances of both species.This has prevented robust mea-surement of the C/O ratio and bulk atmospheric metallicity, potentially key diagnostics to trace the formation and evolution of the HD 189733 system.Additionally, despite the anticipated presence of CH 4 under local thermal equilibrium at the equilibrium temperature of HD 189733 b, none of these studies made a clear CH 4 detection, leaving an important potential tracer of photochemistry effects unconstrained.
To address these gaps in our knowledge of an important benchmark hot Jupiter system, we observed HD 189733 b with Keck/KPIC high-resolution K-band spectroscopy.These observations are part of an ongoing survey to characterize a statistically significant number of hot Jupiter atmospheres with KPIC and constrain the underlying distribution of C/O and metallicity in this population.
Section 2 describes the observations, data reduction procedures, and atmospheric retrieval framework, with an emphasis on the differences from Finnerty et al. (2023a).Section 3 presents the results of the atmospheric retrievals and cross-correlation analysis.Section 4 discusses the retrieved atmospheric properties, comparing the measured abundances with expectations from chemical equilibrium and photochemical models.We summarize our results and draw conclusions in Section 5.

OBSERVATIONS AND DATA REDUCTION
2.1.Observations HD 189733 was observed with Keck/KPIC phase II (Delorme et al. 2021;Echeverri et al. 2022) on 2022 August 13 (UT) from 8:39 UT to 11:59 UT at an airmass of 1.0 to 1.36.KPIC provides a single-mode fiber feed into Keck/NIRSPEC (McLean et al. 1998;Martin et al. 2018;López et al. 2020), enabling stable diffractionlimited high-resolution spectroscopy.The observations began approximately 30 minutes after the expected end of secondary eclipse, covering orbital phases from 0.5285 to 0.5898, during which time we expect the star/planet relative radial velocity to shift from −27 km s −1 to −82 km s −1 .Observations were taken with an ABBA nodding pattern between KPIC science fibers 2 and 4 with 30 second exposures, giving a total of 256 science frames.
The KPIC real-time throughput calculator (included in the KPIC DRP) estimated the 95th-percentile top-ofatmosphere throughput was 3-3.4% at the start of the observations, decreasing to ∼ 2.5% in the second half of the sequence, roughly consistent with the typical performance measured for KPIC phase II under photometric conditions (Echeverri et al. 2022).Weather conditions were clear and stable throughout the observations.The throughput loss is likely a result of an increase in the PSF size due to seeing and AO performance changes as the airmass increased from 1.0 at the start of the observation sequence to 1.36 by the end.The median extracted per-pixel signal-to-noise ratio after coadding frames in a nodding sequence was approximately 180, ranging from ∼ 210 in the bluest order to ∼ 160 in the reddest order.

Data Reduction
The data extraction is identical to that described in Finnerty et al. (2023a), which differs slightly from the current version of the supported KPIC Data Reduction Pipeline (DRP) 1 .HIP 95771 (spectral type M0.5IIIb) 1 https://github.com/kpicteam/kpicpipeline/ was used for wavelength calibration at the start of the night.Of the nine observed NIRSPEC orders, three (orders 37-39) are heavily contaminated by telluric CO 2 lines, while two (orders 35 and 36) are almost entirely lacking in stellar or telluric lines, preventing an accurate wavelength calibration.We therefore use the four reddest orders (orders 31-34) spanning ∼2.2-2.5 µm (with significant gaps) for which we have accurate wavelength solutions and minimal telluric contamination.
In order to speed the log-likelihood calculation in the retrieval, we coadded the extracted spectra from each ABBA sequence after re-interpolating all observations onto the science fiber 2 wavelength solution.Each fiber has a slightly different line-spread function (LSF) and the fiber coupling can vary significantly between frames depending on the AO correction, leading to variability in the coadded LSF.The final time series consists of 64 spectra, with a maximum planetary radial velocity shift between consecutive spectra of approximately 900 m s −1 , roughly 10× smaller than the NIRSPEC resolution (R ∼ 35, 000, ∆v ∼ 8.5 km s −1 for the chosen instrument settings).
The data detrending process for one order is shown in Figure 1.The extracted spectral time-series for each order is first scaled to have a median of 1.0, then divided by the time-series median spectrum using a second-order polynomial to fit continuum variations.This removes most temporally fixed features in the spectrum, leaving time-varying fringing (Finnerty et al. 2022;Finnerty et al. 2023b, Xuan et al. submitted), airmass-driven telluric variations, and the time-varying planet signal.Points differing from the median by more than 6× the median absolute deviation (MAD) are then masked, as well as the first and last fifty wavelength bins in each order.
The time series for each order is then mean-subtracted and the first 4 components of the singular value decomposition (SVD) are projected out.To account for the resulting distortion of the planet signal, the dropped components are saved and added to the forward model during the log-likelihood calculation.A new SVD is then performed on the component-injected forward model in order to replicate the distortion of the planet signal from projecting these components out, similar to Line et al. (2021).In contrast with Finnerty et al. (2023a), we do not negatively inject the forward model prior to removal of the SVD components, which significantly speeds the log-likelihood calculation.This choice makes a stronger assumption that repeating the SVD on the componentinjected model accurately reproduces the distortion of the observed planet signal, but has proven successful in previous work (e.g.Line et al. 2021).After the de- composition, points varying by more than 4× the MAD are masked.Previously, Finnerty et al. (2023a) found several orders of KPIC observations continued to show a fringing pattern that dominated over the Gaussian noise even at this stage of the analysis, which was removed by taking a final median division.In these observations, we do not see any evidence of similar residual fringing above the noise level (for example in the third panel of Figure 1), and therefore we omit this step.

Atmospheric Retrieval
The atmospheric retrieval framework is similar to that described in Finnerty et al. (2023a), though we have implemented several changes to improve the time-series detrending.We briefly summarize this framework for completeness.The fit parameters, priors, and retrieved values are listed in Table 2.
We use the 4-parameter pressure-temperature (P − T ) profile model from Guillot (2010), implemented in petitRADTRANS.From initial tests using the Madhusudhan & Seager (2009) 6-parameter P −T model previously used in Finnerty et al. (2023a), we found that multiple parameters were poorly constrained and showed complex degeneracies.The Guillot (2010) parameterization is more constraining in its physical assumptions, making it better-suited to the limited wavelength range of these observations and enabling looser, physically motivated priors for the P − T parameters.Retrieved molecular abundances were generally consistent between both parameterizations.At all pressures, we require 100 K < T < 3000 K for physical consistency and to avoid hitting the bounds of our opacity tables.We also include a grey cloud deck and use the scattering mode of petitRADTRANS.
We fit for vertically-fixed abundances (not necessarily in chemical equilibrium) of H 2 O, 12 CO, CH 4 , NH 3 , HCN, H 2 S, and the 13 CO/ 12 CO ratio.For H 2 O, 12 CO, and 13 CO we used the high-temperature opacity tables described in Finnerty et al. (2023a).For CH 4 , we used the opacity table computed in Xuan et al. (2022) from the Hargreaves et al. (2020) linelist.For NH 3 (linelist: Yurchenko et al. 2011), HCN (linelist: Harris et al. 2006;Barber et al. 2014), and H 2 S (linelist: Rothman et al. 2013), we used the petitRADTRANS high-resolution opacity tables described in Mollière & Snellen (2019).Collision induced absorption (CIA) from H 2 −H 2 and H 2 −He are also included.We fit for the abundance of H 2 , though we do not include H 2 line opacity.The remaining mass of the atmosphere is assumed to be helium.As we do not include H 2 line opacity, the retrieved H 2 abundance reflects the best-fit mean-molecular weight via the impact on the CIA, rather than strengths of H 2 lines in the planetary spectrum.We use a PHOENIX model for the star with T eff = 5000 K, log g = 4.5, and [Fe/H] = 0.0, which has been broadened to v sin i = 1.8 km s −1 (Polanski et al. 2022).While Finnerty et al. (2023a) fixed the size of the Gaussian broadening kernel to a value found by varying the kernel size for a fixed planet model and maximizing the likelihood, in this work we set the kernel width to σ = 1.2 pixels, slightly smaller than the expected value for the instrument, and apply a rotational broadening kernel to the atmospheric forward model using the fast technique described in Carvalho & Johns-Krull (2023) which accounts for the wavelength-dependence of the rotational kernel.The size of this kernel is a free parameter in the fit, allowing us to determine the best-fit line width accounting for coadd-induced LSF variations and the change in projected planetary velocity within a single coadd.
The nested sampling was performed with dynesty (Speagle 2020) using 1200 live points, running to a ∆ log z = 0.01 stopping criteria.The retrieval took approximately four days with 16 Intel E5-2670 CPU cores, totaling about 60 days of CPU time.The significant increase in runtime compared with Finnerty et al. (2023a) is a result of including scattering in the radiative transfer calculation.

RESULTS
Table 2 presents the priors and results from the retrieval, including both the maximum-likelihood and median retrieved values.We discuss the velocity parameters in Section 3.1, the P − T profile and thermal properties in Section 3.2, and the chemical abundances in Section 3.3.The full corner plot is included in Appendix A.

Velocity and winds
The v sys − K p diagram for the maximum-likelihood planet model is shown in Figure 2. Note that selfdivision of the planet spectrum results in substantial structure far from the planet peak, which prevents an accurate estimation of detection significance from dividing by the off-peak standard deviation.This effect can also be seen in Figure 6 of Finnerty et al. (2023a), and is a general feature of detrending schemes which rely on the velocity shift of planet features that is then compounded by the explicit dependence of the Brogi & Line ( 2019) log-likelihood function on the forward-model variance.As an alternative to estimate the detection strength, we use Wilks' Theorem (Wilks 1938) with 17 free parameters to estimate the significance of the detection compared with a flat planet model to be ∼ 6σ.Alternatively, we note that the off-peak structure is a strong function of K p .Subtracting each column of the v sys − K p by its median, while not statistically robust, mostly removes this structure and allows a more reliable estimate of the off-peak variance.Using the K p < 0 values to estimate the standard deviation after this subtraction gives a detection strength of 7.8σ.More rigorously, the cross-correlation coefficient depends only indirectly on the model variance, and the variance of the crosscorrelation coefficient is therefore a weaker function of K p , though self-division of the model still has some impact at small K p .Computing the v sys − K p diagram using the cross-correlation coefficient and dividing by the standard deviation of the K p < 0 region gives a detection strength of 7.3σ near the retrieved maximumlikelihood velocity parameters.
The retrieved planet velocity is blueshifted by several km s −1 .Based on the orbital phase sampling of the observations, we expect that the previously reported dayto-night wind would appear as a redshift.The apparent blueshift is still consistent with values in the literature (e.g.Boucher et al. 2021;Klein et al. 2023), but may suggest a more complicated circulation pattern than a single day-night wind.Alternatively, a number of other factors could bias our retrieved velocity, which we discuss in the next section.
The broadening velocity reported in Table 2 is substantially larger than values previously reported for the planetary rotational velocity in the literature (e.g.3.4 +1.3 −2.1 km s −1 in Brogi et al. 2016) and is inconsistent with the expectation that HD 189733 b be tidally locked, which would correspond to a rotational velocity of 2.6 km s −1 .This is due to several simplifying assumptions in our handling of the KPIC LSF which are absorbed into a larger broadening velocity.Specifically, we assumed a smaller nominal LSF smaller than Finnerty Table 2. List of parameters, priors, and results for atmospheric retrievals.The error bars on the retrieved medians correspond to the 68%/1σ confidence interval.In addition to these priors, we required both a non-inverted P − T profile and that the atmospheric temperature stay below 3000 K at all pressure levels.The full corner plot is included in Appendix A.
et al. (2023a), allowing the retrieval to freely determine the best-fit line width.This helps account for LSF variations between coadded science spectra that result from the slight differences in LSFs between science fibers and the varying frame-to-frame coupling efficiency, and can also help account for the line smearing due to the planet's motion within a coadd.While this prevents us from actually measuring the physical rotational speed of HD 189733 b, the tightly constrained posterior suggests these data have sufficient sensitivity to the planetary line shapes that future analyses will be able to measure v sin i for hot Jupiters.

Thermal properties
Figure 3 shows the retrieved P − T profile, emission contribution function, maximum-likelihood spectrum, and opacities of the major species.While several of the P −T parameters are poorly constrained in the 1D marginalized posteriors, repeated draws from the posterior yield a consistent set of profiles, suggesting the P −T profile is over-parameterized for the available data, but is well-constrained in the physically significant space.
The retrieved P − T profile is slightly cooler than previous theoretical models for HD 189733 b (e.g.Tsai et al. 2021).This may contribute to the preference for larger values of the overall multiplicative scale factor applied to the model planet spectrum, which results in an overall flux level consistent with the expected equilibrium temperature for HD 189733 b (Figure 3, middle panel).Alternatively, lower cloud pressures could be counteracted by increasing the scale factor in order to maintain the effective strength of the planet lines relative to the stellar continuum.Additionally, Finnerty et al. (2023a) found in simulations that this retrieval framework prefers larger values of the scale factor, likely due to artifacts in the data processing.
Due to this uncertainty, we ran an additional retrieval with a Uniform(0,10) prior on the scale factor.The scale factor was again close to the upper bound, and the retrieved P − T profile was even colder.However, the retrieved abundances were consistent at the 1σ level and the flux level of the maximum-likelihood planet model was consistent to within ∼ 10% for both retrievals.Initial analysis of another hot Jupiter observed with KPIC The detrending process suppresses planet features at small Kp, producing the feature centered on Kp = 0, where the planet model is a flat line after detrending.Compared with this null case, the maximum-likelihood model is preferred by ∆ log L = 39, equivalent to ∼ 6σ using Wilks' Theorem (Wilks 1938) with 17 free parameters.Note that the significant off-peak structure prevents accurate estimation of detection significance from division by the standard deviation far from the planet feature.(Finnerty et al. in prep.) is also showing a similar preference for a higher-than-expected scale factor that is countered by a cooler-than-expected P − T profile to match the expected continuum level.This strongly suggests our free retrieval framework has limited sensitivity to absolute temperature from the K-band data alone, but that this uncertainty does not significantly impact the retrieved atmospheric composition.Finally, we also note that residual continuum slopes or offsets in either the data or the model could also result in a spurious preference for a larger scale factor.In this case, the scale factor parameter is attempting to replicate the continuum offset, while the P − T profile is shifted in order to preserve the correct line strength relative to the scaled continuum.

Chemical composition
Of the included species, we obtain bounded constraints only for CO and H 2 O.Of the remaining species, we obtain upper limits on CH 4 , NH 3 , and HCN, though we note that the marginalized HCN posterior shows weak preference for log HCN ∼ −4.8.These species all have opacities comparable to that of the 2.3 µm CO bandhead in at least one of the observed NIRSPEC orders, suggesting these species would have been detected if present in significant abundances.While H 2 S is effectively unconstrained, the opacity table used has a gap from 2.25−2.35µm, which may preclude detection.The retrieved posterior also prefers a very high 13 CO/ 12 CO ratio, peaking at ∼ 10 −0.8 , though with a substantial tail to lower values.We include detection maps for H 2 O, CO, CH 4 , and 13 CO in Appendix B, as well as a discussion of the challenges associated with accurately estimating single-molecule detection strength.

Comparison to previous results
The extensive existing literature on HD 189733 b provides a basis for comparison with our retrieval results.To more easily facilitate comparison to other results, we note that the median retrieved values values reported in Table 2 are equivalent to K p = 160 ± 8 km s −1 , log CO VMR = −3.3±0.5, and log H 2 O VMR = −2.9±0.4.The remainder of this paper will continue to use massmixing ratios for molecular abundances.Klein et al. (2023) analyzed high-resolution infrared transmission spectra from CFHT/SPIROU (Donati et al. 2020) to perform an atmospheric retrieval of the H 2 O abundance and P −T profile using the same Guillot (2010) parameterization with a grey cloud deck.Klein et al. (2023) report K p , v sys , and T eq in good agreement with the values in Table 2, and log H 2 O = −2.95+0.75  −0.53 , slightly smaller than the value in Table 2, while reporting a larger value for the cloud-top pressure log P cloud = 0.479 +1.02 −1.06 .We note that the corner plots for both retrievals indicate a degeneracy between the cloud pres- The observed NIRSPEC orders are shaded in grey.In addition to the maximum-likelihood and median P − T profiles, the top left also includes the corresponding cloud top pressures as dashed horizontal lines and the P − T profiles from 100 draws from the retrieved posterior.While several parameters are poorly constrained in the corner plots, the actual P − T -profiles follow a tight distribution.The emission contribution function shows the emission mostly arises near 100 mbar, just above the cloud deck, with contribution from higher altitudes in the CO line cores.The dashed blue line plotted with the maximum-likelihood spectrum shows the flux from a 1200 K blackbody, which as expected is comparable to the retrieved planet flux.The slope of the retrieved spectrum differs from a blackbody as a result of CO and especially H2O absorption features at the red end of the K-band.Our retrievals provide only upper limits on NH3 and CH4 abundances despite these species' substantial opacity across the entire observed band, suggesting that the limits indicate a real absence of these species.
sure and molecular abundances, with greater pressures corresponding to lower abundances, potentially explaining the minor discrepancy in these parameters.
The data analyzed in Klein et al. (2023) were originally described in Boucher et al. (2021), who reported a somewhat lower H 2 O abundance, higher cloud pressure, and lower temperature.Their corner plots also show a degeneracy between the abundance, cloud-top pressure, and temperature which could explain the discrepancy.The consistency of this degeneracy, regardless of analysis framework, suggests a need for analyses covering a broader wavelength range that can better constrain cloud properties.
Constraints on the atmospheric CO abundance of HD 189733 b have been more elusive.de Kok et al. (2013) reported the first CO detection for HD 189733 b, also based on K-band emission spectroscopy, but did not report an abundance due to uncertainties concerning the impact of hazes.Brogi & Line (2019)

Winds and circulation
Compared with the nominal values in Table 1, the retrieved ∆K p and ∆v sys yield a blueshift of 8 km s −1 at the start of the observation sequence, increasing to 10 km s −1 by the end.While wind speeds as high as 8 km s −1 have been reported for HD 189733 b from transmission observations of sodium lines (Wyttenbach et al. 2015), other observations in the near-infrared have preferred a windspeed of 2 km s −1 (Brogi et al. 2016), while still other infrared transmission analyses have reported K p and v sys values compatible with those reported in Table 2 (Boucher et al. 2021;Klein et al. 2023).We expect that a day-to-night wind would appear as a redshift in post-eclipse emission observations, rather than a blueshift.Similarly, we would expect the hotter dayside dominating the overall emission signal of the rotating planet to lead to a slight additional redshift.Improved orbital phase coverage may be required to directly con-strain the wind speeds and 3D circulation on HD 189733 b from high-resolution spectroscopy alone.
Factors other than planetary winds could be contributing to the observed blueshift.The velocity shift is too large to be easily explained by an error in the orbital phase, particularly given the precision of the transit time for HD 189733 b.Similarly, we expect our wavelength solution to be reliable to ∼ 100 m s −1 , substantially better than the observed shift.Rasmussen et al. (2023) demonstrated that the shift in planetary lines within a science exposure can produce a blueshift of several km s −1 , though we expect this effect to be minimal given the < 1 km s −1 velocity shift between successive spectra in the time series.However, we did not attempt to account for the wavelength-dependence or non-Gaussian wings of the NIRSPEC line-spread function (López et al. 2020), potentially leading to a model mismatch in the line shapes that could bias the retrieved velocities.A laser frequency comb compatible with KPIC is scheduled for deployment in late 2023, which, when completed, should enable the more robust calibration of the line-spread function required to interpret these types of velocity offsets as physical wind speeds or 3D effects.We note that a similar blueshift was also reported in Finnerty et al. (2023a) for WASP-33 b, which may indicate a systematic effect in the retrieval pipeline.
While we included a rotational broadening kernel in our fit, our approach to the instrumental broadening makes this effectively a nuisance parameter rather than a reliable estimate of the planetary v sin i.As discussed above, we ignore both the wavelength-dependence and non-Gaussian shape of the line-spread function, as well as the impact of line smearing from the change in planet velocity within an exposure.To account for this, we chose a slightly smaller than expected Gaussian kernel for the instrument profile, and allowed a larger rotational kernel to make up the difference based on the actual observed line widths after coadding fibers.Future improvements to our retrieval pipeline will include a wavelength-dependent instrument profile similar to that used by Wang et al. (2021) to measure v sin i for brown dwarf companions and directly-imaged planets observed with KPIC.Combined with the use of the laser frequency comb to constrain relative changes in the LSF across the NIRSPEC detector, this should enable robust measurements of v sin i for hot Jupiters in the future with KPIC.

CH 4 and photochemistry
Similar to Finnerty et al. (2023a), we compare our retrieved abundances with the results of the equilibrium chemistry calculator easyCHEM (Mollière et al. 2017, Lei Figure 4. Retrieved molecular abundances (solid lines, 1σ shaded) compared with the equilibrium abundances estimated from easychem (dashed lines).The equilibrium abundances are consistent with C/O = 0.25 and [M/H] = +0.05,both slightly below the values obtained from the gas-phaseonly retrieval.The equilibrium model predicts H2S and CH4 abundances substantially greater than the retrieved upper limits.These species may be photochemically depleted.
et al. in prep.) for the retrieved P − T profile.Figure 4 plots the equilibrium vertical abundance profiles with the constant vertical abundances from our retrievals.The retrieved CO and H 2 O abundances are a good match to the equilibrium predictions for C/O = 0.25 and [M/H] = 0.05, suggesting the gas-phase retrieval may be slightly overestimating the atmospheric metallicity.Both H 2 O and CO abundances are roughly constant with height in equilibrium, validating the constant-withaltitude assumption in our retrieval and indicating that the potential biases resulting from high-altitude H 2 O depletion in ultra-hot Jupiters (Brogi et al. 2023) are not a factor for HD 189733 b.
The cumulative distribution function (CDF) of the retrieved CH 4 posterior indicates log CH 4 < −4.6 with 99% confidence.However, Figure 4 indicates that in chemical equilibrium for the retrieved P − T profile, the CH 4 mass fraction should be > 10 −4 over nearly the entire pressure range probed by the KPIC observations in order to match the retrieved H 2 O and CO abundances.We caution that the retrieved P − T profile in Figure 3 is somewhat colder than expected in the upper atmosphere, which may lead to an overprediction of CH 4 when assuming chemical equilibrium.While previous attempts at CH 4 detection in hot Jupiters have faced challenges with linelist accuracy, we use an opacity table based on the recent Hargreaves et al. (2020) linelist which has been validated on brown dwarfs of similar ef-fective temperature (Tannock et al. 2022).This opacity table has also previously been used to successfully measure the CH 4 abundance in a brown dwarf companion from KPIC observations at a mixing ratio comparable to that expected in HD 189733 b (Xuan et al. 2022).This suggests that the retrieved upper limit on CH 4 is due to a real absence of the expected CH 4 , rather than a linelist issue.
Depletion of CH 4 could be due to photochemistry.We assessed this possibility using the VULCAN chemical kinetics code (Tsai et al. 2021).We used the default P −T /K zz profile for HD 189733 b included with VULCAN and described in Tsai et al. (2021), which is slightly hotter than our retrieved profile in the upper atmosphere, the SNCHO photochemistry network, and the C/H and O/H set to the retrieved medians.All other settings were unchanged from the VULCAN defaults.VULCAN predicts a roughly constant CH 4 mass fraction of ∼ 10 −5 from at pressures between approximately 1 and 10 −3 bar.Over the same pressure range, the NH 3 mass fraction is predicted to be ∼ 10 −4 .The abundance of both species rapidly drops above 1 mbar, which may bias the abundances obtained from a free retrieval assuming constant abundance with pressure towards lower values, similar to the impact of water depletion in the upper atmospheres of ultra-hot Jupiters discussed in Brogi et al. (2023).
From the retrieved posterior, we obtain 99% upper limits log CH 4 < −4.6 and log NH 3 < −4.5, suggesting the current KPIC observations are not quite sensitive enough to make a definitive (non)detection at the VULCAN-predicted abundances.Both CH 4 and NH 3 have significant opacity throughout the NIR.Incorporating additional data such as high-resolution L-band observations from KPIC phase III may improve sensitivity to these species and enable tests of photochemical models in the future.Additionally, incorporating fluxcalibrated broad-band medium-resolution observations from JWST would significantly improve constraints on absolute abundances compared with high-resolution observations alone, and also enable constraints on species which lack features easy to observe from the ground, such as CO 2 4.4. 13CO enrichment While the CO isotopologue ratio is not wellconstrained in the retrieval, comparing the loglikelihood of the maximum-likelihood model with and without 13 CO indicates the presence of 13 CO is favored at ∼ 3σ.The marginalized posterior shows a clear preference for high levels of 13 CO enrichment compared with both the solar system value of 12 C/ 13 C ∼ 89 and local interstellar value of 12 C/ 13 C ∼ 68 (Woods & Willacy 2009;Milam et al. 2005), but with a long tail that is still consistent with the interstellar value at 1σ.Previously, Finnerty et al. (2023a) and Line et al. (2021) reported indications of 13 CO enrichment in hot Jupiter atmospheres compared with the local ISM.These estimates are substantially less than the 12 CO/ 13 CO ∼ 10 +53 −6 median retrieved for HD 189733 b, and were instead broadly consistent with the level of isotopic fractionation expected in the midplanes of protoplanetary disks (Woods & Willacy 2009).
Carbon isotopologue constraints have also been obtained for widely separated companions.Zhang et al. (2021a) found a roughly solar value for the 12 CO/ 13 CO ratio in a young field brown dwarf and suggested a formation process with minimal preferential accretion of 13 CO-enriched material such as direct gravitational collapse.In contrast, 13 CO enrichment has been reported for a young accreting super-Jupiter (Zhang et al. 2021b), possibly suggesting preferential accretion of 13 C via fractionated ices.The enrichment reported in Zhang et al. (2021b) is similar to the values reported in Line et al. (2021) for WASP-77Ab and Finnerty et al. (2023a) for WASP-33b, still less than the preferred value of ∼ 10 for HD 189733 b.
The potential prevalence of 13 CO enrichment in hot Jupiters observed to-date suggests a common process in the formation or evolution of these planets.Accretion of fractionated ices appears to be consistent with the values reported by Line et al. (2021) and Finnerty et al. (2023a) (Woods & Willacy 2009), but may not be sufficient to explain the 12 CO/ 13 CO ∼ 10 preference seen in our HD 189733 b retrieval.Additional observational and modeling work is required to constrain the 12 CO/ 13 CO ratios of host stars and understand the general prevalence of this 13 CO enrichment, and the potential extreme case of HD 189733 b in particular.

C/O and metallicity
The posteriors for C/O, C/H, and O/H derived from the retrieval are plotted in Figure 5.As in previous high-resolution studies (e.g.Xuan et al. 2022;Finnerty et al. 2023a), the ratio of species is better constrained than the absolute abundances, giving C/O = 0.3 ± 0.1, log C/H = −3.1±0.5, and log O/H = −2.5±0.4.Precise abundances for HD 189733 A were previously reported in Polanski et al. (2022), enabling a robust comparison between the composition of the star and the planetary atmosphere.The CDFs of the distributions in Figure 5 give a sub-stellar C/O at 99.8% confidence, super-stellar C/H at 84% confidence, and super-stellar O/H at 99.4% confidence.Comparison to equilibrium models suggests the gas-only retrieval may be slightly overestimating the total atmospheric metallicity, but sub-solar metallicities are still disfavored.These results indicate the atmosphere of HD 189733 b is compositionally distinct from the host star.
The combination of substellar-C/O and super-stellar metallicity was suggested to be an indicator of significant post-formation ice accretion by Öberg et al. (2011).Madhusudhan et al. (2014) found that low-C/O, metalenriched atmospheres could be indicative of either core accretion followed by disk migration or formation at wide separations via gravitational instability followed by disk-free migration.Subsequently, Madhusudhan et al. (2017) found that such an atmospheric composition could also be the result of core erosion or late planetesimal accretion.
Core accretion should lead to atmospheric metal enrichment in giant planets (Thorngren et al. 2016), including a trend between planet mass and metal enrichment relative to the host star that is consistent with our reported values for HD 189733 b (Cridland et al. 2019).This process is also predicted to produce an inverse correlation between atmospheric metallicity and C/O (Espinoza et al. 2017;Cridland et al. 2019;Khorshid et al. 2022), though this trend may not hold in the presence of pebble drift (Booth et al. 2017).These models can predict the observed composition of HD 189733 b as a result of accreting a significant (∼ 10%) fraction of the total planetary mass as solids relatively late in the formation process, with radial pebble drift not substantially altering the overall C/O of this material.Such a scenario may also be consistent with the extremely high 13 CO enrichment suggested by the retrieval, which could be explained through late accretion of a large amount of highly fractionated ice.

SUMMARY AND CONCLUSIONS
We performed a Bayesian retrieval on high-resolution K-band observations of HD 189733 b from Keck/KPIC, successfully constraining the abundances of CO and H 2 O and placing upper limits on the atmospheric mass fractions of CH 4 , NH 3 , and HCN.The retrieved abundances yield an atmospheric C/O = 0.3 ± 0.1 for HD 189733 b, substantially less than the stellar C/O ratio, and a stellar-to-superstellar atmospheric metallicity.This composition could be explained by the accretion of a large amount of solid material late in the planet formation process, possibly as a result of disk migration, which could also produce a well-aligned orbit with respect to the stellar rotation.
We do not detect CH 4 in the atmosphere of HD 189733 b.While this is incompatible with the retrieved P − T profile in chemical equilibrium, the upper limits we obtained may be consistent with models which incorporate atmospheric photochemistry.Additional observations with a wider wavelength coverage will improve sensitivity to both CH 4 and NH 3 , enabling direct tests of photochemical models, haze formation mechanisms, and measurement of the N/O ratio -providing deeper insight into the formation and evolution of the HD 189733 planetary system.Mackey 2016), VULCAN (Tsai et al. 2021) petitRADTRANS (Mollière et al. 2019(Mollière et al. , 2020) ) APPENDIX A. CORNER PLOTS Figure 6 presents the full corner plot from the retrieval for completeness.We discuss the retrievals in Section 3, including poorly constrained and degenerate parameters.The units, priors, and retrieved quantity are listed in Table 2 B. MOLECULAR DETECTION MAPS Figure 7 presents v sys −K p diagrams for each detected species calculated using a leave-one-out approach.As shown in Figure 3, strong, broad absorption by CO and H 2 O cause significant changes to the overall shape of the planetary SED, particularly beyond ∼ 2.35µm.This presents a challenge when using single-molecule templates to assess detection strength, as the resulting inaccuracies in the continuum level can make it impossible for a single-molecule template to reproduce the observed relative fluxes between lines, particularly if other parameters are held fixed.
A slightly more reliable approach is to calculate v sys − K p diagram for the maximum-likelihood model as well as the v sys − K p diagram for a model omitting the molecule of interest.The difference between these two maps is the change in the log-likelihood that can be attributed to the presence of the omitted species.This approach will still be inaccurate in the case of e.g.H 2 O and CO whose omission has a significant impact on the continuum, but should provide a more accurate estimate for the detection strength of trace species.Taking the difference of the two maps should also largely remove the self-division artifact seen at K p = 0 in Figure 2, though the resulting map will show a bias towards higher values of K p .This bias arises from self-removal of the planet model during detrending and is increased by the explicit dependence of the Brogi & Line (2019) log-likelihood function on the model variance.The detrending process used to remove continuum and telluric features reduces the model variance as K p decreases and the self-division effect of the median division and the self-subtraction effect of the SVD remove more of the planet model.As a result, we expect that the apparent strength of noise features to increase with increasing K p when assuming a constant variance for the entire map, and that planetary features will show a bias towards higher K p in difference maps compared with computing log L directly.
These effects can be seen in Figure 7. H 2 O is detected at ∼ 3.8σ near the expected velocity, while CO is detected at ∼ 5σ but with a shift towards larger K p .CH 4 is not significantly detected, as expected, and a weak (∼ 3σ) feature is present in the 13 CO map near the maximum-likelihood velocities, but is not strong enough to constitute an independent detection, particularly given the presence of large-amplitude features at high K p .The H 2 O and CO detections are weaker than would be expected based off of the all-molecule model, likely due to the impact of these species on the planet continuum.However, 13 CO is significantly less abundant and will have a much smaller impact on the planetary continuum level.Consistent with this, the 13 CO feature in the v sys − K p map is similar in strength to the expectation from Wilks' Theorem in Section 4.4.
Consistent with previous high-resolutions studies of other hot Jupiters (Line et al. 2021;Finnerty et al. 2023a), 13 CO is not independently detected in the v sys − K p space, but the presence of 13 CO is favored by ∼ 3σ over a model without 13 CO.As discussed above, estimating the detection strength from the v sys − K p diagram poses a number of challenges that are particularly acute for trace species which do not produce a strong peak to begin with.The repeated tentative detection of 13 CO across multiple targets highlights the need for improved techniques to quantitatively estimate detection strengths for high-resolution observations.subtracting a map omitting each species from a map made using the maximum-likelihood model, in order to better account for the impact of other atmospheric species on the continuum.Significance was estimated by dividing each map by the standard deviation of the Kp < 0 region.Continuum impacts are still significant, leading to weaker-than-expected detections of H2O and CO and biasing the CO velocity.These effects are less important for lower-abundance species.As expected based on the retrieved posteriors, we do not detect CH4, and 13 CO shows a weak (∼ 3σ) feature near the expected planet velocity.This 13 CO feature does not constitute and independent detection, but the inclusion of 13 CO is favored over a 12 CO-only model at ∼ 3σ, consistent with other results for 13 CO in hot Jupiter atmospheres (Line et al. 2021;Finnerty et al. 2023a) 1. INTRODUCTION arXiv:2312.00141v1[astro-ph.EP] 30 Nov 2023

Figure 1 .
Figure 1.Data detrending process for KPIC HD 189733 observations.The top row shows the raw time series, with significant frame-to-frame flux variations, tellurics, and blaze function effects.The colorbar for the top frame is the counts in the extracted 1D spectrum.The second row shows the same data after scaling each exposure by its median and dividing each frame by the time-series median, with the colorbar indicating the median-normalized and continuum-subtracted flux level.The blaze and telluric signals are almost completely eliminated, but a significant time-varying fringe is clearly present.White pixels indicate data that have been masked.In the bottom panel, the fringe has been effectively removed by dropping the first 4 components of the SVD decomposition of the time series.The bottom colorbar shows the median-normalized, continuum-subtracted flux level after the decomposition.

Figure 2 .
Figure 2. vsys − Kp diagram for the maximum-likelihood model.The planet is clearly detected, though with a significant blueshift compared to the expectation indicated in dashed red.The detrending process suppresses planet features at small Kp, producing the feature centered on Kp = 0, where the planet model is a flat line after detrending.Compared with this null case, the maximum-likelihood model is preferred by ∆ log L = 39, equivalent to ∼ 6σ using Wilks' Theorem(Wilks 1938) with 17 free parameters.Note that the significant off-peak structure prevents accurate estimation of detection significance from division by the standard deviation far from the planet feature.

Figure 3 .
Figure 3. Retrieved P − T profile (top left), maximum-likelihood emission contribution function (top right), maximumlikelihood planet spectrum (middle), and opacities for H2O,CO, NH3, and CH4 (bottom).The observed NIRSPEC orders are shaded in grey.In addition to the maximum-likelihood and median P − T profiles, the top left also includes the corresponding cloud top pressures as dashed horizontal lines and the P − T profiles from 100 draws from the retrieved posterior.While several parameters are poorly constrained in the corner plots, the actual P − T -profiles follow a tight distribution.The emission contribution function shows the emission mostly arises near 100 mbar, just above the cloud deck, with contribution from higher altitudes in the CO line cores.The dashed blue line plotted with the maximum-likelihood spectrum shows the flux from a 1200 K blackbody, which as expected is comparable to the retrieved planet flux.The slope of the retrieved spectrum differs from a blackbody as a result of CO and especially H2O absorption features at the red end of the K-band.Our retrievals provide only upper limits on NH3 and CH4 abundances despite these species' substantial opacity across the entire observed band, suggesting that the limits indicate a real absence of these species.
re-analyzed the de Kok et al. (2013) observations of HD 189733 b to retrieve CO and H 2 O, but obtained only lower and upper limits, respectively.The analysis presented in Section 3 is the first simultaneous determination of CO, CH 4 , and H 2 O abundances in the atmosphere of HD 189733 b suitable for constraining the atmospheric C/O ratio.Several JWST programs have obtained spectroscopy of HD 189733 b in both near and mid infrared wavelengths (GO 1633, PI Deming; GO 2001, PI Min; GO 2021, PI Kilpatrick & Kataria), which should provide a second, independent estimate of the atmospheric metallicity and C/O ratio (including the impact of CO 2 ) to compare with our results.

Figure 5 .
Figure 5. Retrieved C/O, C/H, and O/H posteriors.The median value is indicated in solid red, with the upper and lower bounds of the 68% confidence interval in dotted red.Stellar values from Polanski et al. (2022) are shown in dashdot blue.The retrieved C/O ratio is substantially sub-stellar, while the carbon and oxygen abundances are both superstellar.

Figure 7 .
Figure7.vsys − Kp diagrams for H2O, CO, CH4, and 13 CO.Dashed red lines indicate the nominal values of the planetary velocity parameters, and the red dot indicates the maximum-likelihood values for the full model.These maps were made by subtracting a map omitting each species from a map made using the maximum-likelihood model, in order to better account for the impact of other atmospheric species on the continuum.Significance was estimated by dividing each map by the standard deviation of the Kp < 0 region.Continuum impacts are still significant, leading to weaker-than-expected detections of H2O and CO and biasing the CO velocity.These effects are less important for lower-abundance species.As expected based on the retrieved posteriors, we do not detect CH4, and 13 CO shows a weak (∼ 3σ) feature near the expected planet velocity.This 13 CO feature does not constitute and independent detection, but the inclusion of 13 CO is favored over a 12 CO-only model at ∼ 3σ, consistent with other results for 13 CO in hot Jupiter atmospheres(Line et al. 2021;Finnerty et al. 2023a)

Table 1 .
Stellar and planetary properties for the HD 189733 system.