# Detection and Bulk Properties of the HR 8799 Planets with High-resolution Spectroscopy

Jason J. Wang, Jean-Baptiste Ruffio, Evan Morris, Jacques-Robert Delorme, Nemanja Jovanovic, Jacklyn Pezzato, Daniel Echeverri, Luke Finnerty, Callie Hood, J. J. Zanazzi, et al.

► To cite this version:

Jason J. Wang, Jean-Baptiste Ruffio, Evan Morris, Jacques-Robert Delorme, Nemanja Jovanovic, et al. Detection and Bulk Properties of the HR 8799 Planets with High-resolution Spectroscopy. The Astronomical Journal, 2021, 162, 10.3847/1538-3881/ac1349. insu-03705275

HAL Id: insu-03705275

https://insu.hal.science/insu-03705275

Submitted on 27 Jun 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Detection and Bulk Properties of the HR 8799 Planets with High-resolution Spectroscopy

Jason J. Wang<sup>1,24</sup> (王勁飞<sup>1,24</sup>), Jean-Baptiste Ruffio<sup>1</sup> (Jean-Michel Ruffio<sup>1</sup>), Evan Morris<sup>2</sup>, Jacques-Robert Delorme<sup>1,3</sup> (Jean-Michel Delorme<sup>1,3</sup>), Nemanja Jovanovic<sup>1</sup> (Nemanja Jovanovic<sup>1</sup>), Jacklyn Pezzato<sup>1</sup>, Daniel Echeverri<sup>1</sup>, Luke Finnerty<sup>4</sup> (Luke Finnerty<sup>4</sup>), Callie Hood<sup>2</sup> (Callie Hood<sup>2</sup>), J. J. Zanazzi<sup>5</sup> (J. J. Zanazzi<sup>5</sup>), Marta L. Bryan<sup>6,24</sup> (Martola L. Bryan<sup>6,24</sup>), Charlotte Z. Bond<sup>7</sup>, Sylvain Cetre<sup>3</sup>, Emily C. Martin<sup>2</sup> (Emily C. Martin<sup>2</sup>), Dimitri Mawet<sup>1,8</sup> (Dimitri Mawet<sup>1,8</sup>), Andy Skemer<sup>2</sup> (Andy Skemer<sup>2</sup>), Ashley Baker<sup>1</sup>, Jerry W. Xuan<sup>1</sup> (Jerry W. Xuan<sup>1</sup>), J. Kent Wallace<sup>8</sup>, Ji Wang (王吉)<sup>9</sup>, Randall Bartos<sup>8</sup>, Geoffrey A. Blake<sup>10</sup> (Geoffrey A. Blake<sup>10</sup>), Andy Boden<sup>1</sup>, Cam Buzard<sup>11</sup> (Cam Buzard<sup>11</sup>), Benjamin Calvin<sup>1,4</sup>, Mark Chun<sup>12</sup> (Mark Chun<sup>12</sup>), Greg Doppmann<sup>3</sup> (Greg Doppmann<sup>3</sup>), Trent J. Dupuy<sup>7</sup> (Trent J. Dupuy<sup>7</sup>), Gaspard Duchêne<sup>6,13</sup> (Gaspard Duchêne<sup>6,13</sup>), Y. Katherina Feng<sup>2</sup>, Michael P. Fitzgerald<sup>4</sup> (Michael P. Fitzgerald<sup>4</sup>), Jonathan Fortney<sup>2</sup> (Jonathan Fortney<sup>2</sup>), Richard S. Freedman<sup>14,15</sup> (Richard S. Freedman<sup>14,15</sup>), Heather Knutson<sup>10</sup> (Heather Knutson<sup>10</sup>), Quinn Konopacky<sup>16</sup> (Quinn Konopacky<sup>16</sup>), Scott Lilley<sup>3</sup>, Michael C. Liu<sup>17</sup> (Michael C. Liu<sup>17</sup>), Ronald Lopez<sup>4</sup>, Roxana Lupu<sup>18</sup> (Roxana Lupu<sup>18</sup>), Mark S. Marley<sup>19</sup> (Mark S. Marley<sup>19</sup>), Tiffany Meshkat<sup>20</sup> (Tiffany Meshkat<sup>20</sup>), Brittany Miles<sup>2</sup> (Brittany Miles<sup>2</sup>), Maxwell Millar-Blanchaer<sup>21</sup> (Maxwell Millar-Blanchaer<sup>21</sup>), Sam Ragland<sup>3</sup>, Arpita Roy<sup>22</sup> (Arpita Roy<sup>22</sup>), Garreth Ruane<sup>1,8</sup> (Garreth Ruane<sup>1,8</sup>), Ben Sappey<sup>16</sup>, Tobias Schofield<sup>1</sup>, Lauren Weiss<sup>17</sup>, Edward Wetherell<sup>3</sup>, Peter Wizinowich<sup>3</sup> (Peter Wizinowich<sup>3</sup>), and Marie Ygouf<sup>23</sup>

<sup>1</sup>Department of Astronomy, California Institute of Technology, Pasadena, CA 91125, USA; jwang4@caltech.edu

<sup>2</sup>Department of Astronomy & Astrophysics, University of California, Santa Cruz, CA 95064, USA

<sup>3</sup>W.M. Keck Observatory, 65-1120 Mamalahaia Highway, Kamuela, HI, USA

<sup>4</sup>Department of Physics & Astronomy, 430 Portola Plaza, University of California, Los Angeles, CA 90095, USA

<sup>5</sup>Canadian Institute for Theoretical Astrophysics, University of Toronto, 60 St. George Street, Toronto, ON M5S 3H8, Canada

<sup>6</sup>Department of Astronomy, University of California at Berkeley, CA 94720, USA

<sup>7</sup>Institute for Astronomy, University of Edinburgh, Royal Observatory, Blackford Hill, Edinburgh, EH9 3HJ, UK

<sup>8</sup>Jet Propulsion Laboratory, California Institute of Technology, 4800 Oak Grove Drive, Pasadena, CA 91109, USA

<sup>9</sup>Department of Astronomy, The Ohio State University, 100 W 18th Avenue, Columbus, OH 43210, USA

<sup>10</sup>Division of Geological & Planetary Sciences, California Institute of Technology, Pasadena, CA 91125, USA

<sup>11</sup>Division of Chemistry and Chemical Engineering, California Institute of Technology, Pasadena, CA 91125, USA

<sup>12</sup>Institute for Astronomy, University of Hawaii at Manoa, 640 North A'ohoku Place, Hilo, HI 96720-2700, USA

<sup>13</sup>Université Grenoble-Alpes, CNRS Institut de Planétologie et d'Astrophysique (IPAG), F-38000 Grenoble, France

<sup>14</sup>SETI Institute, Mountain View, CA 94043, USA

<sup>15</sup>NASA Ames Research Centre, MS 245-3, Moffett Field, CA 94035, USA

<sup>16</sup>Center for Astrophysics and Space Sciences, University of California, San Diego, La Jolla, CA 92093, USA

<sup>17</sup>Institute for Astronomy, University of Hawaii, 2680 Woodlawn Drive, Honolulu, HI 96822, USA

<sup>18</sup>BAER Institute, NASA Research Park, Moffett Field, CA 94035, USA

<sup>19</sup>Department of Planetary Sciences and Lunar and Planetary Laboratory, University of Arizona, Tucson, AZ 85721, USA

<sup>20</sup>IPAC, California Institute of Technology, M/C 100-22, 1200 East California Boulevard, Pasadena, CA 91125, USA

<sup>21</sup>Department of Physics, University of California, Santa Barbara, Santa Barbara, CA, USA

<sup>22</sup>Space Telescope Science Institute, Baltimore, MD 21218, USA

<sup>23</sup>NASA Exoplanet Science Institute, IPAC, Pasadena, CA 91125, USA

Received 2021 April 20; revised 2021 June 18; accepted 2021 July 9; published 2021 September 17

## Abstract

Using the Keck Planet Imager and Characterizer, we obtained high-resolution (R ~ 35,000) K-band spectra of the four planets orbiting HR 8799. We clearly detected H₂O and CO in the atmospheres of HR 8799 c, d, and e, and tentatively detected a combination of CO and H₂O in b. These are the most challenging directly imaged exoplanets that have been observed at high spectral resolution to date when considering both their angular separations and flux ratios. We developed a forward-modeling framework that allows us to jointly fit the spectra of the planets and the diffracted starlight simultaneously in a likelihood-based approach and obtained posterior probabilities on their effective temperatures, surface gravities, radial velocities, and spins. We measured v sin(i) values of 10.1±2.8 km s⁻¹ for HR 8799 d and 15.0±2.3 km s⁻¹ for HR 8799 e, and placed an upper limit of <14 km s⁻¹ of HR 8799 c. Under two different assumptions of their obliquities, we found tentative evidence that rotation velocity is anticorrelated with companion mass, which could indicate that magnetic braking with a circumplanetary disk at early times is less efficient at spinning down lower-mass planets.

*Unified Astronomy Thesaurus concepts: Exoplanet atmospheres (487); Exoplanet formation (492); High angular resolution (2167); High resolution spectroscopy (2096)*

## 1. Introduction

In the past two decades, direct imaging has discovered several dozen substellar companions with masses from 1 to 70 M<sub>Jup</sub> with orbital separations from ~3 au out to ~1000 au

(for a review, see Bowler 2016). The occurrence rates of giant planets and brown dwarfs beyond ~10 au have begun to show that multiple formation channels are responsible for the current population of imaged substellar companions (Nielsen et al. 2019; Vigan et al. 2020). Between 10 and 100 au, Nielsen et al. (2019) showed that giant planets between 5 and 13 M<sub>Jup</sub> have a higher occurrence rate compared to their brown dwarf counterparts (13–80 M<sub>Jup</sub>) and preferentially occur around higher-mass stars, indicating the known exoplanet companions likely formed as the high-mass tail of planet formation through core accretion (Pollack et al. 1996), whereas brown dwarf

<figure><img src="image_3.png" /></figure>

<footnote n="24"><sup>24</sup> 51 Pegasi b Fellow.</footnote>

<footnote>Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.</footnote>
companions formed like binary stars through gravitational instability (Boss 1998, 2001). Vigan et al. (2020) found a similar dichotomy looking at the mass ratios between the companions and their host stars. They inferred that companions around lower-mass stars with mass ratios closer to unity formed like binary stars whereas more extreme mass ratio companions around more massive stars formed like planets.

Because the directly imaged companions are amenable to spectroscopy and astrometric monitoring, we can use population-level characteristics beyond detection to study this population and understand how they formed. Bowler et al. (2020) reinforced the finding of multiple formation channels by showing evidence that giant planets at wide separations (5–100 au) had an eccentricity distribution similar to that of close-in (<1 au) giant planets, whereas the brown dwarf eccentricity distribution resembled the stellar-binary population. Measurements of the spin of planetary-mass companions have pointed to magnetic braking quickly slowing down the spin rate of planets after formation (Bryan et al. 2020b). While there has not yet been a population-level study of atmospheric compositions, compositional studies of individual objects are able to contribute evidence to discerning their formation channels (Konopacky et al. 2013; Barman et al. 2015; Gravity Collaboration et al. 2020; Mollière et al. 2020; Wilcomb et al. 2020).

High-resolution spectroscopy of directly imaged companions allows us to characterize their orbits, spin, and compositions. The Doppler shift of molecular absorption lines in the planetary atmosphere allows us to measure the radial velocity (RV) of the planet, which is useful to break degeneracies between orbital inclination and eccentricity for companions with limited orbital coverage (Snellen et al. 2014; Schwarz et al. 2016; Ruffio et al. 2019). The rotational broadening of these absorption lines allows for a direct measurement of planetary spin (Snellen et al. 2014; Schwarz et al. 2016; Bryan et al. 2018). The detection of molecular signatures through cross-correlation methods takes advantage of the fact that the planet and star have different spectral features, enables the detection of trace molecular species, and allows for the inference of planetary composition (Konopacky et al. 2013; Barman et al. 2015; Brogi & Line 2019; Wilcomb et al. 2020).

However, up until now, slit spectrographs assisted by adaptive optics (AO) have been observationally limited to bright companions that are at relatively large angular separations from their host stars. For planetary-mass companions within 1" (only β Pic b and HR 8799 c have been detected at a spectral resolution $R > 10,000$ (Snellen et al. 2014; Wang et al. 2018a). The main difficulty is that the bright glare of the host star overwhelms the signal of the planet. The glare of the star can be mitigated by combining high-contrast imaging techniques with high-resolution spectroscopy (Snellen et al. 2015). The combination of these two techniques, which we term high-dispersion coronagraphy (HDC), drives the design of the Keck Planet Imager and Characterizer (KPIC; Mawet et al. 2017a; Jovanovic et al. 2019; Delorme et al. 2021). KPIC combines starlight suppression using the Keck AO system and single-mode fibers with the NIRSPEC high-resolution spectrograph to enable high-resolution spectroscopy of fainter and closer-in directly imaged planets. Similar instrument designs are also being pursued by Subaru/REACH (Jovanovic et al. 2017; Kotani et al. 2020) and VLT/HIRISE (Vigan et al. 2018; Otten et al. 2021).

In this paper, we present the first science demonstration of KPIC with observations of the four planets orbiting HR 8799.

HR 8799 is a notable planetary system, as it is the only known system with four directly imaged exoplanets (Marois et al. 2008, 2010). The four planets are either near or in mean-motion resonance, and dynamical modeling of their orbits has constrained their masses to be $7.2 \pm 0.7 M_{Jup}$ for the inner three planets and $5.8 \pm 0.5 M_{Jup}$ for planet b (Wang et al. 2018b, also see Goździewski & Migaszewski 2020). Since their discovery, the planets have been characterized extensively in the 1–5 µm range using broadband photometry and low- and medium-resolution spectroscopy (Bowler et al. 2010; Barman et al. 2011; Konopacky et al. 2013; Currie et al. 2014; Ingraham et al. 2014; Skemer et al. 2014; Bonnefoy et al. 2016; Zurlo et al. 2016; Greenbaum et al. 2018; Mollière et al. 2020). At medium resolution, the individual molecular lines begin to be resolved spectrally, allowing for detection of molecular signatures as well as the RV of the planet (Barman et al. 2011; Konopacky et al. 2013; Ruffio et al. 2019). Measurements of elemental abundances in these planetary atmospheres have shown deviations from the stellar abundances, which have been interpreted to mean that these planets formed via core accretion rather than gravitational instability (Konopacky et al. 2013; Barman et al. 2015; Lavie et al. 2017; Mollière et al. 2020).

