A Precise Carbon and Oxygen Abundance Determination in a Hot Jupiter Atmosphere


 The origins of gas giant planets orbiting close to their host stars (``hot Jupiters'') remain a mystery despite more than a quarter-century of study (Fortney et al. 2021). The atmospheric compositions of these planets are highly sought after to provide insight to their formation location in protoplanetary disks, how they migrated to be so close to their host stars, and the relative role of solid versus gas accretion during their assembly (Madhusudhan 2019). However, simultaneous, bounded constraints on both carbon and oxygen abundances, which are key for understanding giant planet formation (Oeberg et al. 2011, Mordasini et al. 2016, Madhusudhan et al. 2017,Cridland et al. 2016), have been elusive (Kreidberg et al. 2014,Wakeford et al. 2018,Pelletier et al. 2021). Here, we report precise abundance measurements of both water and carbon monoxide in a hot Jupiter atmosphere via ground-based, high resolution spectroscopy. From these constraints on the primary carbon- and oxygen-bearing molecules, paired with upper limits on other minor volatile elemental carriers, we are able to derive the atmospheric elemental metal enrichment (metallicity) and the carbon-to-oxygen ratio (C/O). The inferred atmospheric metallicity is slightly sub-stellar (-0.48$+0.15/-0.13) and the C/O is consistent with stellar (0.59 ±0.08). The former is suggestive of a metal-depleted atmosphere relative to expectations based on extrapolation from the solar system, indicative of a greater partitioning of metals within the core vs the atmosphere. The C/O constraint rules out gas-dominated accretion followed by disk free migration. Taken together in the context of past inferences, these results point to a diversity of planetary atmospheric compositions in addition to the observed diversity of planetary system architectures.

quence on 14 December, 2020 with the Immersion GRating INfrared Spectrometer (IGRINS 11 ) 48 at the Gemini-South (GS) Observatory located on Cerro Pachón, Chile. Owing to the broad 49 wavelength range (1.43 -2.42 µm over 54 spectral orders), high spectral resolution (R∼45,000), 50 and sensitivity (SNR∼180-270/resolution-element), IGRINS on GS is particularly sensitive to the 51 molecular lines from multiple carbon, nitrogen, oxygen, and sulfur bearing species (see Methods).  Removal of the telluric contamination also removes any continuum level information in the 63 planet-to-star flux ratio 14 . In order to extract meaningful information from data processed in this 64 way 15,16 , we must first cross-correlate (CC) the data with model templates. Using a set of represen-65 tative thermal emission models that include the dominant absorbers expected at these temperatures 66 and over the IGRINS wavelength range (primarily, H 2 O and CO), we cross-correlate as a func-67 tion of velocity against each processed spectrum. Provided the model is an adequate template, 68 the CC function (CCF) for each spectrum reaches its maximum at the planetary velocity (a sum 69 of the system velocity and orbital velocity) at that specific orbital phase, and hence, trace out a 70 CC trail in velocity 17 . The CC trail is clearly visible in Fig. 1b, corresponding to the appropri-71 ate planetary velocity components, demonstrating that we are detecting the planetary atmosphere 72 as the planet orbits the star. We further leverage the circular orbital geometry, which predicts 73 the phase-dependent line-of-sight velocity/Doppler shift, to determine the total atmospheric signal 74 detection by summing over the CCF at each phase 15,16 (see Methods), for each pair of the plan-75 etary orbital velocity (K p ) and system velocity (∆V sys , relative to the reported literature value of 76 +1.685±0.0004 km/s, 10 ). Fig. 1c shows the total atmospheric thermal emission cross-correlation 77 signal-to-noise, peaking at an S/N=12.8 very near (offset by ∼-7 km s −1 ∆V sys , see Methods) the 78 anticipated pair of velocities, clearly indicating a strong detection of atmospheric thermal emission. 79 The next step is to identify the specific trace molecular species (the bulk atmosphere is pre- 80 dominately H 2 /He) present in the spectrum and retrieve their absolute abundances. To do this 81 we perform an automated procedure via a Bayesian inference (retrieval) scheme 14 (see Meth-82 2 ods). This approach simultaneously optimizes the volume mixing ratios for each trace species 83 (log 10 (n i ), i=H 2 O, CO, CH 4 , H 2 S, NH 3 , and HCN as well as a CO isotopic abundance parameter, 84 log 10 ( 13 C 16 O/ 12 C 16 O) where 12 C 16 O is the main isotopologue-see Methods), the vertical temper-85 ature structure, the planetary orbital and system velocities, and nuisance parameters to account for 86 uncertainties in the reported transit timing and possible signal stretching due the PCA analysis (see 87 Methods). This method accounts for all of the degeneracies that arise amongst the the multiple 88 overlapping molecular lines and absorption strength with atmospheric temperature gradient and 89 permits for absolute molecular abundance determinations 14,18 . 90 We achieve bounded constraints for log(n H2O ) and log(n CO ) (Fig. 2a) 93 1c). We are also able to retrieve a bounded constraint on the CO isotopic abundance ratio (see 94 Methods for a discussion). To quantify the significance of each species, we perform a nested model 95 comparison via a Bayes factor analysis 20 whereby each individual gas is removed and the retrieval 96 is re-run to compute 21 and compare the model evidences. The resulting log-Bayes factors amongst 97 the nested models are 298.2 for H 2 O (strongly preferred), 55.6 for CO (strongly preferred), -1.8