Section 2 details the observations of the HR 8799 planets made with KPIC. Section 3 describes the initial data reduction steps. The detection of molecular features through cross-correlation as well as fitting atmospheric models of rotating planets directly to the data is discussed in Section 4. We obtained the first spin measurements for these exoplanets, and we put them, as well as our measured orbital velocities and bulk atmospheric properties, in context in Section 5. We summarize our work and discuss future avenues both in obtaining better data and utilizing better models to study these planets in Section 6.

## 2. Observations

### 2.1. Instrument Description

KPIC consists of a series of upgrades for the Keck II AO system and its two facility instruments: the NIRC2 infrared imager and the upgraded NIRSPEC infrared high-resolution spectrograph (McLean et al. 1998; Martin et al. 2018; López et al. 2020). As part of this project, instrument upgrades included a new vortex coronagraph for NIRC2 (Vargas Catalán et al. 2016; Serabyn et al. 2017) and an infrared pyramid wave-front sensor operating in the *H* band for the AO system (Bond et al. 2020). Particularly relevant for high-resolution spectroscopy of directly imaged exoplanets, a fiber injection unit (FIU) following the concept presented in Mawet et al. (2017a) was deployed in 2018. The FIU benefits from the NIRSPEC detector upgrade that allows KPIC to reach $R \sim 35,000$ (Martin et al. 2014, 2018; López et al. 2020). We point the reader to Delorme et al. (2021) for detailed information about the instrumentation.

Here, we provide a brief summary of the relevant instrumentation. The pyramid wave-front sensor drives the facility deformable mirrors in the AO system to compensate for atmospheric turbulence. In addition to imaging exoplanets using NIRC2, KPIC can also send the *K*-band light of the system to the FIU to spectroscopically characterize these planets. Located after the AO system, the FIU steers the light of a planet into one of five single-mode fibers represented by circles in Figure 1. These fibers are connected to NIRSPEC to spectrally disperse the light injected into the fibers (see right panel of Figure 1). Because the
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Schematic of KPIC observations of the HR 8799 planets. On the left are the locations of the fibers relative to a Keck/NIRC2 PSF of the star for an observing sequence on HR 8799 c. The red circle marks the location of the fiber centered on the planet, while the yellow circles mark the positions of the three other fibers that were in use. A fifth fiber is not shown in this image as it was blocked by the NIRSPEC slit and was not used. Note that the NIRC2 image was not taken during the observing sequence and is just to serve as an example to show the dynamic range of the stellar halo. On the right is a small portion (~360 × 370 pixels) of the NIRSPEC detector image in this configuration. A portion of two echelle orders is shown. In each order, four spectral traces are seen corresponding to the placement of the four fibers shown on the left image.</figcaption></figure>

<table><caption>**Table 1**
HR 8799 Planet Observations</caption><thead><tr><th>Date</th><th>Target</th><th>Integration Time (minutes)</th><th>Airmass</th><th>Throughput (%)</th></tr></thead><tbody><tr><td>2020-07-01</td><td>HR 8799 c</td><td>230</td><td>1.0–1.8</td><td>1.7–1.9</td></tr><tr><td>2020-07-02</td><td>HR 8799 d</td><td>230</td><td>1.0–1.7</td><td>1.9–2.4</td></tr><tr><td>2020-07-03</td><td>HR 8799 e</td><td>130</td><td>1.1–2.0</td><td>2.2–2.4</td></tr><tr><td>2020-09-29</td><td>HR 8799 b</td><td>160</td><td>1.0–1.1</td><td>1.8–2.4</td></tr></tbody></table>

star is bright, diffracted starlight (stellar “speckles”) leaks into all of the fibers, with the amount of starlight in each fiber depending on the angular distance between the fiber position on the sky and the star. In Mawet et al. (2017b), we demonstrated that the use of single-mode fibers provides the following advantages: a well-defined Gaussian-like line spread function (LSF) that is independent of the shape of the input wave-front aberrations, and, on average, has 3× suppression of the underlying stellar speckle flux at the location of the fiber. Note that the current KPIC FIU does not utilize a coronagraph to suppress diffracted starlight, so the utilization of AO with single-mode fibers is the main starlight suppression mechanism.

The fiber extraction unit (FEU) reimages the light from the single-mode fibers onto the NIRSPEC spectrograph slit, where the light is then dispersed inside the spectrograph. The resulting NIRSPEC data have up to five fibers illuminated with light in each echelle order (only four fibers were imaged onto the detector for the HR 8799 observations). Each fiber roughly subtends an angular diameter of ~50 mas. The projected separation between two consecutive fibers is ~800 mas for the four fibers shown. Because the slit sees the FEU and not the sky, the slit background is dominated primarily by the thermal emission from the warm room-temperature optics of the FEU, with a lesser contribution from the thermal instrument background of NIRSPEC. The sky background only appears in the fibers but is generally well below the thermal background of the FEU, with only weak OH line emission seen in long exposures.

## 2.2. HR 8799 Observations

We listed the epochs of observations for each planet in Table 1. Each planet was observed with a similar observing sequence. NIRSPEC was set up to use the “Thin” filter, a thin piece of clear PK50 glass that blocks light at wavelengths

longer than ~2.5 µm, along with a custom pupil stop designed for the KPIC FEU and the NIRSPEC 0"/679 × 1"/13 slit to maximize the signal to noise of the light from each fiber relative to the instrument background. The NIRSPEC echelle grating was set to 63°/0 and the cross disperser was set to 35°/76 to allow us to obtain nine spectral orders ranging from approximately 1.94 to 2.49 µm.

We started each observing sequence by placing the host star on each fiber and took at least one exposure (30–60 s) of the host star in all four fibers for data calibration purposes. We then designated one of the fibers as the primary science fiber, based on the fiber that had the best throughput in daytime testing. We used the tip/tilt mirror in the FIU to offset the star from the fiber bundle such that the planet of interest is placed on the primary science fiber. The science fiber received a combination of planet light plus diffracted residual starlight. The offset amplitude and direction were computed using the where-istheplanet<fnref n="25" /> orbit prediction tool (Wang et al. 2021) based on the dynamically stable and coplanar orbit fit from Wang et al. (2018b). From preliminary instrument characterization efforts, we estimated that the offset accuracy is 10 mas, which corresponds to a <10% loss in throughput due to fiber misalignment. Although the fibers in the bundle were fixed in a linear arrangement relative to one another, the AO field rotator (K-mirror) was used to rotate the field of view such that the star was coupled into as many of the other fibers as possible so that we could obtain simultaneous stellar spectra.

We then took 600 s exposures with NIRSPEC using the MCDS-16 detector readout mode in this configuration where one fiber had the light of the planet and at least two other fibers were transporting significant amounts of starlight. This integration time was chosen to be long enough so that read

<footnote n="25">²⁵ http://whereistheplanet.com/</footnote>
noise is negligible (thermal background noise ~10× larger), but short enough to tolerate any bad frames due to technical issues such as the AO loops opening. Every hour or so, we placed the host star on each science fiber and took a short exposure (30–60 s) for calibration. The open-shutter time obtained on each planet is listed in Table 1.

Using the short exposures on the star, we calculated the end-to-end throughput from the top of the atmosphere to the detector. We reported the end-to-end throughput measured between 2.29 and 2.34 µm in Table 1 as a metric that combines both instrument performance and atmospheric conditions during the observations. We note that this wavelength range is not the best-performing wavelength, as Delorme et al. (2021) showed that a throughput of over 3% has been achieved at shorter wavelengths. However, our data analysis in the following sections is focused on this wavelength range because it coincides with the CO bandhead in the K band, so this is the most relevant throughput metric. Overall, we see that the conditions between the four nights are pretty comparable, with 2020 July 1 having slightly poorer performance due to issues in fiber injection that were fixed after the night.

## 2.3. Calibration Data

In addition to the spectra of the star we obtained during the observing sequence, we also observed an M giant (HIP 81497) for wavelength solution calibration and a telluric standard star for the wavelength calibrator (ups Her) each night. We took five 1.5 s exposures of HIP 81497 and three 30 s exposures of ups Her on axis in each of the four fibers. After the night, we took thermal background frames looking at the FEU with no light injected into the fibers at each of the exposure times used to model the thermal background of the instrument.

# 3. Data Reduction

## 3.1. Raw Data Reduction

The process of going from the original detector images to extracted 1D spectra was the same for all data regardless of the object being observed. First, the images were background subtracted using the instrument thermal background frames taken when no light was being injected into the fibers during the daytime. The FEU typically was at a different temperature during the day so these thermal frames did not perfectly subtract the data and left some residual background, which we modeled during the extraction step. Thermal images taken during the observing sequence would have provided better background subtraction (i.e., nodding), but were not acquired as we had not developed an efficient way to nod the planet light on the detector.

For each night of observation, the trace of each of the four science fibers in each of the nine orders was determined by using the data on the telluric standard star ups Her. As the point-spread function (PSF) of monochromatic light coming from a single-mode fiber is nearly a 2D Gaussian, we just needed to measure the position and standard deviation of the PSF at each wavelength. To do that, we fit a 1D Gaussian to the cross section of the trace in each column of each order to determine the position and standard deviation of the Gaussian PSF at that column. For each of the nine orders, we recorded the position and standard deviation of the Gaussian PSF for each of the four fibers in each column of the detector (2048 in total). To mitigate measurement noise and biases from telluric

lines, we smoothed the measurements by fitting a cubic spline to the measured positions and standard deviations. These PSF standard deviations will also be used for estimating the LSF width in Sections 4.2 and 4.3.

For each exposure, we then extracted the 1D spectra of each of the four fibers in each column of each order. Due to the imperfect background subtraction, we used pixels that are at least 5 pixels away from the center of any fiber to estimate the residual background level in each column and subtracted the median of these pixels from every pixel in the column. Then, for each fiber, we used optimal extraction to measure the flux using a 1D Gaussian profile as the PSF with the positions and standard deviations fixed to the values measured on the telluric standard star (Horne 1986). The total integrated flux of the 1D Gaussian is then the flux we extracted for that fiber in that column. The uncertainty in the optimal extraction was used as the uncertainty in the flux measurement (Horne 1986).

## 3.2. Wavelength Calibration

We observed the M giant HIP 81497 in each fiber each night to determine the wavelength solution from the stellar and telluric spectral lines. The wavelength solution was modeled as a spline using six nodes per order (i.e., piecewise third-order polynomials). We modeled the data as the multiplication of a stellar spectrum, a telluric transmission term, and the transmission of the telescope and the instrument, complemented by an additive linear background term. The star HIP 81497, with an M2.5III spectral type, was chosen from the Gaia RV standard catalog (Soubiran et al. 2018). It was modeled with a PHOENIX stellar spectrum (Husser et al. 2013) assuming a temperature of 3600 K, surface gravity of log(g) = 1, solar metallicity, and a fixed known RV of −55.567 ± 0.0011 km s⁻¹. The telluric transmission of Earth’s atmosphere was modeled from a linearly interpolated grid of 25 ATRAN<fnref n="26" /> models (Lord 1992) over the water vapor overburden (500, 1000, 5000, 10,000, and 20,000 µm) and zenith angle (0°, 25°, 45°, 68°, and 89°). The water vapor overburden and the zenith angle were fit as nuisance parameters. The transmission of the telescope and instrument varies with wavelength primarily due to the efficiency of the spectrograph to disperse light at each wavelength (i.e., blaze function). To model the spectrally dependent transmission, we used a piecewise-linear function that divided each order into five pieces. The best-fit wavelength solution was found using the Nelder–Mead optimization implemented in the `scipy.opt-imize.minimize` routine that jointly fit for the wavelength solution, telluric parameters, and instrument and telescope transmission (Virtanen et al. 2020). The search was initialized around the minimum of a coarse grid search to avoid local minima. This technique was found to be precise to the 0.1 km s⁻¹ level (Morris et al. 2020), which will be sufficient for the following data reduction. Morris et al. (2020) also found that the wavelength solution was stable within a night to the same level of precision as long as optics inside the spectrograph were not moved.

# 4. Fitting High-resolution Spectra

## 4.1. Forward Modeling the Planet Spectra

In this work, we build up a full forward model of the spectrum we recorded, which consists of both planet plus

<footnote n="26">²⁶ https://atran.arc.nasa.gov/cgi-bin/atran/atran.cgi</footnote>
starlight and is similar to the approach from Ruffio et al. (2019) for fitting medium-resolution integral field spectroscopy data of imaged exoplanets. The signal obtained from the fiber placed on one of the planets can be deconstructed into the following components:

$$D_p(\lambda) = \alpha_p(\lambda)T(\lambda)P_{\text{LSF}}(\lambda) + \alpha_s(\lambda)T(\lambda)S_{\text{LSF}}(\lambda) + n(\lambda). \quad (1)$$

Here, $D_p$ is the signal from the planet fiber, $T$ is the transmission of the optical system (atmosphere, telescope, and instrument) excluding fiber coupling efficiency, $P_{\text{LSF}}$ is the spectrum from the planet after it has been convolved by the instrumental LSF, $\alpha_p$ is a scaling term for the planet brightness accounting for fiber coupling efficiency, $S_{\text{LSF}}$ is the spectrum from the star after it is convolved by the instrumental LSF (this will be modeled by the empirical spectra of the star), $\alpha_s$ is a scaling term for the brightness of stellar speckles accounting for its fiber coupling efficiency, and $n$ is the noise (instrumental thermal background noise generally dominates). For simplicity, we will drop the $(\lambda)$ notation in the rest of the paper, but we note that all of these parameters will remain wavelength-dependent implicitly.

The transmission of the optical system was calculated using on-axis observations of the star HR 8799. In these observations, the signal $D_s$ can simply be written as

$$D_s = T S_{\text{LSF}} + n. \quad (2)$$

We approximated $S_{\text{LSF}}$ using a model PHOENIX spectrum with an effective temperature of 7200 K (Husser et al. 2013). Given that HR 8799 is an F0 star (Gray et al. 2003), it has nearly no spectral lines in the K band, which mitigates any errors due to an imperfect stellar spectrum. Additionally, we did not fit the data near the Brγ line in our analysis. The signal-to-noise ratio (S/N) per spectral channel of the stellar spectra is ~400, so $n$ was assumed to be negligible. The transmission is then simply obtained by dividing the data by the model of the star: $T = D_s/S_{\text{LSF}}$. We used $T$ solely to model the transmission of the planet light from the top of the atmospheres to the detector. To model the stellar light that leaked into the planet fiber, we simply used the on-axis observations of the star, $D_s$, to account for $T(\lambda) \times S(\lambda)$ simultaneously. This has the added benefit that the exact shape of the LSF only matters for the planet signal and is not used in fitting the stellar spectrum.

Thus, we constructed a model for the data obtained on the planet fiber ($M_p$) as

$$M_p = \alpha_p T P_{\text{LSF}} + \alpha_s D_s, \quad (3)$$

where we can use planetary atmosphere models for $P_{\text{LSF}}$ to find the model atmosphere parameters that best fit the data.

The terms $\alpha_p$ and $\alpha_s$ vary slowly as a function of wavelength. The wavelength dependence of $\alpha_p$ is dominated by differential atmospheric refraction changing the apparent sky position of the planet, and this effect changes slowly as a function of wavelength. The use of an atmospheric dispersion corrector in KPIC Phase II will mitigate this effect on $\alpha_p$ (Jovanovic et al. 2020; Wang et al. 2020b). We did not find chromatic optical aberrations in the system to be measurable for inclusion in $\alpha_p$. On similar observations of a speckle field with a single-mode fiber, Gravity Collaboration et al. (2020) showed that $\alpha_s$ can be approximated across the entire K band as a low-order polynomial. In some preliminary analyses to characterize $\alpha_p$ and $\alpha_s$ from observations

of standard stars, we confirmed this and found that they change on ~0.01 µm scales (~400 spectral channels). This means that high-pass filtering the data would be adequate for continuum subtraction and will remove chromatic modulations in the continuum due to $\alpha_p$ and $\alpha_s$.

To do this high-pass filter, we median-filtered the spectrum with a 200 pixel (~0.004 µm) box to compute the moving median and subtracted the median-filtered spectrum off from the original spectrum. We note that this continuum subtraction does not account for the fact that $\alpha_p$ and $\alpha_s$ also induce chromatic modulations in the line depths, which we will ignore in this work as it is a smaller effect. In principle, this effect can be modeled in our KPIC data using the low-pass-filter approach discussed in Ruffio et al. (2019).

We applied high-pass filtering to both the data and model spectra. The model spectrum after high-pass filtering can be written as

$$\mathcal{H}[M_p] \approx \mathcal{H}[\bar{\alpha}_p T P_{\text{LSF}} + \bar{\alpha}_s D_s], \quad (4)$$

where $\mathcal{H}$ denotes high-pass filtering, and $\bar{\alpha}_p$ and $\bar{\alpha}_s$ are wavelength-independent versions of $\alpha_p$ and $\alpha_s$. They represent achromatic scaling terms for the mean planetary and stellar speckle flux levels in the data. The high-pass filter of the data would be written in the same way.

In Section 4.2 (but not in Section 4.3), we made one further approximation assuming that the cross-talk between the various components is negligible when high-pass filtering using a median-filter implementation, allowing us to break down the model to filter the individual components:

$$\mathcal{H}[M_p] \approx \bar{\alpha}_p \mathcal{H}[T P_{\text{LSF}}] + \bar{\alpha}_s \mathcal{H}[D_s]. \quad (5)$$

We pulled $\bar{\alpha}_p$ and $\bar{\alpha}_s$ outside of $\mathcal{H}$ because we have assumed they are constant values after we removed the slowly varying chromatic continuum with the high-pass filter. We note that a linear implementation of the high-pass filter such as the Fourier-based approach from Ruffio et al. (2019) would make Equation (5) exact. However, we found that a median-filter implementation modeled the continuum better than the Fourier-based ones we tried (simple cutoff in frequency or a Gaussian filter in frequency space). Removing the chromatic continuum of $\alpha_p$ and $\alpha_s$ was more important than ensuring perfect linearity in the high-pass filter.

## 4.2. Detection of Molecules

First, we assessed the detection of all four planets by looking for the signature of molecular absorption lines from their planetary atmospheres as is common to do in the literature (e.g., Konopacky et al. 2013; Ruffio et al. 2019; Xuan et al. 2020). As shown in Ruffio (2019), cross-correlation techniques estimate the maximum likelihood value for the planet flux as a function of RV shift for a given planet template $P_{\text{LSF}}$. As we have spatially resolved data that allow us to construct the planet spectrum and stellar spectrum separately, we performed a modified cross-correlation as in Ruffio (2019) where we estimated the maximum likelihood value for both the planet and star flux as a function of RV shift for a given planet template. For a given RV shift, we can rewrite Equation (5) in matrix form to solve for the planet and star flux to best match
the high-pass-filtered data $\mathcal{H}[D_p]$:

$$ \left( \begin{array}{c} \vdots \\ \mathcal{H}[D_p] \\ \vdots \end{array} \right) = \left( \begin{array}{c} \vdots \\ \mathcal{H}[TP_{\text{LSF}}(\text{RV})] \\ \vdots \\ \mathcal{H}[D_s] \\ \vdots \end{array} \right) \left( \begin{array}{c} \bar{\alpha}_p \\ \bar{\alpha}_s \end{array} \right). \quad (6) $$

Here, the data on the left-hand side is a column vector with a length equal to the number of spectral channels, $N_\lambda$. The model on the right-hand side consists of an $N_\lambda \times 2$ matrix where the first column is the model flux of the planet and the second column is the model flux of the star (which we modeled empirically with the on-axis stellar spectrum) and a column vector of length 2 corresponding to the flux of the planet and star. The two unknowns are $\bar{\alpha}_p$ and $\bar{\alpha}_s$, making it straightforward to calculate their maximum likelihood values through linear least-squares optimization techniques. The cross-correlation function (CCF) for a given planet model $P_{\text{LSF}}$ is then the calculated values of $\bar{\alpha}_p$ at every requested RV. Note that we recomputed $P_{\text{LSF}}$ at every RV shift considered. When combining data from multiple echelle orders, the data and model columns can simply be extended so that $N_\lambda$ is the total number of spectral channels from all of the orders considered in the fit. While not included in this analysis, each data point can also be weighted by its respective error when estimating $\bar{\alpha}_p$ by using equation D6 of Ruffio et al. (2019) for a more optimal “matched filtering” detection metric. In the following analysis in this section, such error weighting made negligible changes to the CCF (<10%), so we have presented values without using it for simplicity.

For $P_{\text{LSF}}$, we generated molecular templates from the cloudless Sonora-Bobcat model set (Marley et al. 2018, 2021). Although this model lacked clouds, we found that the CCF signal from the cloudless Sonora-Bobcat CO+H$_2$O template was within 10% of the CCF signal from the cloudy BT-Settl models used in Section 4.3 that contained all molecular opacities, indicating that the strength of our molecular detections does not strongly depend on cloud assumptions. Note that molecular templates constructed from the BT-Settl models are not available. We used the same atmospheric model parameters for all four planets: an effective temperature of 1200 K and a surface gravity of $\log(g) = 4.0$ based on the latest values for HR 8799 e from Mollière et al. (2020). The detections were relatively insensitive to the exact choice of effective temperature and surface gravity: using a 1400 K model, which is representative of the scatter in effective temperature between different models (Greenbaum et al. 2018), changed the CCF signal by 5%.27 Using the temperature structure and molecular abundance profile from this Sonora-Bobcat model and following the procedure in Morley et al. (2015), we postprocessed the atmosphere model profile to compute emergent spectra at $R = 200,000$, repeating the process to only consider the opacity of the molecule or molecules of interest each time. The model opacities used were from Freedman et al. (2008, 2014), which utilized the CH$_4$ line lists from Yurchenko & Tennyson (2014) and H$_2$O line lists from Barber et al. (2006). The templates were then convolved to instrumental resolution assuming the spatial size of the fiber

trace measured in Section 3.1 is equivalent to the width of the LSF.

We computed CCFs using CO, H$_2$O, CH$_4$, and CO+H$_2$O molecular templates. Spectra of each planet from echelle orders 31 (2.44–2.49 µm), 32 (2.36–2.41 µm), and 33 (2.29–2.34 µm) were used, as these three orders had the best calibration. Including additional orders of the data or additional opacity sources in the model templates did not significantly increase the CCF signal. For each planet, we also computed CCFs for “noise” spectra that were taken contemporaneously. These noise spectra were either extracted from other fibers with similar amounts of stellar flux leaking into them as the science fiber but were not pointed at the planet or from regions of the detector that imaged the slit without any fibers and thus were dominated by thermal background noise. The CCFs were computed with velocity offsets between –500 and 500 km s$^{-1}$ from the solar system barycenter. The CCFs for each planet were normalized by dividing by the standard deviation of the CCFs from all the noise spectra, resulting in CCF S/N functions. The CCF S/N functions for each planet and for each molecular species are plotted in Figure 2. We note that the CCFs of the noise spectra may not be Gaussian-distributed, so there may be biases when quantifying the false-alarm probabilities of these molecular detections.

We found strong detections of CO for HR 8799 c, d, and e with S/N between 7 and 11. For these three planets, we also have a strong detection of H$_2$O with S/N > 5. For HR 8799 b, we have a weak detection with S/N ≈ 3 when combining the CO and H$_2$O templates, but there is no significant detection of any individual molecule alone. This is likely due to the fact HR 8799 b is about twice as faint as the other planets and had one of the shortest exposure times. Longer exposure times are needed to obtain confident detections of molecular signatures in HR 8799 b. For all four planets, we did not detect CH$_4$ at a significant level. The nondetection of CH$_4$ despite the strong CO and H$_2$O detections is consistent with molecular detections of these planets at medium resolution with OSIRIS (Konopacky et al. 2013; Barman et al. 2015; Petit dit de la Roche et al. 2018). Our KPIC detections are the first detections of HR 8799 b, d, and e at high spectral resolution ($R > 10,000$), where many molecular absorption features start becoming spectrally resolved. Our $6\sigma$ detection of H$_2$O in HR 8799 c in the K band is better than the previous $4.6\sigma$ L-band detection of H$_2$O by NIRSPAO that used $3.5 \times$ more integration time (Wang et al. 2018a).

### 4.3. Bayesian Inference of Planetary Parameters

We put our forward-modeling approach into a Bayesian inference framework to fit our extracted spectra directly, retrieve atmospheric parameters, and assess how complex a model is required to fit the data. Our likelihood-based framework is similar to the one developed by Ruffio et al. (2019) for medium-resolution spectroscopy of the HR 8799 planets, although we did not analytically marginalize over any of the parameters in our likelihood. Unlike Brogi & Line (2019), we did not use the CCF in the likelihood, as we did not assume that each spectral channel has the same noise level. Our method is similar to the method used by Gibson et al. (2020) to characterize the Fe absorption on an ultrahot Jupiter, except we simultaneously fit the star and planet together, which minimizes overfitting of the planet signal when subtracting off the star using techniques such as principal component analysis

<footnote n="27">27 Note, however, that the change in the likelihood of a model also depends on the number of data points so a small change in the CCF signal can significantly alter the likelihood of a model (Brogi & Line 2019).</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Cross-correlation functions with individual molecular species for all four HR 8799 planets using the technique described in Section 4.2. For each planet's spectrum (one planet per row), the S/Ns of the CCF are plotted for molecules templates using CO, H₂O, CH₄, and H₂O+CO (one per column). The solid color line is the cross-correlation function of each molecular template with the data on the planet. The solid gray lines are the cross-correlation function of each template with extracted spectra that contained no planet signal. The vertical dashed line marks the approximate systemic radial velocity as calculated by Ruffio et al. (2019). HR 8799 c, d, and e have strong detections of CO and H₂O individually, whereas HR 8799 b is only weakly detected when combining H₂O and CO molecular templates.</figcaption></figure>

(Pueyo 2016; Ruffio et al. 2019). We further note that the relative flux ratio between the planet and the stellar components in the data are much less extreme in our case, which likely makes it easier to fit both components to the data simultaneously.

We were interested in constraining the atmospheric and bulk parameters of the planets. For the planet spectrum, we fit the BT-Settl-CIFIST model grid to the data, varying both effective temperature ($T_{eff}$) and surface gravity (log(g) in cgs units; Allard et al. 2012). We chose the BT-Settl models as it was the only publicly available grid of models that are available at a spectral resolution $R > 35,000$ and includes clouds, which have been shown to be important in shaping these planets' spectra (e.g., Mollière et al. 2020). In particular for our high-resolution spectra, the depths of molecular absorption lines change when clouds are included (Hood et al. 2020). Rather than just stepping through RV shifts as in Section 4.2, we fit for the planetary RVs relative to the solar system barycenter. Note that the stellar RV is not well determined, making relative RV measurements challenging (Ruffio et al. 2019). Similar to Section 4.2, we can fit for the planet fluxes, which control the depth of the planetary atmosphere lines compared to the stellar and telluric line depths. This can also be directly translated to K-band flux ratios of the planets, which we can check against literature photometric measurements. We also fit for the

rotational broadening of the planetary spectra ($v \sin(i)$) by using the `fastRotBroad` function in `PyAstronomy` to broaden our planet atmosphere models (Czesla et al. 2019).

We also improved on the stellar model by considering multiple $D_s$ exposures. In Section 4.2, we averaged all of the on-axis images of the star in time. Here, we fit a set of linear coefficients $c_i$ to optimally weigh the individual on-axis exposures of the star, $D_{s,i}$, to create the master stellar spectrum that best fits the data:

$$D_s = \sum_i c_i D_{s,i} / \sum_i c_i. \qquad (7)$$

We note that this is similar to the LOCI technique in high-contrast imaging (Lafrenière et al. 2007), except that we optimized the coefficients while simultaneously fitting the planet model to the data, like what has been done in medium-resolution spectroscopy (Ruffio et al. 2019). In our fits, we assumed the $c_i$ coefficients to be unchanged across orders, but the overall spectrum $D_s$ can be scaled by a different flux value to account for the chromaticity of the stellar speckle flux.

Accurate estimation of errors is required for any robust statistical analysis. Preliminary analysis of the data indicated that the residuals to the forward-model fits are dominated by uncorrelated noise. However, we found that the amplitude of
the uncorrelated noise was higher than what the formal extraction errors predicted. This may be due to an underestimation of extraction errors or due to unaccounted-for noise terms by ~30%. Future work that performs more thorough analysis of the instrumental noise and data pipeline could help identify the source of this noise. For the purpose of this work, we simply fit for this excess uncorrelated noise, assuming it is Gaussian. Thus, the total error $\sigma_{tot}^2 = \sigma_{pipe}^2 + \sigma_{fit}^2$, where $\sigma_{pipe}$ is the nominal extraction error from our pipeline and $\sigma_{fit}$ is the error term we fit for. We used a separate $\sigma_{fit}$ for each spectral order but assumed that $\sigma_{fit}$ is constant within an order. This seemed to be a suitable approximation based on our analysis of the residuals to the fit.