98
(not preferred, see Methods) for the combination of CH 4 +H 2 S+NH 3 +HCN, and 3.0 (moderately 99 preferred) for the 13 C isotope. This quantification of individual species detection is different than 100 the classic 22,23 frequentist gas detection methods (see Methods) that compare the on-peak CCF 101 (nominal K p , ∆V sys ) to an off-peak region. We prefer the former as it fully incorporates the 102 change in likelihood and prior volume that arises from the removal of a gas, whereas the latter 103 is only relevant to a specific model choice out of an ensemble of possible models. These data 104 also prefer a monotonically decreasing relatively cool temperature profile (Fig. 2d), suggestive of 105 an atmosphere that either has a fairly efficient day-to-night atmospheric circulation (the predicted 106 day-side temperature for poor circulation planets should follow the hotter inverted predicted profile 107 in Fig. 2d) or lacks high altitude UV/optical absorbers (e.g., metal hydrides and oxides), possibly 108 indicative of night-side condensation (cold-trapping) of refractory species 24 .

109
The intrinsic elemental abundances in a planetary atmosphere are illuminating quantities be-110 cause they are diagnostic of both atmospheric chemical processes (e.g., droplet condensation and 111 sedimentation) and formation conditions. Furthermore, C and O account for ∼70% 25 of the total 112 "metals" (e.g., any species heavier than H, He) in a typical solar-like composition gas, and are 113 hence, good tracers for the metal enrichment of an atmosphere. Since H 2 O and CO are the domi-114 nant C and O bearing molecules in this atmosphere (with relatively low abundance upper limits on 115 the other major C and O bearing molecules-see Methods), are expected to be largely unperturbed 116 by disequilibrium chemistry mechanisms at these temperatures 19 , and are expected to be homoge-117 neous with altitude over the pressures probed by typical observations 19 (Fig. 2c), we can convert 118 them directly into the elemental oxygen (n O = n CO +n H 2 O ) and carbon (n C = n CO ) abundances.