In the analysis until now, we assumed that the width of the Gaussian trace of the fiber in the spatial dimension could be used as the LSF in the dispersion dimension. This assumption is not perfect as the NIRSPEC spectrograph was designed with a difference in focal lengths in the spatial and dispersion directions by a factor of 1.13 (Robichaud et al. 1998). In a preliminary analysis of the OH sky lines, which are unresolved at our resolution, we found that the LSF in the dispersion direction is indeed 1.12 ± 0.02 times wider than the width of the Gaussian profile measured in the spatial dimension. To conservatively account for any systematics in this preliminary measurement, we will allow the LSF width to vary between 1.0 and 1.2 times the width we measured in the spatial direction in Section 3.1 (i.e., the aspect ratio of the 2D LSF). This corresponds to varying the resolution from ~35,000 to ~29,000. We note that because the stellar spectrum is built using empirical data, the LSF size only affects the broadening of the planetary model, $P_{LSF}$.

We defined the log-likelihood to be

$$ \ln(\mathcal{L}) = -\frac{1}{2} \sum \left( \frac{(\mathcal{H}[D_p] - \mathcal{H}[M_p])^2}{\sigma_{tot}^2} + \ln(2\pi\sigma_{tot}^2) \right), \quad (8) $$

where we are summing over each spectral channel in the data. $\mathcal{H}[M_p]$ was constructed using Equation (4). Note that implicit in $M_p$ are the parameters of the planet, the instrumental LSF, and the nuisance parameters of the stellar speckle spectrum. For each planet, we considered three different models: a forward model that only contains the stellar speckle spectrum and no planet signature (i.e., $\bar{\alpha}_p = 0$ and all planet parameters are fixed), which we call “No Planet,” a model containing a planet with no discernible rotation ($v \sin(i) = 0$ is fixed) which we call “No Rotation,” and a rotating planet model where all planet parameters are allowed to vary which we call “Full Model.”

For HR 8799 b, c, and d, we used echelle orders 31 (2.44–2.49 µm), 32 (2.36–2.41 µm), and 33 (2.29–2.34 µm) to fit our model to the data just as was done in Section 4.2. We did not use orders 34–39 (1.95–2.29 µm) as they either had poor wavelength calibration (uncertainties larger than a spectral channel) or had strong CO$_2$ telluric absorption that was difficult to forward model accurately. For HR 8799 e, we only used orders 32 and 33 as we found that the strong telluric absorption features in order 31 could not be fully modeled with correlated residuals of comparable amplitude to the planet signal. We omitted this order to mitigate the effect of residual telluric lines from biasing our fit of the planetary atmosphere. Future work to marginalize over localized telluric residuals (Czekala et al. 2015) or improvements in modeling

the stellar speckle spectrum could make this order useful in the fit, but we chose to omit this order for now.

For the free parameters in each model, we used the following priors. For the planet properties, a uniform prior on $T_{eff}$ between 1000 and 1800 K, a uniform prior on $\log(g)$ between 3.5 and 5.5, a uniform prior on the planetary RV between −150 and 150 km s$^{-1}$, a uniform prior on $v \sin(i)$ in the Full Model between 0 and 60 km s$^{-1}$, and a uniform prior on the planet flux from 0 to 25 DN (data numbers). For each order, we included a parameter for the stellar speckle flux and nuisance parameters to fit for systematics in the data: a uniform prior on the stellar speckles flux between 0 and 500 DN, a uniform prior on a residual background term between −10 and 10 DN due to imperfect background subtraction, and a log-uniform prior on the $\sigma_{fit}$ between 0.1 and 30 DN. In the Full Model and No Rotation models, we also included for each order a multiplicative factor to account for any broadening of the LSF beyond what we measured as the spatial width of the fiber trace with a uniform prior between 1.0 and 1.2. This term was not in the No Planet model because that model consisted solely of empirical data where we did not need to broaden anything to instrumental resolution. To fit the stellar speckle spectrum using a linear combination of on-axis stellar spectra, we used a uniform prior between 0 and 2 for each $c_i$ term. If we term $N_s$ as the number of stellar spectra, $D_{s,i}$, used to compute $D_s$ and $N_{order}$ as the number of spectral orders used in the fit, the No Planet model has $3N_{order} + N_s$ free parameters, the No Rotation model has $4 + 4N_{order} + N_s$ parameters, and the Full Model has $5 + 4N_{order} + N_s$ free parameters. In the end, this resulted in ~20–30 free parameters, of which, we were mostly interested in the planetary parameters.

We sampled the posterior using the nested sampling algorithm (Skilling 2004, 2006) implemented in `pymultin-est`, which allowed us to both perform parameter estimation and compute the model evidence (Buchner et al. 2014). The evidence, $P(D|M)$, is used to compute the Bayes factor, $B$, that can be used to assess the relative probability of model $M_1$ compared to $M_2$. If $P(M_1|D)/P(M_2|D)$ expresses the relative probability of the two models given the data, then

$$ B \equiv \frac{P(D|M_1)}{P(D|M_2)} = \frac{P(M_1|D)P(M_2)}{P(M_2|D)P(M_1)}, \quad (9) $$

where $P(M)$ is the prior probability of that model. As we assumed each model has an equal prior probability, then $B$ is equivalent to the relative probability of two models. We used the Bayes factor to compare the simpler models to the Full Model to determine whether the additional free parameters are justified. Thus, our $M_2$ will always be the Full Model. We listed the estimates on the planet parameters and the Bayes factor of each model compared to the Full Model of that planet in Table 2. For the Full Model fits, we also plotted the posterior distributions of the key planet parameters in Figure 3. For the Full Model fits, the strongest covariance is between $T_{eff}$ and $\log(g)$, which we discuss in Section 5.1. We show corner plots of the posterior distributions from the Full Model fits in Appendix A.

For each planet, we can “decisively” reject any model with a $B$ that is more than 100 times smaller than the model with the highest $B$ based on the interpretation of the Bayes factor suggested by Jeffreys (1983). In Table 2, we see that the No Planet model is ruled out for c, d, and e, but remains 9% as
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Posterior distributions of the main planetary parameters for each planet in the Full Model fit. In the bottom-left panel, the gray shaded region represents the 68% credible interval of the radial velocity of the host star relative to the solar system barycenter.</figcaption></figure>

<table><caption>**Table 2**
Model Fits to HR 8799 KPIC Data</caption><thead><tr><th>Planet</th><th>Model</th><th>T<sub>eff</sub> (K)</th><th>log(g)</th><th>RV (km s<sup>-1</sup>)</th><th>v sin(i) (km s<sup>-1</sup>)</th><th>Planet Flux (DN)</th><th>B</th></tr></thead><tbody><tr><td>b</td><td>Full Model</td><td>1423.3<sup>+278.4</sup>(+362.3)<br/>-212.6(-354.0)</td><td>4.8<sup>+0.4</sup>(+0.7)<br/>-0.8(-1.2)</td><td>-15.6<sup>+6.5</sup>(+12.5)<br/>-7.1(-14.9)</td><td>31.2<sup>+11.9</sup>(+22.5)<br/>-11.7(-22.4)</td><td>5.4<sup>+1.9</sup>(+3.6)<br/>-1.7(-3.0)</td><td>1.0</td></tr><tr><td>b</td><td>No Rotation</td><td>1448.5<sup>+275.4</sup>(+339.1)<br/>-269.1(-399.4)</td><td>4.7<sup>+0.5</sup>(+0.8)<br/>-0.7(-1.1)</td><td>-9.8<sup>+3.2</sup>(+6.6)<br/>-5.9(-16.3)</td><td>...</td><td>3.4<sup>+1.4</sup>(+2.8)<br/>-1.1(-2.0)</td><td>0.13</td></tr><tr><td>b</td><td>No Planet</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>0.092</td></tr><tr><td>c</td><td>Full Model</td><td>1482.6<sup>+26.2</sup>(+52.2)<br/>+41.7(-118.0)</td><td>5.4<sup>+0.1</sup>(+0.1)<br/>-0.2(-0.6)</td><td>-12.4<sup>+0.8</sup>(+1.6)<br/>-1.1(-2.0)</td><td>8.1<sup>+3.8</sup>(+7.3)<br/>-4.0(-7.3)</td><td>9.5<sup>+1.1</sup>(+4.0)<br/>-0.8(-1.5)</td><td>1.0</td></tr><tr><td>c</td><td>No Rotation</td><td>1474.4<sup>+24.4</sup>(+38.9)<br/>-36.3(-102.6)</td><td>5.4<sup>+0.1</sup>(+0.1)<br/>-0.2(-0.4)</td><td>-12.4<sup>+0.9</sup>(+1.5)<br/>-1.0(-1.9)</td><td>...</td><td>9.0<sup>+0.7</sup>(+1.5)<br/>-0.6(-1.3)</td><td>2.1</td></tr><tr><td>c</td><td>No Planet</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>1.6 &times; 10<sup>-41</sup></td></tr><tr><td>d</td><td>Full Model</td><td>1558.8<sup>+50.9</sup>(+105.6)<br/>+81.4(-238.7)</td><td>5.1<sup>+0.3</sup>(+0.4)<br/>-0.4(-1.2)</td><td>-14.1<sup>+1.1</sup>(+2.2)<br/>-1.2(-2.2)</td><td>10.1<sup>+2.8</sup>(+5.7)<br/>-2.7(-6.6)</td><td>14.8<sup>+2.7</sup>(+4.8)<br/>-3.8(-6.4)</td><td>1.0</td></tr><tr><td>d</td><td>No Rotation</td><td>1642.2<sup>+108.8</sup>(+150.5)<br/>-127.1(-211.8)</td><td>5.4<sup>+0.1</sup>(+0.1)<br/>-0.2(-0.6)</td><td>-15.2<sup>+1.4</sup>(+2.7)<br/>-1.0(-1.8)</td><td>...</td><td>11.2<sup>+1.8</sup>(+3.8)<br/>-2.2(-4.1)</td><td>0.032</td></tr><tr><td>d</td><td>No Planet</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>3.1 &times; 10<sup>-29</sup></td></tr><tr><td>e</td><td>Full Model</td><td>1345.6<sup>+57.0</sup>(+124.8)<br/>-53.3(-99.5)</td><td>3.7<sup>+0.3</sup>(+0.9)<br/>-0.1(-0.2)</td><td>-12.3<sup>+1.2</sup>(+2.5)<br/>-1.4(-2.9)</td><td>15.0<sup>+2.3</sup>(+4.6)<br/>-2.6(-6.1)</td><td>16.2<sup>+4.5</sup>(+7.7)<br/>-3.6(-6.3)</td><td>1.0</td></tr><tr><td>e</td><td>No Rotation</td><td>1323.0<sup>+124.3</sup>(+318.3)<br/>-87.6(-163.8)</td><td>4.3<sup>+0.7</sup>(+1.1)<br/>-0.6(-0.8)</td><td>-14.6<sup>+1.1</sup>(+2.6)<br/>-1.1(-1.9)</td><td>...</td><td>9.2<sup>+1.7</sup>(+5.4)<br/>-1.3(-2.3)</td><td>0.049</td></tr><tr><td>e</td><td>No Planet</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>4.5 &times; 10<sup>-20</sup></td></tr></tbody></table>

Note. For each parameter, the median value is listed, with the subscript and superscript values representing the range of the central 68% credible interval with equal probability above and below the median (the central 95% credible interval is listed in parentheses).

likely as the Full Model for HR 8799 b. This finding is consistent with the marginal $3\sigma$ detection of b using the template cross-correlation in Section 4.2, but we are now able to assign relative probabilities to the cases. The Bayes factor between a planet and no planet model offers an alternative way to determine detection significance rather than using CCF S/N, where the false-alarm probability is unclear. Bayesian hypothesis testing methods have been used previously for source detection in cosmological data sets (Carvalho et al. 2009) and for exoplanet direct imaging (Golomb et al. 2019).

The data is less definitive in distinguishing between No Rotation and Full Model. In all cases, the No Rotation has a >1% probability compared to the Full Model. We know that it is unphysical to assume these planets have no spin, but the No Rotation model is a good approximation for a planet with a spin that remains undetectable. Thus, we interpret this result as the current data does not provide a definitive detection for rotational broadening. This could be due to a low S/N detection like in the case of HR 8799 b, but it could also be the difficulty of measuring small $v \sin(i)$ values, especially given
<figure><img src="image_7.png" /><figcaption>**Figure 4.** The 1D extracted spectra (black) from the fiber placed on each planet and the best-fitting forward model from echelle order 33. The forward model (blue) has been deconstructed into its two constituent parts: the stellar model (cyan) built from a linear combination of on-axis stellar spectra and the planet model generated from the BT-Settl atmospheric models (red). The residuals to the fit are plotted as gray circles and appear to be dominated by uncorrelated noise. A zoomed-in version of this plot is available in Appendix B.</figcaption></figure>

the near-face-on orbital configuration of the system. Still, HR 8799 d and e have a No Rotation model with $B < 0.05$, which Jeffreys (1983) interprets as “very strong” evidence for rotational broadening, the Full Planet model being 30 times more likely for d and 20 times more likely of e. We note that $v \sin(i)$ values from the Full Model for both HR 8799 d and e are inconsistent with 0 by >3σ, but the *B* metric downweighs this to account for the addition of this free parameter that could cause overfitting. *B* provides a more straightforward and more conservative assessment of the detection of rotational broadening rather than determining how far $v \sin(i)$ must deviate from 0 to be quantified as a detection because the median $v \sin(i)$ will always be nonzero by construction. For HR 8799 b

and c, we derive a 95% upper limit on the $v \sin(i)$ of 51 and 14 km s⁻¹, respectively, based on their marginalized 1D posteriors plotted in Figure 3. Thus, in this work, we were able to make the first measurements or constraints of the rotation velocities of the HR 8799 planets. In the future, higher-S/N detections or detections at higher spectral resolution would enable more definitive detections of rotation.

To assess the quality of our fits, we plotted the best-fitting set of parameters from the Full Model fit to one order of the data for all four planets in Figure 4. Visual inspection indicated our forward models fit the data adequately, with the residuals appearing to be dominated by uncorrelated noise. We verified this by computing the autocorrelation function (ACF) of the
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Atmospheric constraints on the HR 8799 planets. On the left are 1D and 2D marginalized posteriors of $T_{\text{eff}}$ and log(g) for HR 8799 c (red), d (green), and e (yellow). The 2D plot shows the 68% (solid) and 95% (dotted-dashed) credible region for each planet. A strong positive correlation between $T_{\text{eff}}$ and log(g) is observed. Black points and dotted lines correspond to representative models, which are plotted on the right.</figcaption></figure>

residuals and found that the ACF is well approximated by a delta function, with the wings of the ACF having an amplitude of ≲5% of the peak (see Appendix B). This confirms that an uncorrelated noise model is sufficient, as we are likely dominated by the thermal background noise of the instrument in our observations.

# 5. Discussion

## 5.1. Atmospheric Parameters

The detection of CO and H₂O but not CH₄ in our high-resolution spectra is consistent with previous atmospheric studies of the HR 8799 planets. Our CCF detections of HR 8799 c and d agree well with previous molecular cross-correlation analyses at lower resolution (Konopacky et al. 2013; Wang et al. 2018a). Due to the weak detection of HR 8799 b in our limited observations, we were only able to marginally detect the combined signal of CO and H₂O and did not have enough S/N to address previous disagreements on the amount of methane in its atmosphere (Barman et al. 2015; Petit dit de la Roche et al. 2018). Longer integration times and improvements to instrument performance will improve on these data. Regardless, this is the first time HR 8799 e is studied at a high spectral resolution, which demonstrates the high-contrast capabilities of KPIC. Previously, the highest-resolution spectrum for HR 8799 e was the R ~ 500 GRAVITY spectrum (Gravity Collaboration et al. 2019). For HR 8799 e, we found detections of CO and H₂O of similar strength as planets c and e and a similar nondetection of CH₄. This is consistent with the picture that the inner three planets have similar spectral signatures and luminosities based on lower-resolution spectroscopy.

We compared the planet fluxes measured in our forward-model fits and reported in Table 2 to the expected fluxes of the planets using the end-to-end system throughputs reported in Table 1, the exposure times of each frame, and the gain of the NIRSPEC detector (3.03 e-/DN; López et al. 2020). Using the photometry of the planets in the SPHERE K2 band from Zurlo et al. (2016), we expected to measure a flux of 4, 13, 16, and 17 DN for HR 8799 b, c, d, and e, respectively. These values are mostly consistent with our measured fluxes, although it appears that our inferred fluxes for HR 8799 c, d, and e are lower on average. This could be due to losses in flux due to a misaligned fiber when the blind offset was performed. It is not entirely clear this is the case, because we do not see a similar effect for

HR 8799 b, which had the largest blind offset. Still, the 95% credible intervals of the fluxes are consistent with the literature photometry for all four planets. The fact that the planet fluxes measured from our high-resolution data agree with broadband photometry increases our confidence that we have reliably separated the planet signal from the star.

The bulk atmospheric properties (effective temperature and surface gravity) we obtained from our forward-model fits somewhat agree with previous work at lower resolutions. Our $T_{\text{eff}}$ between 1400 and 1700 K for HR 8799 c and d agrees quite well with those obtained by Bonnefoy et al. (2016) and Greenbaum et al. (2018) in their BT-Settl fits, but HR 8799 e with a $T_{\text{eff}}$ between 1200 and 1500 K is lower than those works by 300 K (only 0.02% of our posterior have effective temperature >1650 K). We note that they obtained nonphysical radii in their BT-Settl fits, which pinpoint issues with fitting the BT-Settl grid to the HR 8799 spectral data (we did not perform absolute flux calibration so we did not constrain their radii). The log(g) values for HR 8799 e agree well with Bonnefoy et al. (2016) and Greenbaum et al. (2018). The posteriors for HR 8799 c and d favor log(g) > 4 and peak at the upper bound of the prior (i.e., upper bound of the model grid), which is at odds with those works that prefer log(g) < 4 (only 0.05% and 3% of our posteriors for HR 8799 c and d, respectively, had log(g) < 4). When looking beyond the BT-Settl model fits, our effective temperatures are systematically higher than the ~1100 K derived from other atmospheric models or predicted by evolutionary models, and the log(g) values for HR 8799 c and d continue to remain higher than literature values between 3.5 and 4 (Marois et al. 2008, 2010; Konopacky et al. 2013; Currie et al. 2014; Ingraham et al. 2014; Skemer et al. 2014; Bonnefoy et al. 2016; Greenbaum et al. 2018; Mollière et al. 2020). This indicates there may be some systematic errors in the BT-Settl models because both high- and low-resolution data found higher-$T_{\text{eff}}$ values for the BT-Settl model than all other models. One reason could be that BT-Settl models use equilibrium chemistry, and there has been evidence for disequilibrium chemistry in the HR 8799 planets causing stronger CO lines than expected and could make equilibrium chemistry models favor higher temperatures (Skemer et al. 2014; Mollière et al. 2020). However, we were not able to fit our data to other publicly available model grids with clouds as they are not generated at a sufficiently high spectral resolution to match our data.
HR 8799 c, d, and e are thought to have similar bulk atmospheric properties (e.g., luminosity, effective temperature, and radius) based on broadband spectroscopy and photometry (e.g., Marois et al. 2010; Bonnefoy et al. 2016). This appears to be at odds with our marginalized 1D posteriors, which prefer different $T_{\text{eff}}$ and log(g) values for HR 8799 e. Partially, this is due to a strong correlation between log(g) and $T_{\text{eff}}$: lower $T_{\text{eff}}$ values in our fits also prefer a lower log(g) (left panel of Figure 5). The contours for all three planets lie on nearly the same plane, suggesting that the differences in both $T_{\text{eff}}$ and log(g) may stem from additional parameters that we were not able to adjust in the model grid fits (e.g., cloud properties or disequilibrium chemistry). The need to adjust more parameters has been seen in BT-Settl fits to substellar companions where additional extinction parameters have been used to augment the default cloud prescription (Marocco et al. 2014; Bonnefoy et al. 2016; Delorme et al. 2017; Ward-Duong et al. 2021).

Looking more closely at why these fits to just our high-resolution spectra prefer these values, we plot how the BT-Settl spectra change as we change surface gravity and effective temperature in this correlated way in the right panel of Figure 5. We observed that the spectra are quite similar visually. The $T_{\text{eff}} = 1300$ K, log(g) = 3.5 model preferred by HR 8799 e has similar CO line depths in many of the lines as the $T_{\text{eff}} = 1650$ K, log(g) = 5.5 model preferred by HR 8799 d, whereas the intermediate-temperature models have the deepest CO lines. Because the BT-Settl models assumed solar abundances in carbon and oxygen, such changes in the H₂O and CO line depths could also be due to differences in abundance, given the broadband spectral data indicates the planets should have similar temperatures. Performing abundance measurements on this high-resolution data in the future will offer independent and sensitive measurements of the C/O ratios of the planets.

While we did not observe that planetary RV and $v \sin(i)$ are strongly correlated to $T_{\text{eff}}$ and log(g) (see Appendix A for corner plots), any inconsistencies in these bulk atmospheric properties could bias our inferred planetary RV and $v \sin(i)$ values. This should partially be mitigated by the fits having broad priors in $T_{\text{eff}}$ and log(g), which allowed us to marginalize over a large range of possible models in deriving 1D posteriors for planetary RV and $v \sin(i)$. However, there could be remaining systematic errors in the model that are not accounted for. To examine if there could be additional biases, we fixed $T_{\text{eff}}$ and log(g) to the literature best-fit values for the BT-Settl grid from Greenbaum et al. (2018) for HR 8799 c, d, and e and reran the fits. We found that planetary RV changed by <1σ and $v \sin(i)$ changed by <1.5σ for all three cases. We also assessed the sensitivity to cloud assumptions by rerunning the fits using the cloudless Sonora-Bobcat models instead and found that both planetary RV and $v \sin(i)$ changed by <1σ. These changes are comparable to our statistical errors, so any systematic biases should not significantly alter our results on RV and spin. We note that the B of these fits were <10⁻³ compared to the Full Model, so cloudless models are not good fits to the data, even if the planetary RV and $v \sin(i)$ appear to be consistent. Any smaller tweaks to the spectra (e.g., changing the C/O ratio) should have even smaller effects on the inferred planetary RV and $v \sin(i)$.

## 5.2. Orbital Velocity

The planetary RVs we measured are not precise enough to add significant constraints on their orbits but are consistent in amplitude and sign with previous work (Wang et al. 2018b;

Ruffio et al. 2019). Using the dynamically stable solutions from Wang et al. (2018b) and the measurement that the position angle of the ascending node is less than 180° (Ruffio et al. 2019), we expect an RV offset relative to the star of 1.9 ± 0.2 km s⁻¹ for planet b, 0.0 ± 0.2 km s⁻¹ for planet c, -2.7 ± 0.3 km s⁻¹ for planet d, and -2.1 ± 0.3 km s⁻¹ for planet e. We adopted the convention that negative RVs indicate the planet is moving toward us. The stellar RV is somewhat uncertain with the latest estimate being -10.5±0.5 km s⁻¹ with respect to the solar system barycenter (Ruffio et al. 2019).

Combining these two terms together, we have predicted planetary RVs relative to the solar system barycenter of -8.6±0.5 km s⁻¹ for planet b, -10.5±0.6 km s⁻¹ for planet c, -13.2±0.6 km s⁻¹ for planet d, and -12.6±0.7 km s⁻¹ for planet e. Although within the 95% credible interval of all our measurements, our values for HR 8799 b, c, and d are systematically more negative by ~1σ (~1 km s⁻¹). This could indicate an error in the absolute wavelength solution of our instrument or systematics in computing the stellar RV in previous analyses at the 1 km s⁻¹ level. Overall, these results confirm the ability to measure planetary RVs with KPIC.

## 5.3. Planetary Spin

We used our measured $v \sin(i)$ values to compare the rotational velocities for HR 8799 c, d, and e to the population of substellar companions and free-floating objects with measured spins. Although the spin measurements for HR 8799 c and d are not definitive, the strong evidence for rotation still gave us confidence to assess what these measurements would imply about their rotational evolution. We excluded HR 8799 b in the analysis because its $v \sin(i)$ is unconstrained given the tentative detection, and any preference in $v \sin(i)$ values may be spurious (we still reported some numbers taking the posterior at face value).

The bulk of the companion population was studied in Bryan et al. (2020b). We also did not consider HD 106906 b and AB Pic b in our subsequent analyses, because they only have tenuous rotation measurements (Zhou et al. 2019, 2020). We included the $v \sin(i)$ measured for GQ Lup b (Schwarz et al. 2016), the rotation period measured for HN Peg b (Zhou et al. 2016), and the rotation period of Jupiter (Dessler 1983) and Saturn (Helled et al. 2015) to extend the mass range of substellar companions considered. For the free-floating objects, we included the planetary-mass objects from Bryan et al. (2020b) as well as variable brown dwarfs with photometrically measured rotation periods compiled in Vos et al. (2020) and supplemented with measurements from Tannock et al. (2021). For the variable brown dwarfs, we only used those with parallax measurements from the literature (Faherty et al. 2012; Smart et al. 2013; Tinney et al. 2014; Liu et al. 2016; Gaia Collaboration et al. 2018, 2021; Theissen 2018; Best et al. 2020) in order to convert their K-band fluxes (Cutri et al. 2003) to masses using evolutionary tracks (Baraffe et al. 2003). We note that this sample contains objects formed through a variety of formation mechanisms: free-floating objects likely are the low-mass tail of cloud fragmentation (e.g., Luhman 2012), while planets like Jupiter and Saturn formed via core accretion (Pollack et al. 1996). Thus, depending on the dominant processes that produced the spins we measured, these populations may or may not obey similar spin trends. However, Bryan et al. (2020b) recently argued for a single spin regulation mechanism for both free-floating objects and companions.
<table><caption>**Table 3**

Derived Rotation Rate Values for the Hr 8799 Planets</caption><thead><tr><th rowspan="2">Planet</th><th rowspan="2">v<sub>break</sub></th><th colspan="4">Aligned to Orbit</th><th colspan="4">Random</th></tr><tr><th>v (km s<sup>-1</sup>)</th><th>Period (hr)</th><th>v/v<sub>break</sub></th><th>v<sub>final</sub>/v<sub>break</sub></th><th>v (km s<sup>-1</sup>)</th><th>Period (hr)</th><th>v/v<sub>break</sub></th><th>v<sub>final</sub>/v<sub>break</sub></th></tr></thead><tbody><tr><td>b</td><td>93<sup>+6(+13)</sup><sub>-6(-12)</sub></td><td>&lt;131</td><td>&gt;1.1</td><td>&lt;1.4</td><td>&lt;1.5</td><td>&lt;108</td><td>&gt;1.4</td><td>&lt;1.2</td><td>&lt;1.3</td></tr><tr><td>c</td><td>103<sup>+6(+14)</sup><sub>-7(-13)</sub></td><td>&lt;37</td><td>&gt;4.0</td><td>&lt;0.36</td><td>&gt;0.39</td><td>&lt;28</td><td>&gt;5.4</td><td>&lt;0.27</td><td>&gt;0.31</td></tr><tr><td>d</td><td>103<sup>+6(+14)</sup><sub>-7(-13)</sub></td><td>25<sup>+8(+17)</sup><sub>-7(-16)</sub></td><td>6.0<sup>+2.5(+1.3)</sup><sub>-1.5(-2.5)</sub></td><td>0.24<sup>+0.08(+0.17)</sup><sub>-0.07(-0.16)</sub></td><td>0.27<sup>+0.08(+0.18)</sup><sub>-0.08(-0.18)</sub></td><td>13<sup>+7(+36)</sup><sub>-4(-8)</sub></td><td>12.0<sup>+5.3(+22.6)</sup><sub>-4.4(-8.9)</sub></td><td>0.12<sup>+0.07(+0.34)</sup><sub>-0.04(-0.08)</sub></td><td>0.13<sup>+0.08(+0.39)</sup><sub>-0.04(-0.09)</sub></td></tr><tr><td>e</td><td>103<sup>+6(+14)</sup><sub>-7(-13)</sub></td><td>37<sup>+8(+18)</sup><sub>-7(-15)</sub></td><td>4.1<sup>+1.1(+3.3)</sup><sub>-0.8(-1.4)</sub></td><td>0.36<sup>+0.08(+0.18)</sup><sub>-0.08(-0.16)</sub></td><td>0.39<sup>+0.09(+0.19)</sup><sub>-0.08(-0.17)</sub></td><td>18<sup>+10(+47)</sup><sub>-4(-8)</sub></td><td>8.4<sup>+2.3(+6.4)</sup><sub>-3.0(-6.1)</sub></td><td>0.17<sup>+0.10(+0.47)</sup><sub>-0.04(-0.07)</sub></td><td>0.19<sup>+0.10(+0.53)</sup><sub>-0.04(-0.08)</sub></td></tr></tbody></table>