119
It is customary to normalize the elemental abundances relative to hydrogen (n i /n H ), relative to 120 that in the sun ([X/H]:=log 10 ((n X /n H )/(n X /n H ) sun )) 25 to facilitate comparisons with other as-121 trophysical bodies in a common abundance reference frame. We find the elemental abundances in and a ratio of carbon to oxygen, C/O=0.59±0.08 (the Solar value is 0.55) (Fig. 2b). We also 124 3 retrieve a sub-terrestrial 12 C/ 13 C abundance ratio (10. formation than for our own solar system giants.

141
It is with the sheer numbers of exoplanets that we can quantitatively test specific formation- 164 Improvement in our understanding of how atmospheres came to be and how they evolve will con-165 4 tinue as we push towards higher precision abundance measurements of more targets from both 166 ground-and space-based platforms, ultimately paving the way for understanding our own Solar 167 system's formation history in the galactic context.  b) Cross correlation coefficient as a function of orbital phase/spectrum and planet velocity using a model template from the Bayesian inference procedure described in the text. The white trail corresponds to higher cross-correlation values (hence, atmospheric signal) and is consistent with the predicted velocity trail given the planetary orbital velocity and system velocity (light yellow dashed line). c) Atmospheric day-side thermal flux detection signal-to-noise (detection of absorption due to H 2 O and CO, see text) as a function of the planetary orbital velocity, K p , and the relative system velocity (∆V sys ) (see Text). The significance was computed by subtracting off the mean of the cross-correlation map and then dividing through by the standard deviation of a box far from the planet velocity pair. White dashed lines indicate the known 10 velocities (K p =192.06, ∆V sys =0 km s −1 ) and the "×" denotes the location of the peak signal (S/N=12.8)  The solar abundance value 25 is given as the black point and the solar abundance value accounting for oxygen sequestration due to potential condensate rain out 30 on the night side 24 is shown as the red point. The 1(39.3%)-2(86.4%)-and 3(98.9%)σ joint probability contours are indicated in both a and b and the numerical values above each histogram are the marginalized median and 68% confidence interval range. c) Vertical abundance profiles for the major species predicted with both equilibrium (dashed) and disequilibrium (vertical transport and photochemistry, solid) chemistry (see Methods). d) Retrieved vertical temperature structure (magenta, 68 and 95% confidence intervals) compared to 1D radiative-convective equilibrium models with the coldest resulting from efficient day-to-night heat transport, the hottest poor heat transport, and the middle, poor heat transport but with nightside condensation of refractory species (see Methods). Comparison of the IGRINS WASP-77Ab abundance constraints with the solar system planets (a, adapted from 31 ), with other exoplanets, and interior structure based envelope predictions (b). In (b), the gray points are from a uniform HST transmission spectrum water retrieval analysis by 28 . The solar system planets (black diamonds) follow a decreasing trend (dotted line) with increasing mass 7,28 . The light blue and green dots are the predicted envelope enrichments for the gas rich planet population based upon their mass and radius measurements 32 . The blue dots assume a coreless planet with all metals (e.g., C and O) uniformly mixed throughout the gas, whereas the green dots assume that 99% of the metals are in a solid planetary core (1% in the envelope).  The 1D atmosphere is parameterized with constant-with-altitude volume mixing ratios for 235 the aforementioned gases (6 gases with He+H 2 as the remainder with n He /n H 2 =0.176 ) and the 236 9 6-parameter temperature profile scheme described in 48 . For completeness we also include as a free 237 parameter the CO isotopic abundance ratio, 13 C 16 O/ 12 C 16 O (log 10 relative to the terrestrial value 238 of 1:89, more on this below). For HRCCS specific retrievals, we must also include the planet 239 Keplerian and system velocities (∆K p and ∆V sys relative to the literature 10 reported values of 240 192±11 and 1.6845±0.0004 km s −1 , respectively). Finally, we include as nuisance parameters 241 a stretching term to the planet flux to account for uncertainties in the reported planet/star radius 242 or data reduction induced stretching and a phase offset term to account for errors in the reported 243 ephemera. Extended data Table 1 summarizes the parameters and their uniform prior ranges.

244
The model planet spectrum is convolved with both a planetary rotational (vsin(i)=4.52 km/s) 245 and an instrumental broadening (assumed to be Gaussian) kernel followed by an interpolation ref. 14 except that we use the PCA instead of the airmass detrending method. To do so, we save 252 the N comp discarded eigenvectors from the SVD to reconstruct the telluric/systematic data matrix 253 followed by a multiplicative injection of the model (Doppler shifted matrix of 1+(F p /F star )(λ, φ)).

254
The PCA/SVD is then re-applied to the model injected data for each order. This matrix is then    O enrichments and carbon-to-oxygen ratio (main text Fig. 1b). The total C abundance is given (0.29 -0.54×Solar) and C/O=0.46±0.08. We use the "non-rainout" abundances in the main-text 319 discussion and Fig. 3.

320
The C and O abundances for WASP-77Ab are interpreted through the lens of the solar sys-321 tem abundance determinations, the representative exoplanet population abundances as measured 322 with low spectral resolution platforms (e.g., HST), and theoretical models in main text Fig. 3.

323
To compare to the solar system (main text Fig. 3b) we use the abundances given in Table 1  Ref. 32 provide predictive models for the maximum metal enrichment (based upon O) for the 342 exoplanet population given their measured mass and radius for a "core-less" planet, e.g., the metals 343 and gas are well mixed throughout the entire planet. This is clearly extreme as these values (blue 344 dots, main text Fig. 3b) vastly overshoot the measured Jupiter and Saturn envelope enrichment's, 345 suggesting a large fraction of metals must be sequestered into a solid core (on the order of 90%).

346
To match the retrieved depletion for WASP-77Ab, approximately 99% of the accreted metals must To test the robustness of the abundance and temperature profile constraints, we perform a battery 354 of tests that explore the impact of data processing and modeling assumptions on the retrieved H 2 O,

356
The first test is used to evaluate the influence of the TP-profile parameterization ("Tempera- Finally, we under go an independent Bayesian/retrieval analysis using an entirely indepen-386 dent tool/code (HyDRA-H 18,72 ), but also utilizing the log-likelihood mapping from 14 and PCA for 387 airmass detrending. A comparison of a subset of common parameters is shown in ED Fig. 6. In 388 this comparison, HyDRA-H retrieves for identically the same parameters as described in ED Table   389 1 with the following differences: the NH 3 , HCN, H 2 S, and CH 4 gas mixing ratios are not included decreases with decreasing pressure (e.g., no thermal inversion). 396 We thus conclude that the resulting constraints, and subsequent derivatives there-of, pre- . For these reasons, we include this ratio as a free parameter (ED Table   405 1).

406
Surprisingly (ED Fig. 3 ED Fig. 1: Summary of the data and PCA procedure. a) shows the median per-resolution element signal-to-noise for each order for the night (in red). The blue curve is the median SNR in both time and over an individual order. b) shows example raw data cubes (top row)-spectra vs. time/frame for representative two orders (25, 5). Stationary tellurics show up as vertical dark streaks. Wavelength dependent gradient is due to the echelle blaze throughput. The PCA/SVD method can remove these stationary features, leaving behind the planetary signal buried in the noise (bottom row). We use these "Post-PCA" frames for the subsequent cross-correlation/retrieval analysis (repeated for all 43 use orders  ED Fig. 6) The inset shows the molecular components of the maximum likelihood model spectrum.  Fig. 4: Classic cross-correlation analysis data products (about the maximum likelihood solution/spectrum from the retrieval analysis). The left column illustrates the gas detection's (all gases, H 2 O, CO, and other-NH 3 +H 2 S+HCN+CH 4 ) in the standard K p -V sys plane, with a slice in V sys along the literature reported K p at the bottom. The detection maps are constructed by subtracting the mean total CC, then dividing by an "off peak" (a boxed region in the lower left corner of each panel) CC standard deviation. Using this method, only H 2 O is strongly detected, with a hint of CO present at the expected velocities. The right column reproduces analogous products using the log-likelihood formalism 14 (∆logL relative to the minimum), resulting in a stronger presence of CO. We emphasize, that while such maps may be instructive for "detecting" species or "atmosphere", they do not marginalize over all of the degeneracy, nor do they maximize the information content in the data. This is why in our analysis, we focus on the the results arising from the more comprehensive log-likelihood/retrieval formalism.   Fig. 7: Carbon isotopic abundance analysis. The top row of histograms compares the constraints from a nominal simplified retrieval model applied to the the true data (red) and an the reverse-injected data re-injected with 13 C isotope removed model (black). The upper limit on the simulated data and bounded constraint arising from the true dataset suggests that there is indeed isotopic information in these IGRINS data. The bottom panel compares the retrieved 12-C to 13-C ratio (red) to common solar system bodies (blue) and various reference values (black dashed lines). WASP-77Ab sits anomalously low (enhanced 13-C) compared to most solar system objects.