Note. For nondetections of rotation, the 95% upper limit or 5% lower limit is reported. Otherwise, the values reported follow the convention in Table 2.

Converting $v \sin(i)$ measurements requires accounting for the inclination degeneracy, and for most of these objects, including the HR 8799 planets, there are no empirical constraints on the spin axis inclination, and thus obliquity. While Vos et al. (2020) found a correlation between near-infrared color anomaly with spin axis inclination for low-gravity brown dwarfs, this relation is not particularly constraining for the HR 8799 planets: they have $(J-K) \sim 2.5$ mag (Zurlo et al. 2016), which corresponds to a $(J-K)$ color anomaly of $\sim -0.1$ mag (Liu et al. 2016), which lands in a region where inclinations between 20° and 90° were measured (Vos et al. 2020). Directly measuring the inclination of the spin axis typically requires a photometrically derived rotation period and a rotational broadening measurement, and an obliquity constraint needs the spin axis inclination and the orientation of the orbital plane. The planetary-mass companion 2M0122b has the only exoplanetary obliquity measurement to date and its obliquity may be high (Bryan et al. 2020a). A scenario that can produce this companion and give it a misaligned obliquity is formation via gravitational instability, where gravitoturbulence in the disk can torque fragment spin axes out of alignment (Bryan et al. 2020a). For the HR 8799 planets, formation via core accretion is preferred based on occurrence rate statistics and atmospheric compositions (Konopacky et al. 2013; Barman et al. 2015; Nielsen et al. 2019; Mollière et al. 2020). In this case, obliquity excitation may arise from chaotic spin-orbit dynamics due to the near commensurability of nodal and spin-precession frequencies (e.g., Neron de Surgy & Laskar 1997; Li & Batygin 2014; Shan & Li 2018; Saillenfest et al. 2019). Specifically, consider the planet $i=b, c, d$, and $e$, with mass $m_i$, radius $R_i$, spin frequency $\Omega_i = \bar{\Omega}_i/(Gm_i/R_i^3)^{1/2}$, and semimajor axis $a_i$ orbiting a host star of mass $M_*$. Because the nodal precession rate of planet $i$ due to the other planets is of order (e.g., Pu & Lai 2018)

$$ \omega_i \approx \sum_{j \neq i} \omega_{ij}, \omega_{ij} \approx \frac{3Gm_j a_i^2}{4a_i^3 \sqrt{GM_* a_i}}, \quad (10) $$

where $a_i = \min(a_i, a_j)$ and $a_i = \max(a_i, a_j)$, while the spin-precession rate of planet $i$ due to the host star's tidal torque is (assuming a fully convective planet, e.g., Lai 2014)

$$ \alpha_i \approx \frac{3GM_*}{4a_i^3} \left( \frac{R_i^3}{Gm_i} \right)^{1/2} \bar{\Omega}_i, \quad (11) $$

and all spin frequencies have magnitudes $\bar{\Omega}_i \sim 0.3$, we have, for the HR 8799 planets, $\omega_b/\omega_a \sim 8.1$, $\omega_c/\omega_a \sim 5.9$, $\omega_d/\omega_a \sim 2.8$, and $\omega_e/\omega_a \sim 0.6$, so all ratios of $\omega_i/\omega_a$ are of order unity. Hence, this system may be susceptible to secular spin-orbit resonances (e.g., Millholland & Laughlin 2019; Bryan et al. 2020b). The coupled

inclination and obliquity dynamics of the HR 8799 system will be investigated in a future work.

In the following analysis, we considered two bounding assumptions: the spin axes being randomly oriented in space as was done in Bryan et al. (2020b) and spin axes being aligned with the orbital planes based on orbital inclination posteriors derived from the literature (Ginski et al. 2014; Bryan et al. 2016, 2020a; Wang et al. 2018b; Pearce et al. 2019; Bowler et al. 2020; Nowak et al. 2020). We considered these to be bounding assumptions as reality likely does not have perfectly aligned obliquities nor a large fraction of retrograde obliquities. Two companions (SR 12 c and 2M0249b) did not have measurements of orbital motion or photometrically derived rotation periods, so the priors on their spin axes were assumed to be isotropic in both cases. Similarly, the free-floating objects were assumed to have isotropic spin axes orientations.

In both cases, we computed the rotational velocities ($v$), with a preference for using photometrically derived rotation periods over $v \sin(i)$ measurements because typically radius uncertainties are smaller than the inclination uncertainties. For the HR 8799 planets, we used an orbital inclination of $i = 24^\circ \pm 3^\circ$ for the zero-obliquity case (Wang et al. 2018b). We compared the rotation velocities against their breakup velocities, which we define simply as

$$ v_{\text{break}} = \sqrt{GM/R}, \quad (12) $$

where $M$ is the mass of an object and $R$ is its radius. This ignores the effect of oblateness, which will decrease the breakup velocity, with the exact decrease depending on the moment of inertia of the bodies (Marley & Sengupta 2011; Tannock et al. 2021). To derive breakup velocities for the HR 8799 planets, we used masses of $7.2 \pm 0.7 M_{\text{Jup}}$ for the inner three planets, a mass of $5.8 \pm 0.5 M_{\text{Jup}}$ for HR 8799 b, and radii of $1.2 \pm 0.1 R_{\text{Jup}}$ (Marois et al. 2008, 2010; Wang et al. 2018b). The values of $v_{\text{break}}$ for each planet are listed in Table 3. Under both obliquity assumptions, we computed $v$, rotation period, and $v/v_{\text{break}}$, and listed the values for the HR 8799 planets in Table 3. Note that we only reported 95% upper limits for HR 8799 b and c. Although no rotation periods have been measured from photometric monitoring of the HR 8799 planets (Apai et al. 2016; Biller et al. 2021), our inferred rotation periods are $>3$ hr for HR 8799 d and e under both obliquity assumptions. This is consistent with the rotation periods of PSO J318.5-22 and 2M1207b, the only objects similar in mass and age with photometrically measured rotation periods (Zhou et al. 2016; Biller et al. 2018). The rotation periods are also consistent with the picture that predicted magnetic braking from the circumplanetary disk regulating the spin of the planet
<figure><img src="image_9.png" /><figcaption>**Figure 6.** The final rotation rates relative to their breakup velocities as a function of mass for substellar objects. The left column assumes companion spins are aligned with their orbital planes, and the right column assumes randomly oriented spin axes. The top row only considers substellar companions excluding Jupiter and Saturn, the middle row includes Jupiter and Saturn, and the bottom row also includes free-floating substellar objects. In each plot, 50 possible power-law relations using the fits from Table 4 are plotted as purple lines with their dispersions shaded in blue. HR 8799 d is green and e is yellow. HR 8799 c is in red, but we note that one should consider this an upper limit. Other substellar companions with robust spin measurements are plotted in black. Free-floating objects are plotted in gray. Rotational velocities derived from photometrically measured rotation periods are unchanged between the two plots.</figcaption></figure>

at very early times ($\lesssim$1 Myr): Bryan et al. (2018, 2020b) found spins between 5% and 20% of breakup for planetary-mass companions, Batygin (2018) predicted a terminal rotation period prediction of $\sim$9 hr, and Ginzburg & Chiang (2020) placed a maximum rotation period of 29–43 hr.

To compare these rotation rates against the whole population of substellar objects with measured rotation rates, we needed to account for the fact that the rotation speed of an object depends on age. In our sample, ages ranged from $10^6$ to $10^9$ yr old. In the absence of outside influence, an isolated body will spin up as its radius contracts due to radiative cooling with $v \propto R^{-1}$ and $v/v_{\text{break}} \propto R^{-1/2}$. The current population of planetary-mass objects has been shown to be consistent with this trend (Bryan et al. 2020b). We evolved the radii of all companions from their current radii to the radii predicted by hot-start evolutionary models at 5 Gyr (Baraffe et al. 2003) and computed their final spin velocities relative to their final breakup velocities, $v_{\text{final}}/v_{\text{break}}$. The values of $v_{\text{final}}/v_{\text{break}}$ for the HR 8799 planets are listed in Table 3.

For the youngest objects in the sample, they may still host a circumplanetary disk (CPD), which combined with planetary magnetic fields is thought to regulate planetary spin (Batygin 2018; Ginzburg & Chiang 2020). CPD lifetimes are poorly constrained, so we do not know for sure which objects

are currently experiencing magnetic breaking. With evidence that CPDs and accretion occur for the PDS 70 planets (Haffert et al. 2019; Isella et al. 2019), which are $8 \pm 1$ Myr (Wang et al. 2020a), we assumed CPDs and regulation of planetary spin by magnetic braking occur up to 10 Myr. For companions younger than 10 Myr, we made one modification and assumed their rotational velocities are constant until 10 Myr. We note that magnetic braking is actually spinning down the planets, although the time dependence is weak (see Figure 1 of Ginzburg & Chiang 2020), so this is a rough approximation.

We compared the $v_{\text{final}}/v_{\text{break}}$ values of the HR 8799 planets against those of the rest of the gas giant and brown dwarf population in Figure 6 under our bounding obliquity assumptions and with three different subsets of the data (substellar companions excluding Jupiter and Saturn, substellar companions including Jupiter and Saturn, and all objects including free-floating substellar objects). Regardless of obliquity assumption and data subset, Figure 6 appears to show a trend of increasing rotational velocity (relative to breakup velocity) with decreasing companion mass. A similar trend in $v$ and mass for substellar objects has been pointed out in previous works (Snellen et al. 2014; Scholz et al. 2018), but no significant trend was found when looking at 1–20 $M_{\text{Jup}}$ planetary-mass objects previously (Bryan et al. 2018; Zhou et al. 2019; Xuan et al. 2020). Batygin (2018) and Ginzburg &
<table><caption>**Table 4**
Power-law Fits to the Rotation Rate of Gas Giant Planets</caption><thead><tr><th>Case</th><th>C</th><th>β</th><th>σ<sub>f</sub></th></tr><tr><th>Prior</th><th>LogUniform(0.01, 0.9)</th><th>Uniform(-3, 3)</th><th>LogUniform(0.01, 0.5)</th></tr></thead><tbody><tr><td><strong>Aligned to Orbit</strong></td><td></td><td></td><td></td></tr><tr><td>Substellar Companions (no Jupiter/Saturn)</td><td>0.20<sup>+0.03</sup>(+0.06)<br>0.02(-0.04)</td><td>-0.56<sup>+0.32</sup>(+0.66)<br>0.09(-0.21)</td><td>0.39<sup>+0.07</sup>(+0.10)<br>0.08(-0.14)</td></tr><tr><td>Substellar Companions (with Jupiter/Saturn)</td><td>0.20<sup>+0.03</sup>(+0.07)<br>0.02(-0.04)</td><td>-0.18<sup>+0.10</sup>(+0.18)<br>0.11(-0.24)</td><td>0.42<sup>+0.05</sup>(+0.08)<br>0.07(-0.14)</td></tr><tr><td>Including Free-floating Objects</td><td>0.19<sup>+0.02</sup>(+0.03)<br>0.01(-0.03)</td><td>-0.15<sup>+0.05</sup>(+0.09)<br>0.05(-0.11)</td><td>0.47<sup>+0.02</sup>(+0.03)<br>0.03(-0.07)</td></tr><tr><td><strong>Random</strong></td><td></td><td></td><td></td></tr><tr><td>Substellar Companions (no Jupiter/Saturn)</td><td>0.15<sup>+0.02</sup>(+0.04)<br>0.02(-0.03)</td><td>-0.61<sup>+0.25</sup>(+0.53)<br>0.27(-0.57)</td><td>0.33<sup>+0.09</sup>(+0.16)<br>0.07(-0.12)</td></tr><tr><td>Substellar Companions (with Jupiter/Saturn)</td><td>0.15<sup>+0.01</sup>(+0.03)<br>0.01(-0.03)</td><td>-0.30<sup>+0.07</sup>(+0.14)<br>0.09(-0.21)</td><td>0.33<sup>+0.08</sup>(+0.15)<br>0.06(-0.11)</td></tr><tr><td>Including Free-floating Objects</td><td>0.19<sup>+0.02</sup>(+0.03)<br>0.01(-0.02)</td><td>-0.14<sup>+0.05</sup>(+0.09)<br>0.05(-0.10)</td><td>0.47<sup>+0.02</sup>(+0.03)<br>0.03(-0.07)</td></tr></tbody></table>

**Note.** The values reported follow the convention in Table 2.

Chiang (2020) suggested that lower-mass planets may be less effective at ionizing their CPDs, making magnetic braking less effective in spinning down the planets and creating a spin speed that is mass dependent. To quantify the influence of mass on rotation rates, we followed the methodology from Wolfgang et al. (2016) where they used hierarchical Bayesian modeling to derive a probabilistic mass-radius relationship for transiting exoplanets. This approach accounts for uncertainties in both the independent and dependent variables by using their posteriors in the fit. This is important for our study as the constraints on mass and rotation speed can be weak and non-Gaussian, so point estimates such as the median and 68% credible interval do not provide the full picture. We fit a model that is a power law with intrinsic scatter of the following functional form:

$$ \frac{v_{\text{final}}}{v_{\text{break}}} \sim \text{Normal} \left( \mu = C \left( \frac{M}{10M_{\text{Jup}}} \right)^{\beta}, \sigma = \sigma_{f} \mu \right). \quad (13) $$

Following the notation of Wolfgang et al. (2016), the above indicates that the rotation rate is “drawn from the normal distribution” with a mean ($\mu$) and standard deviation ($\sigma$) that varies with companion mass ($M$). We fit for the power-law constant ($C$), the power-law index ($\beta$), and the fractional dispersion of the distribution at a given mass ($\sigma_f$). For both spin axis inclination assumptions and for the three different data subsets, we estimated these population-level parameters and marginalized over the priors on rotation velocities and masses for the individual companions using the *dynesty* nested sampling package (Speagle 2020) with multiple bounding ellipsoids (Feroz et al. 2009) and random slice sampling (Neal 2003; Handley et al. 2015a, 2015b). The priors and inferred values of the population-level parameters are listed in Table 4 and the fits excluding the solar system gas giants are plotted in Figure 6.

In all cases, we found a negative power-law index was preferred. With the exception of the case where we looked at the substellar companions excluding Jupiter and Saturn and assuming zero obliquity, the 95% credible interval of $\beta$ was entirely negative, although the mass dependence is weak: the 95% credible interval for $\beta$ is between -0.5 and 0 when including the solar system gas giants, and between -1.2 and +0.1 when excluding them. This could be a tentative sign that

the degree of ionization of the CPD by the planet itself affects the final rotation rate of a planet. However, there are other ways to create a mass dependence in $v/v_{\text{break}}$: for example, the effect of oblateness would cause $v_{\text{break}}$ to decrease by a factor that depends on the object's moment of inertia, which itself is also mass dependent (Marley & Sengupta 2011; Tannock et al. 2021). Even if it is true, the large scatter in the relation ($\sigma_f$ ranging from 25% to 50%) indicates there would still be other variables at play. For example, because the planet is free to spin up again after the CPD disperses in the magnetic braking scenario, different CPD lifetimes could also contribute to the scatter in the population (Bryan et al. 2020b). Such a scatter is seen in young free-floating brown dwarfs where a rotation period dispersion was linked to the presence or absence of a disk at ~8 Myr ages (Scholz et al. 2018; Moore et al. 2019). In our fits, including the free-floating objects caused $\sigma_f$ to reach the upper bound of the prior (50% dispersion), which may indicate that mass is not a driving factor in shaping the spins of the free-floating objects, or that they are affected in a different way, because they formed through a different process than planetary companions. As the spins of more companions in the 1–10 $M_{\text{Jup}}$ range are measured with instruments like KPIC, REACH, and HiRISE, we will be able to more confidently assess whether there is a trend in their rotation velocity with mass and if there are other factors at play.

## 6. Conclusions

We obtained high-spectral-resolution ($R \sim 35,000$) K-band spectra of all four HR 8799 planets taken as part of the science demonstration of KPIC, a new instrument that combines high-contrast imaging with high-resolution spectroscopy to enable HDC techniques. By cross-correlating the spectra with molecular templates, we detected both CO and H$_2$O at high significance (>10$\sigma$ combined) for HR 8799 c, d, and e, and tentatively detected HR 8799 b at $3\sigma$. These are the first detections of HR 8799 b, d, and e at high spectral resolution ($R \gtrsim 10,000$), where we can fully resolve the individual molecular lines in the planets' atmospheres. The detection of HR 8799 c has an S/N that is 2× better while using 3.5× less exposure time than the previous NIRSPAO detection (Wang et al. 2018a). With KPIC, we are able to access closer-in planets such as HR 8799 e (385 mas from its star), which has previously never been detected at either medium or high
resolution, demonstrating the HDC capabilities of KPIC. These
are the most challenging directly imaged exoplanets character-
ized with high-resolution spectroscopy to date given their flux
ratios and angular separations: while β Pic b (Snellen et al.
2014) is slightly closer in than HR 8799 e, it is an order-of-
magnitude brighter in the K band, which makes it easier to
isolate from the diffracted stellar light.

Beyond molecular detections via cross-correlation, we
developed a likelihood-based approach to jointly fit the stellar
speckles and planet light in our high-resolution data. We used
the Bayes factor computed from our framework to assess
detection probabilities that corroborated the molecular CCF
analysis: definitive detections of HR 8799 c, d, and e, but a
tentative detection of HR 8799 b where the No Planet model is
still 9% as likely as the best model fit with a planet. Using our
high-resolution data alone, we inferred bulk properties from
their atmospheres using the BT-Settl atmospheric models
(Allard et al. 2012). The measured RVs of the four planets are
consistent with previous orbital constraints (Wang et al. 2018b;
Ruffio et al. 2019). We found HR 8799 c, d, and e to have T<sub>eff</sub>
constrained to between 1200 and 1600 K, which is somewhat
consistent with previous works that used the BT-Settl models
(Bonnefoy et al. 2016; Greenbaum et al. 2018), but
significantly higher than other literature fits that used different
model grids. We also favored values of log(g) above 4.5 for
HR 8799 c and d, which are inconsistent with fits at low
spectral resolution and predictions of mass and radius from
evolutionary models. We speculated that deficiencies in the
BT-Settl models could be partially responsible for our inferred
parameters. Future work to study these planets with more
accurate atmospheric models would be able to determine
whether our results are due to model or data systematics, but
we were not able to perform such an analysis due to the limited
number of publicly available forward-model grids that both
model clouds and support a spectral resolution of 35,000.

We found strong evidence for rotational broadening in
our spectra of HR 8799 d and e, although more data are
needed to make them decisive detections. Taken at face value,
our inferred v sin(i) = 10.1<sup>+2.8</sup><sub>-2.7</sub> km s<sup>-1</sup> for HR 8799 d and
v sin(i) = 15.0<sup>+2.3</sup><sub>-2.6</sub> km s<sup>-1</sup> for HR 8799 e correspond to
0.24<sup>+0.08</sup><sub>-0.07</sub> v<sub>break</sub> and 0.36<sup>+0.08</sup><sub>-0.08</sub> v<sub>break</sub> assuming zero obliquities
and 0.12<sup>+0.07</sup><sub>-0.04</sub> v<sub>break</sub> and 0.17<sup>+0.10</sup><sub>-0.04</sub> v<sub>break</sub> assuming random obli-
quities. We found a rotational velocity upper limit for HR 8799
c of <14 km s<sup>-1</sup> corresponding to <0.36v<sub>break</sub> in the zero-
obliquity case, but the spin of HR 8799 b was essentially
unconstrained given the tentative detection. These spin
measurements are consistent with the picture of magnetic
braking regulating the spins of gas giant planets during the
planet formation process. When combining these rotation
measurements with literature rotation measurements of sub-
stellar companions as well as our solar system gas giants, we
found a tentative trend of increasing rotation rate with
decreasing planet mass that could point to magnetic braking
being less efficient for lower-mass planets.

6.1. *Future Prospects*

More spin measurements of 1–10 $M_{\text{Jup}}$ companions will test
whether this mass–spin relation holds, or if there are other key
parameters to control the spin rate of giant planets. In our
analysis, we assumed two bounding assumptions on the

obliquity of the planets to derive true rotation rates. Obliquity
constraints on the HR 8799 planets themselves may be possible
in the future if their periods can be derived from rotational
modulations in their light curves. This could be possible with
new hardware designed to improve the photometric calibration
of coronagraphic systems (Bos 2020).

Future analysis work that moves beyond molecular detection and characterization of bulk parameters will provide more insights into the HR 8799 planets. Retrievals of the chemical abundances of their atmospheres have been applied to broad-band photometry and lower-resolution spectra (Lavie et al. 2017; Mollière et al. 2020). The likelihood framework we constructed to fit the high-spectral-resolution KPIC data would be straightforward to include in these Bayesian retrievals, allowing us to make use of the information at all spectral resolutions in understanding the chemical makeup of these planets. KPIC is also being used to conduct a spectroscopic survey of substellar companions that were previously too close to their bright host stars to study at high spectral resolution with slit spectrographs. Measuring planet characteristics like spin, orbital parameters, and composition across this population will allow us to look for trends that would pinpoint the dominant processes in their formation and evolution.

Lastly, future planned upgrades to KPIC (Pezzato et al.
2019; Jovanovic et al. 2020) are designed to improve the
coupling of planet light, reduce the coupling of starlight, and
mitigate the instrument thermal background. As we are limited
by thermal background noise currently, increasing the amount
of planet light reaching the detector and reducing the amount of
thermal background photons will directly translate into
improvements in the S/N of KPIC observations. Reobserving
the HR 8799 planets with these upcoming upgrades would give
us more precise planetary spin measurements and improved
sensitivity to any trace methane abundances in their atmo-
spheres to constrain the level of disequilibrium chemistry. The
initial phase of KPIC is the bare-bones hardware necessary to
enable HDC techniques, and refined HDC techniques will
allow for detailed characterization of all of the directly imaged
exoplanets.

We thank Sivan Ginzburg, Konstantin Batygin, and Eve Lee
for useful discussions. We thank the Keck staff for helping
support the deployment of KPIC during the global pandemic.
J.J.W. is supported by the Heising-Simons Foundation 51
Pegasi b Fellowship. This work was supported by the Heising-
Simons Foundation through grants #2015-129, #2017-318,
and #2019-1312. This work was supported by the Simons
Foundation. E.C.M. is supported by an NSF Astronomy and
Astrophysics Postdoctoral Fellowship under award AST-
1801978. Part of the computations presented here were
conducted on the Caltech High Performance Cluster, partially
supported by a grant from the Gordon and Betty Moore
Foundation. Part of this work was carried out at the Jet
Propulsion Laboratory, California Institute of Technology,
under contract with NASA. Data presented in this work were
obtained at the W. M. Keck Observatory, which is operated as
a scientific partnership among the California Institute of
Technology, the University of California, and the National
Aeronautics and Space Administration. The Observatory was
made possible by the generous financial support of the W. M.
Keck Foundation. We wish to recognize and acknowledge the
very significant cultural role and reverence that the summit of
Maunakea has always had within the indigenous Hawaiian
community. We are most fortunate to have the opportunity to
conduct observations from this mountain.

Facility: Keck (KPIC).

Software: astropy (Astropy Collaboration et al. 2018), scipy (Virtanen et al. 2020), PyAstronomy (Czesla et al. 2019), pymultinest (Buchner et al. 2014), dynasty (Speagle 2020), whereistheplanet (Wang et al. 2021), corner (Foreman-Mackey 2016).

**Appendix A**

Corner Plots for Full Model Fits

We present corner plots for the Full Model fits for each
planet: HR 8799 b in Figure 7, HR 8799 c in Figure 8, HR
8799 d in Figure 9, and HR 8799 e in Figure 10. For each
planet, we only showed six parameters, which are the main
parameters of interest. These included the planet bulk proper-
ties T<sub>eff</sub>, log(g), planetary RV, v sin(i), and planetary flux. We
also showed the measured speckle flux in Order 33. We did not
show the other orders for simplicity as they have the same
behavior. We also did not show the error inflation, offset, LSF size

<figure><img src="image_10.png" /><figcaption>**Figure 7.** Corner plot for HR 8799 b for the posterior distribution of key parameters in the Full Model fit.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>**Figure 8.** Corner plot for HR 8799 c for the posterior distribution of key parameters in the Full Model fit.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>**Figure 9.** Corner plot for HR 8799 d for the posterior distribution of key parameters in the Full Model fit.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>Figure 10. Corner plot for HR 8799 e for the posterior distribution of key parameters in the Full Model fit.</figcaption></figure>

scaling, or linear coefficients $c_i$ for combining together the stellar spectra as they are nuisance parameters. The inclusion of all of these parameters would have limited the readability of these plots.

## Appendix B
### Forward-model Fit and Residuals

To better assess the fits in detail, we created a version of Figure 4 that is zoomed in to a 7 nm chunk of the Full Model fit in Order 33 (Figure 11). We can see good agreement between the data and forward model on a channel-to-channel basis.

Another way to assess the fits are satisfactory is to look for any correlated structure in the residuals. This could be due to model mismatch or correlated noise sources that we had not accounted for. Although we did not visually see any structures in the residuals, the fact there are ~2000 data points per order means we can look for lower levels of correlated residual structure by computing the ACF:

$$ACF(x) = \frac{\sum_{i} R_{i} R_{i+x}}{\sum_{i} R_{i} R_{i}}. \qquad (B1)$$
<figure><img src="image_14.png" /><figcaption>**Figure 11.** A version of Figure 4, except zoomed in to a 7 nm chunk of the fit. The gaps in the data and models correspond to channels that were masked due to bad pixels.</figcaption></figure>

Here, *R* is the 1D residual vector of length equal to the number of data points in the order, *x* is the lag in pixels (lag of zero results in ACF ≡ 1), and Σᵢ iterates over all elements of *R*. For an infinitely long sequence of uncorrelated Gaussian noise, the ACF is the Kronecker delta function. Any correlated structure

in the residuals would result in nonzero values in the wings of
the ACF.

From the fits in Section 4.3, we computed the fit residuals for
the parameter set from the Full Model fit of each planet with the
highest likelihood. For each echelle order used, we computed
<figure><img src="image_15.png" /><figcaption>**Figure 12.** ACF of the fit residuals for the best-fitting model for each planet. For each planet, all of the orders we fit to are plotted. The lack of discernible power in the wings of ACF for each order indicates we are dominated by uncorrelated noise.</figcaption></figure>

the ACF of these residuals and plotted them in Figure 12. From visual inspection of the ACFs, we see that they are almost entirely consistent with a delta function. For HR 8799 e, we saw some power in the wings, indicating small correlated structures in the residual (≤5% of the total scatter in the residuals).

## ORCID iDs

* Jason J. Wang (王劲飞) https://orcid.org/0000-0003-0774-6502

* Jean-Baptiste Ruffio https://orcid.org/0000-0003-2233-4821

* Jacques-Robert Delorme https://orcid.org/0000-0001-8953-1008

* Nemanja Jovanovic https://orcid.org/0000-0001-5213-6207

* Luke Finnerty https://orcid.org/0000-0002-1392-0768

* Callie Hood https://orcid.org/0000-0003-1150-7889

* J. J. Zanazzi https://orcid.org/0000-0002-9849-5886

* Marta L. Bryan https://orcid.org/0000-0002-6076-5967

* Emily C. Martin https://orcid.org/0000-0002-0618-5128

* Dimitri Mawet https://orcid.org/0000-0002-8895-4735

* Andy Skemer https://orcid.org/0000-0001-6098-3924

* Jerry W. Xuan https://orcid.org/0000-0002-6618-1137

* Geoffrey A. Blake https://orcid.org/0000-0003-0787-1610

* Cam Buzard https://orcid.org/0000-0002-9943-6124

* Mark Chun https://orcid.org/0000-0002-8462-0703

* Greg Doppmann https://orcid.org/0000-0001-9823-1445

* Trent J. Dupuy https://orcid.org/0000-0001-9823-1445

* Gaspard Duchêne https://orcid.org/0000-0002-5092-6464

* Michael P. Fitzgerald https://orcid.org/0000-0002-0176-8973

* Jonathan Fortney https://orcid.org/0000-0002-9843-4354

* Richard S. Freedman https://orcid.org/0000-0001-9333-4306

* Heather Knutson https://orcid.org/0000-0002-5375-4725

* Quinn Konopacky https://orcid.org/0000-0002-9936-6285

* Michael C. Liu https://orcid.org/0000-0003-2232-7664

* Roxana Lupu https://orcid.org/0000-0003-3444-5908

* Mark S. Marley https://orcid.org/0000-0002-5251-2943

* Tiffany Meshkat https://orcid.org/0000-0001-6126-2467

* Brittany Miles https://orcid.org/0000-0002-5500-4602

* Maxwell Millar-Blanchaer https://orcid.org/0000-0001-6205-9233

* Arpita Roy https://orcid.org/0000-0001-8127-5775

* Greeth Ruane https://orcid.org/0000-0003-4769-1665

* Peter Wizinowich https://orcid.org/0000-0002-1646-442X

## References

* Allard, F., Homeier, D., & Freytag, B. 2012, RSPTA, 370, 2765

* Apai, D., Kasper, M., Skemer, A., et al. 2016, ApJ, 820, 40

* Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123

* Baraffe, I., Chabrier, G., Barman, T. S., Allard, F., & Hauschildt, P. H. 2003, A&A, 402, 701

* Barber, R. J., Tennyson, J., Harris, G. J., & Tolchenov, R. N. 2006, MNRAS, 368, 1087

* Barman, T. S., Konopacky, Q. M., Macintosh, B., & Marois, C. 2015, ApJ, 804, 61
Barman, T. S., Macintosh, B., Konopacky, Q. M., & Marois, C. 2011, ApJ, 733, 65

Batygin, K. 2018, AJ, 155, 178

Best, W. M. J., Liu, M. C., Magnier, E. A., & Dupuy, T. J. 2020, AJ, 159, 257

Biller, B. A., Apai, D., Bonnefoy, M., et al. 2021, MNRAS, 503, 743

Biller, B. A., Vos, J., Buenzli, E., et al. 2018, AJ, 155, 95

Bond, C. Z., Cetre, S., Lilley, S., et al. 2020, JATIS, 6, 039003

Bonnefoy, M., Zurlo, A., Baudino, J. L., et al. 2016, A&A, 587, A58

Bos, S. P. 2020, A&A, 638, A118

Boss, A. P. 1998, ApJ, 503, 923

Boss, A. P. 2001, ApJL, 551, L167

Bowler, B. P. 2016, PASP, 128, 102001

Bowler, B. P., Blunt, S. C., & Nielsen, E. L. 2020, AJ, 159, 63

Bowler, B. P., Liu, M. C., Dupuy, T. J., & Cushing, M. C. 2010, ApJ, 723, 850

Brogi, M., & Line, M. R. 2019, AJ, 157, 114

Bryan, M. L., Benneke, B., Knutson, H. A., Batygin, K., & Bowler, B. P. 2018, NatAs, 2, 138

Bryan, M. L., Bowler, B. P., Knutson, H. A., et al. 2016, ApJ, 827, 100

Bryan, M. L., Chiang, E., Bowler, B. P., et al. 2020a, AJ, 159, 181

Bryan, M. L., Ginzburg, S., Chiang, E., et al. 2020b, ApJ, 905, 37

Buchner, J., Georgakakis, A., Nandra, K., et al. 2014, A&A, 564, A125

Carvalho, P., Rocha, G., & Hobson, M. P. 2009, MNRAS, 393, 681

Currie, T., Burrows, A., Girard, J. H., et al. 2014, ApJ, 795, 133

Cutri, R. M., Skrutskie, M. F., van Dyk, S., et al. 2003, yCat, II, 246

Czekala, I., Andrews, S. M., Mandel, K. S., Hogg, D. W., & Green, G. M. 2015, ApJ, 812, 128

Czesla, S., Schröter, S., Schneider, C. P., et al. 2019, PyA: Python Astronomy-Related Packages, Astrophysics Source Code Library, ascl:1906.010

Delorme, J.-R., Jovanovic, N., Echeverri, D., et al. 2021, JATIS, 7, 035006

Delorme, P., Schmidt, T., Bonnefoy, M., et al. 2017, A&A, 608, A79

Dessler, A. J. 1983, Coordinate Systems (Cambridge: Cambridge Univ. Press), 498

Faherty, J. K., Burgasser, A. J., Walter, F. M., et al. 2012, ApJ, 752, 56

Feroz, F., Hobson, M. P., & Bridges, M. 2009, MNRAS, 398, 1601

Foreman-Mackey, D. 2016, JOSS, 1, 24

Freedman, R. S., Lustig-Yaeger, J., Fortney, J. J., et al. 2014, ApJS, 214, 25

Freedman, R. S., Marley, M. S., & Lodders, K. 2008, ApJS, 174, 504

Gaia Collaboration, Brown, A. G. A., & Vallenari, A. 2018, A&A, 616, A1

Gaia Collaboration, Brown, A. G. A., & Vallenari, A. 2021, A&A, 649, A1

Gibson, N. P., Merritt, S., Nugroho, S. K., et al. 2020, MNRAS, 493, 2215

Ginski, C., Schmidt, T. O. B., Mugrauer, M., et al. 2014, MNRAS, 444, 2280

Ginzburg, S., & Chiang, E. 2020, MNRAS, 491, L34

Golomb, J., Rocha, G., Meshkat, T., et al. 2019, arXiv:1912.01232

Goździewski, K., & Migaszewski, C. 2020, ApJL, 902, L40

Gravity Collaboration, Lacour, S., Nowak, M., et al. 2019, A&A, 623, L11

Gravity Collaboration, Nowak, M., Lacour, S., et al. 2020, A&A, 633, A110

Gray, R. O., Corbally, C. J., Garrison, R. F., McFadden, M. T., & Robinson, P. E. 2003, AJ, 126, 2048

Greenbaum, A. Z., Pueyo, L., Ruffio, J.-B., et al. 2018, AJ, 155, 226

Haffert, S. Y., Bohn, A. J., de Boer, J., et al. 2019, NatAs, 3, 749

Handley, W. J., Hobson, M. P., & Lasenby, A. N. 2015a, MNRAS, 450, L61

Handley, W. J., Hobson, M. P., & Lasenby, A. N. 2015b, MNRAS, 453, 4384

Helled, R., Galanti, E., & Kaspi, Y. 2015, Natur, 520, 202

Hood, C. E., Fortney, J. J., Line, M. R., et al. 2020, AJ, 160, 198

Horne, K. 1986, PASP, 98, 609

Husser, T. O., Wende-von Berg, S., Dreizler, S., et al. 2013, A&A, 553, A6

Ingraham, P., Marley, M. S., Saumon, D., et al. 2014, ApJL, 794, L15

Isella, A., Benisty, M., Teague, R., et al. 2019, ApJ, 879, L25

Jeffreys, H. 1983, Theory of Probability (Oxford: Clarendon), 459

Jovanovic, N., Calvin, B., Porter, M., et al. 2020, Proc. SPIE, 11447, 114474U

Jovanovic, N., Delorme, J.-R., Bond, C. Z., et al. 2019, Proc. SPIE, 11117, 111170T

Jovanovic, N., Guyon, O., Kotani, T., et al. 2017, arXiv:1712.07762

Konopacky, Q. M., Barman, T. S., Macintosh, B. A., & Marois, C. 2013, Sci, 339, 1398

Kotani, T., Kawahara, H., Ishizuka, M., et al. 2020, Proc. SPIE, 11448, 1144878

Lafrenière, D., Marois, C., Doyon, R., Nadeau, D., & Artigau, E. 2007, ApJ, 660, 770

Lai, D. 2014, MNRAS, 440, 3532

Lavie, B., Mendonça, J. M., Mordasini, C., et al. 2017, AJ, 154, 91

Li, G., & Batygin, K. 2014, ApJ, 790, 69

Liu, M. C., Dupuy, T. J., & Allers, K. N. 2016, ApJ, 833, 96

López, R. A., Hoffman, E. B., Doppmann, G., et al. 2020, Proc. SPIE, 11447, 114476B

Lord, S. D. 1992, A New Software Tool for Computing Earth's Atmospheric Transmission of Near- and Far-Infrared Radiation, Technical Report, 103957

Luhman, K. L. 2012, ARA&A, 50, 65

Marley, M., Saumon, D., Morley, C., & Fortney, J. 2018, Sonora 2018: Cloud-free, Solar Composition, Solar C/O Substellar Atmosphere Models and Spectra, nc_m+0.0_co1.0_v1.0, Zenodo, doi:10.5281/zenodo.1309035

Marley, M. S., Samon, D., Visscher, C., et al. 2021, arXiv:2107.07434

Marley, M. S., & Sengupta, S. 2011, MNRAS, 417, 2874

Marocco, F., Day-Jones, A. C., Lucas, P. W., et al. 2014, MNRAS, 439, 372

Marois, C., Macintosh, B., Barman, T., et al. 2008, Sci, 322, 1348

Marois, C., Zuckerman, B., Konopacky, Q. M., Macintosh, B., & Barman, T. 2010, Natur, 468, 1080

Martin, E. C., Fitzgerald, M. P., McLean, I. S., et al. 2014, Proc. SPIE, 9147, 914781

Martin, E. C., Fitzgerald, M. P., McLean, I. S., et al. 2018, Proc. SPIE, 10702, 107020A

Mawet, D., Delorme, J. R., Jovanovic, N., et al. 2017a, Proc. SPIE, 10400, 1040029

Mawet, D., Ruane, G., Xuan, W., et al. 2017b, ApJ, 838, 92

McLean, I. S., Becklin, E. E., Bendiksen, O., et al. 1998, Proc. SPIE, 3354, 566

Millholland, S., & Laughlin, G. 2019, NatAs, 3, 424

Mollière, P., Stolker, T., Lacour, S., et al. 2020, A&A, 640, A131

Moore, K., Scholz, A., & Jayawardhana, R. 2019, ApJ, 872, 159

Morley, C. V., Fortney, J. J., Marley, M. S., et al. 2015, ApJ, 815, 110

Morris, E. C., Wang, J. J., Ruffio, J.-B., et al. 2020, Proc. SPIE, 11447, 1144761

Neal, R. M. 2003, AnSta, 31, 705

Neron de Surgy, O., & Laskar, J. 1997, A&A, 318, 975

Nielsen, E. L., De Rosa, R. J., Macintosh, B., et al. 2019, AJ, 158, 13

Nowak, M., Lacour, S., Lagrange, A.-M., et al. 2020, A&A, 642, L2

Otten, G. P. P. L., Vigan, A., Muslimov, E., et al. 2021, A&A, 646, A150

Pearce, L. A., Kraus, A. L., Dupuy, T. J., et al. 2019, AJ, 157, 71

Petit dit de la Roche, D. J. M., Hoeijmakers, H. J., & Snellen, I. A. G. 2018, A&A, 616, A146

Pezzato, J., Jovanovic, N., Mawet, D., et al. 2019, Proc. SPIE, 11117, 111170U

Pollack, J. B., Hubickyj, O., Bodenheimer, P., et al. 1996, Icar, 124, 62

Pu, B., & Lai, D. 2018, MNRAS, 478, 197

Pueyo, L. 2016, ApJ, 824, 117

Robichaud, J. L., Zellers, B., Philippon, R., McLean, I. S., & Figer, D. F. 1998, Proc. SPIE, 3354, 1068

Ruffio, J.-B. 2019, PhD thesis, Stanford University, http://purl.stanford.edu/ yq071tj0740

Ruffio, J.-B., Macintosh, B., Konopacky, Q. M., et al. 2019, AJ, 158, 200

Saillenfest, M., Laskar, J., & Boué, G. 2019, A&A, 623, A4

Scholz, A., Moore, K., Jayawardhana, R., et al. 2018, ApJ, 859, 153

Schwarz, H., Ginski, C., de Kok, R. J., et al. 2016, A&A, 593, A74

Serabyn, E., Huby, E., Matthews, K., et al. 2017, AJ, 153, 43

Shan, Y., & Li, G. 2018, AJ, 155, 237

Skemer, A. J., Marley, M. S., Hinz, P. M., et al. 2014, ApJ, 792, 17

Skilling, J. 2004, in AIP Conf. Ser. 735, Bayesian Inference and Maximum Entropy Methods in Science and Engineering: 24th Int. Workshop on Bayesian Inference and Maximum Entropy Methods in Science and Engineering, ed. R. Fischer, R. Preuss, & U. V. Toussaint (Melville, NY: AIP), 395

Skilling, J. 2006, BayAn, 1, 833

Smart, R. L., Tinney, C. G., Bucciarelli, B., et al. 2013, MNRAS, 433, 2054

Snellen, I., de Kok, R., Birkby, J. L., et al. 2015, A&A, 576, A59

Snellen, I. A. G., Brandl, B. R., de Kok, R. J., et al. 2014, Natur, 509, 63

Soubiran, C., Jasnowicz, G., Chemin, L., et al. 2018, A&A, 616, A7

Speagle, J. S. 2020, MNRAS, 493, 3132

Tannock, M. E., Metchev, S., Heinze, A., et al. 2021, AJ, 161, 224

Theissen, C. A. 2018, ApJ, 862, 173

Tinney, C. G., Faherty, J. K., Kirkpatrick, J. D., et al. 2014, ApJ, 796, 39

Vargas Catalán, E., Huby, E., Forsberg, P., et al. 2016, A&A, 595, A127

Vigan, A., Fontanive, C., Meyer, M., et al. 2020, A&A, 651, A72

Vigan, A., Otten, G. P. P. L., Muslimov, E., et al. 2018, Proc. SPIE, 10702, 1070236

Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, NatMe, 17, 261

Vos, J. M., Biller, B. A., Allers, K. N., et al. 2020, AJ, 160, 38

Wang, J., Mawet, D., Fortney, J. J., et al. 2018a, AJ, 156, 272

Wang, J. J., Ginzburg, S., Ren, B., et al. 2020a, AJ, 159, 263

Wang, J. J., Graham, J. R., Dawson, R., et al. 2018b, AJ, 156, 192

Wang, J. J., Kulikauskas, M., & Blunt, S. 2021, Whereistheplanet: Predicting Positions of Directly Imaged Companions, Astrophysics Source Code Library, ascl:2101.003
Wang, J. J., Wallace, J. K., Jovanovic, N., et al. 2020b, Proc. SPIE, 11447, 1144754

Ward-Duong, K., Patience, J., Follette, K., et al. 2021, AJ, 161, 5

Wilcomb, K. K., Konopacky, Q. M., Barman, T. S., et al. 2020, AJ, 160, 207

Wolfgang, A., Rogers, L. A., & Ford, E. B. 2016, ApJ, 825, 19

Xuan, J. W., Bryan, M. L., Knutson, H. A., et al. 2020, AJ, 159, 97

Yurchenko, S. N., & Tennyson, J. 2014, MNRAS, 440, 1649

Zhou, Y., Apai, D., Bedin, L. R., et al. 2020, AJ, 159, 140

Zhou, Y., Apai, D., Lew, B. W. P., et al. 2019, AJ, 157, 128

Zhou, Y., Apai, D., Schneider, G. H., Marley, M. S., & Showman, A. P. 2016, ApJ, 818, 176

Zurlo, A., Vigan, A., Galicher, R., et al. 2016, A&A, 587, A57