# A Clear View of a Cloudy Brown Dwarf Companion from High-resolution Spectroscopy

Jerry Xuan, Jason Wang, Jean-Baptiste Ruffio, Heather Knutson, Dimitri Mawet, Paul Mollière, Jared Kolecki, Arthur Vigan, Sagnick Mukherjee, Nicole Wallack, et al.

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

► To cite this version:

Jerry Xuan, Jason Wang, Jean-Baptiste Ruffio, Heather Knutson, Dimitri Mawet, et al. A Clear View of a Cloudy Brown Dwarf Companion from High-resolution Spectroscopy. The Astrophysical Journal, 2022, 937 (2), pp.54. 10.3847/1538-4357/ac8673. hal-03808127

HAL Id: hal-03808127

https://hal.science/hal-03808127

Submitted on 31 Mar 2023

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_2.png" /></figure>
<figure><img src="image_3.png" /></figure>

# A Clear View of a Cloudy Brown Dwarf Companion from High-resolution Spectroscopy

Jerry W. Xuan<fnref n="1" />, Jason Wang<fnref n="1" />,<fnref n="2" />,<fnref n="17" />, Jean-Baptiste Ruffio<fnref n="1" />, Heather Knutson<fnref n="3" />, Dimitri Mawet<fnref n="1" />,<fnref n="4" />, Paul Mollière<fnref n="5" />, Jared Kolecki<fnref n="6" />, Arthur Vigan<fnref n="7" />, Sagnick Mukherjee<fnref n="8" />, Nicole Wallack<fnref n="3" />, Ji Wang<fnref n="6" />, Ashley Baker<fnref n="1" />,<fnref n="17" />, Randall Bartos<fnref n="4" />, Geoffrey A. Blake<fnref n="3" />, Charlotte Z. Bond<fnref n="9" />, Marta Bryan<fnref n="10" />,<fnref n="17" />, Benjamin Calvin<fnref n="1" />,<fnref n="11" />, Sylvain Cetre¹², Mark Chun<fnref n="13" />, Jacques-Robert Delorme<fnref n="1" />,¹², Greg Doppmann¹², Daniel Echeverri<fnref n="1" />, Luke Finnerty<fnref n="11" />, Michael P. Fitzgerald<fnref n="11" />, Katelyn Horstman<fnref n="1" />, Julie Inglis<fnref n="3" />, Nemanja Jovanovic<fnref n="1" />, Ronald López<fnref n="11" />, Emily C. Martin<fnref n="14" />, Evan Morris<fnref n="14" />, Jacklyn Pezzato<fnref n="1" />, Sam Ragland¹², Bin Ren<fnref n="1" />, Garreth Ruane<fnref n="1" />,<fnref n="4" />, Ben Sappey<fnref n="15" />, Tobias Schofield<fnref n="1" />, Andrew Skemer<fnref n="14" />, Taylor Venenciano<fnref n="16" />, J. Kent Wallace<fnref n="4" />, and Peter Wizinowich¹²

## Abstract

Direct imaging studies have mainly used low-resolution spectroscopy ($R \sim 20–100$) to study the atmospheres of giant exoplanets and brown dwarf companions, but the presence of clouds has often led to degeneracies in the retrieved atmospheric abundances (e.g., carbon-to-oxygen ratio, metallicity). This precludes clear insights into the formation mechanisms of these companions. The Keck Planet Imager and Characterizer (KPIC) uses adaptive optics and single-mode fibers to transport light into NIRSPEC ($R \sim 35,000$ in the K band), and aims to address these challenges with high-resolution spectroscopy. Using an atmospheric retrieval framework based on petitRADTRANS, we analyze the KPIC high-resolution spectrum (2.29–2.49 µm) and the archival low-resolution spectrum (1–2.2 µm) of the benchmark brown dwarf HD 4747 B ($m = 67.2 \pm 1.8 M_{Jup}$, $a = 10.0 \pm 0.2$ au, $T_{eff} \approx 1400$ K). We find that our measured C/O and metallicity for the companion from the KPIC high-resolution spectrum agree with those of its host star within $1\sigma–2\sigma$. The retrieved parameters from the K-band high-resolution spectrum are also independent of our choice of cloud model. In contrast, the retrieved parameters from the low-resolution spectrum are highly sensitive to our chosen cloud model. Finally, we detect CO, H$_2$O, and CH$_4$ (volume-mixing ratio of log(CH$_4$) = −4.82 ± 0.23) in this L/T transition companion with the KPIC data. The relative molecular abundances allow us to constrain the degree of chemical disequilibrium in the atmosphere of HD 4747 B, and infer a vertical diffusion coefficient that is at the upper limit predicted from mixing length theory.

*Unified Astronomy Thesaurus concepts: Brown dwarfs (185); Atmospheric composition (2120); High resolution spectroscopy (2096); High angular resolution (2167); Atmospheric clouds (2180); Exoplanet atmospheric composition (2021)*

## 1. Introduction

The Keck Planet Imager and Characterizer (KPIC) is a new suite of instrument upgrades at Keck II, including a single-mode fiber injection unit (FIU; Delorme et al. 2021; Mawet et al. 2017) that feeds light into the upgraded NIRSPEC

(Martin et al. 2018; López et al. 2020), enabling high-resolution spectroscopy (HRS<fnref n="18" />) at $R \sim 35,000$ in the K band. By using single-mode fibers to inject light from planets and brown dwarfs at high contrast, KPIC provides suppression of the stellar point-spread function (PSF) at the fiber input and a stable line-spread function (LSF) that is independent of incoming wave-front aberrations (Mawet et al. 2017; Wang et al. 2021b). By observing at high resolution, we can further distinguish between star and planet light from their spectral

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

<footnote n="1">¹Department of Astronomy, California Institute of Technology, Pasadena, CA 91125, USA; jxuan@astro.caltech.edu</footnote>

<footnote n="2">²Center for Interdisciplinary Exploration and Research in Astrophysics (CIERA) and Department of Physics and Astronomy, Northwestern University, Evanston, IL 60208, USA</footnote>

<footnote n="3">³Division of Geological & Planetary Sciences, California Institute of Technology, Pasadena, CA 91125, USA</footnote>

<footnote n="4">⁴Jet Propulsion Laboratory, California Institute of Technology, 4800 Oak Grove Drive, Pasadena, CA 91109, USA</footnote>

<footnote n="5">⁵Max-Planck-Institut für Astronomie, Königstuhl 17, D-69117 Heidelberg, Germany</footnote>

<footnote n="6">⁶Department of Astronomy, The Ohio State University, 100 W 18th Avenue, Columbus, OH 43210 USA</footnote>

<footnote n="7">⁷Aix Marseille Univ, CNRS, CNES, LAM, Marseille, France</footnote>

<footnote n="8">⁸Department of Astronomy and Astrophysics, University of California, Santa Cruz, CA 95064, USA</footnote>

<footnote n="9">⁹UK Astronomy Technology Centre, Royal Observatory, Edinburgh EH9 3HJ, UK</footnote>

<footnote n="10">¹⁰Department of Astronomy, University of California at Berkeley, CA 94720, USA</footnote>

<footnote n="11">¹¹Department of Physics & Astronomy, 430 Portola Plaza, University of California, Los Angeles, CA 90095, USA
¹²W.M. Keck Observatory, 65-1120 Mamalahoa Highway, Kamuela, HI, USA</footnote>

<footnote n="13">¹³Institute for Astronomy, University of Hawaii at Hilo, 640 N Aohoku Place, Hilo, HI 96720, USA</footnote>

<footnote n="14">¹⁴Department of Astronomy & Astrophysics, University of California, Santa Cruz, CA 95064, USA</footnote>

<footnote n="15">¹⁵Center for Astrophysics and Space Sciences, University of California, San Diego, La Jolla, CA 92093, USA</footnote>

<footnote n="16">¹⁶Physics and Astronomy Department, Pomona College, 333 N. College Way, Claremont, CA 91711, USA</footnote>

<footnote>Received 2022 May 18; revised 2022 July 19; accepted 2022 August 1; published 2022 September 27</footnote>

<footnote n="17">¹⁷51 Pegasi b Fellow.</footnote>

<footnote n="18">¹⁸We will use HRS to abbreviate both high-resolution spectroscopy (the technique) and high-resolution spectrum (the data) in this paper. The same is true for LRS, low-resolution spectrum or spectroscopy.</footnote>
differences (Mawet et al. 2017; Wang et al. 2017). Recently, Wang et al. (2021b) published the KPIC detections of HR 8799 c, d, e, demonstrating the ability of KPIC to detect molecular lines and measure the rotational line broadening of planets at high contrast ($\Delta K \approx 11$) and small separations ($\approx 0.4''$) from their host star.

The atmospheric composition of a substellar object holds a wealth of information about its formation, accretion, and evolutionary history, as well as fundamental physical processes that shape its atmosphere. It is therefore important to assess how well KPIC and other fiber-fed, high-resolution spectrographs (e.g., Subaru/REACH; Kotani et al. 2020; and Very Large Telescope/HiRISE; Otten et al. 2021) can measure the atmospheric compositions of directly imaged planets and brown dwarfs. Specifically, previous studies of gas giant planet formation have highlighted the carbon-to-oxygen ratio (C/O) and metallicity (e.g., [C/H]) of the atmosphere as informative probes of formation history (e.g., Öberg et al. 2011; Madhusudhan 2012; Piso et al. 2015). To first order, a companion with a C/O and metallicity similar to that of its host star is consistent with formation via gravitational instability in a protostellar disk or fragmentation in a molecular cloud, akin to binary star formation (Bate et al. 2002). On the other hand, differences between the companion and stellar C/O are suggestive of core accretion (Pollack et al. 1996) as the likely formation mechanism, and in that scenario, could be used to constrain where the companion formed in the disk relative to ice lines of major C- and O-bearing molecules (e.g., H₂O, CO₂, and CO). This picture can be complicated by a variety of effects such as the relative amount of solids incorporated into the planet's atmosphere (e.g., Madhusudhan et al. 2014; Mordasini et al. 2016; Öberg & Bergin 2016; Gravity Collaboration et al. 2020; Pelletier et al. 2021).

So far, atmospheric characterization of directly imaged companions has mostly relied on low-resolution spectroscopy (LRS) with resolving powers of $R \approx 20–100$. LRS is sensitive to continuum emission originating from the deepest observable layer of the atmosphere and modified by opacity sources farther up. Many of these companions have temperatures warm enough for silicate clouds to condense in their atmospheres (Marley & Robinson 2015), and there is much evidence that cloud opacity plays an important role in the spectrum of directly imaged companions and brown dwarfs with L or L/T transition spectral types (e.g., Skemer et al. 2014; Burningham et al. 2017; Gravity Collaboration et al. 2020). However, due to our limited knowledge of cloud physics, a reliable assessment of atmospheric abundances from LRS could be fraught with degeneracies between clouds, the pressure-temperature (P-T) profile, and chemical abundances (e.g., Burningham et al. 2017). In addition, the retrieval results can also be highly sensitive to systematics in different data sets that are combined to obtain a wider wavelength coverage (Wang et al. 2020a). More encouragingly, Mollière et al. (2020) report atmospheric abundances that are relatively robust to clouds and model choices, though Burningham et al. (2021) show that issues such as an unphysically small radius could persist despite improvements in cloud modeling and extensive wavelength coverage (1–15 µm).

Recently, Wang et al. (2022) presented the first atmospheric free retrievals at high resolution for a directly imaged companion. They studied the L-type brown dwarf HR 7672 B ($T_{eff} \approx 1800$ K) using KPIC HRS and near-infrared

photometry, and measured carbon and oxygen abundances that are consistent within $<1.5\sigma$ to that of its host star. In this paper, we present a detailed atmospheric study of HD 4747 B using both KPIC HRS (K band) and archival low-resolution spectra (LRS) from 1–2.2 µm that we re-extract in a uniform manner. While the KPIC HRS resolves individual molecular lines and contains direct information about a companion's atmospheric abundances, LRS provides spectral shape and luminosity measurements, which has the potential to complement the HRS.

Compared to HR 7672 B, HD 4747 B is a colder L/T transition object ($T_{eff} \approx 1400$ K) with strong evidence for clouds and a similar color to directly imaged planets such as HR 8799 c, d, e (Crepp et al. 2018; Peretti et al. 2019). Like HR 7672 B, the wealth of prior knowledge available for HD 4747 B makes it a valuable benchmark object to test whether we can make robust inferences with spectroscopic data. First, we are able to precisely measure the dynamical mass of HD 4747 B (Section 2.2). Mass is a fundamental quantity that is poorly constrained for most directly imaged companions (Bowler 2016). Furthermore, given its high mass, HD 4747 B is expected to have formed via direct gravitational collapse in the same cloud or disk as its host star, which means that we can assume chemical homogeneity: the brown dwarf and primary star should share the same chemical composition. Finally, with the companion mass, observed luminosity, and stellar age, we can independently estimate the brown dwarf's radius from evolutionary models.

In this paper, we use the open-source radiative transfer code petitRADTRANS (Mollière et al. 2019, 2020) to fit the HRS and LRS for HD 4747 B in a “retrieval” framework. The main goals of our study are to measure the atmospheric composition of this brown dwarf companion using both the HRS and LRS, and to present a detailed characterization of its atmosphere, including constraints on clouds, chemical equilibrium or disequilibrium, and the detection of CH₄. In this process, we also explore the relative advantages and disadvantages of HRS versus LRS.

This paper is organized as follows: in Section 2, we summarize the system properties including our mass measurement for HD 4747 B. Our spectroscopic data and data reduction procedure is described in Section 3. We then discuss our spectral analysis framework in Section 4. We present individual and joint retrievals of the HRS and LRS in Sections 5, 6, and 7, respectively. We summarize the lessons learned in Section 8, and conclude in Section 9.

## 2. System Properties

### 2.1. Host Star

In this section, we summarize relevant properties of the host star. HD 4747 is a main-sequence, solar-type star located $\approx 19$ pc away based on its Gaia eDR3 parallax (Brown et al. 2021). Chromospheric emission in the Ca II H&K lines is visible in the stellar spectrum ($\log R_{HK} = -4.72 \pm 0.02$), which Peretti et al. (2019) used to derive an age of $2.3 \pm 1.4$ Gyr from the age-log $R_{HK}$ calibration of Mamajek & Hillenbrand (2008). This agrees with the gyro-chronological age estimate of $3.3_{-1.9}^{+2.3}$ Gyr from Crepp et al. (2018). These studies also converged on $T_{eff}$ around 5300–5400 K, and a surface gravity log(g) of 4.5–4.65. Of particular relevance to this study are the C/O ratio and metallicity of the host star, since we expect these to be roughly
similar to those of the brown dwarf. HD 4747 is found to have a sub-solar metallicity, with $[Fe/H] = -0.23 \pm 0.05$ from Peretti et al. (2019) and $[Fe/H] = -0.22 \pm 0.04$ from Crepp et al. (2018). Previous studies including Brewer et al. (2016) and Peretti et al. (2019) also measured the elemental abundances for the star, but either did not take into account non-local thermodynamic equilibrium (LTE) effects on their oxygen abundances (Amarsi et al. 2019) or do not quote error bars. We instead carry out a new analysis using the method described in Kolecki et al. (2021) to derive the abundances for different elements, and correct the results to account for 3D non-LTE effects (Amarsi et al. 2019) on the results. For this analysis, we used an archival spectrum from FEROS (Kaufer et al. 1997), which covers 350–920 nm at R = 48,000. Using this spectrum, we measure the equivalent widths of absorption lines and compare them to model stellar atmospheres in an iterative approach using the MOOG code (Sneden 1973). From our derived carbon and oxygen abundances, we find $C/O = 0.48 \pm 0.08$. The iron abundance is $[Fe/H] = -0.30 \pm 0.5$, in agreement with previous studies.

Since Fe condenses out for temperatures below ≈1800 K (Marley & Robinson 2015), it is not a relevant gaseous absorber in the photosphere of HD 4747 B. Therefore, the more useful metrics for comparison are C and O. From our analysis above, we find $[C/H] = -0.08 \pm 0.06$ and $[O/H] = -0.02 \pm 0.04$ for the host star. $[C/H]$ is defined as $\log_{10}(N_C/N_H)_{\text{star}} - \log_{10}(N_C/N_H)_{\text{sun}}$, where $N_C$ and $N_H$ are the number fraction of C and H, respectively. $[O/H]$ is defined similarly. We adopt Asplund et al. (2009) as our solar reference in order to be consistent with petitRADTRANS, which we use to model the atmosphere of HD 4747 B.

## 2.2. Orbit and Dynamical Mass

The orbit and mass of HD 4747 B have been measured by several studies using relative astrometry from Keck/NIRC2, host star radial velocities (RVs) from Keck/HIRES, and Gaia-Hipparcos absolute astrometry (Brandt et al. 2019; Xuan & Wyatt 2020). Here, we take advantage of 23 yr of RV observations published in Rosenthal et al. (2021) and the improved precision of the Hipparcos-Gaia Catalog of Accelerations (Brandt 2021) based on Gaia eDR3 (Brown et al. 2021) to update the orbit and mass of HD 4747 B. HD 4747 B shows significant proper motion anomalies (PMa) in both the Gaia and Hipparcos epochs, with signal-to-noise ratios (S/Ns) of 77.2 and 9.1, and the position angle and amplitude of the PMa is consistent with being induced by the brown dwarf companion. For the relative astrometry, we use data points tabulated in Brandt et al. (2019), except for the two GPI epochs measured by Crepp et al. (2018), which we replaced with our new measurements from Section 3.2. We choose not to use the companion RV as measured by KPIC for this fit, because it does not appreciably improve our already well-constrained orbital solution.

To fit the relative astrometry, radial velocity, and absolute astrometry from Gaia and Hipparcos together, we use the orvar package (Brandt et al. 2021b), which is designed to jointly fit these types of data and takes into account the Gaia and Hipparcos astrometry at the epoch astrometry level using htof (Brandt et al. 2021a). We use the priors listed in Table 4 of Brandt et al. (2021b) for the fitted parameters. The posterior is sampled using the parallel-tempering Markov Chain Monte Carlo sampler (Vousden et al. 2016), a fork of emcee

<table><caption>**Table 1**
Selected Parameters from Orbit Fit</caption><thead><tr><td>Parameter</td><td>Value</td></tr></thead><tbody><tr><td><i>M</i> (<i>M</i><sub>☉</sub>)</td><td>0.85 ± 0.04</td></tr><tr><td><i>m</i> (<i>M</i><sub>Jup</sub>)</td><td>67.2 ± 1.8</td></tr><tr><td><i>a</i> (au)</td><td>10.0 ± 0.2</td></tr><tr><td>inclination (deg)</td><td>48.0 ± 0.9</td></tr><tr><td>ascending node (deg)</td><td>89.4 ± 1.1</td></tr><tr><td>period (yr)</td><td>33.2 ± 0.4</td></tr><tr><td>argument of periastron (deg)</td><td>267.2 ± 0.5</td></tr><tr><td>eccentricity</td><td>0.7317 ± 0.0014</td></tr><tr><td>Epoch of periastron (JD)</td><td>2462615 ± 155</td></tr></tbody></table>

**Note.** The dynamical mass of the host star, which is fit as a free parameter, agrees well with isochrone-derived masses from Peretti et al. (2019) and Crepp et al. (2018).

(Foreman-Mackey et al. 2013) over 50,000 steps with 10 temperatures and 100 walkers. The fits converged as determined by visual inspection of the chains, and we discarded the first 10% as burn-in. In orvar, the system parallax and other linear parameters are analytically marginalized out to speed up the fits.

The resulting orbit and mass measurements are tabulated in Table 1, while the model fits are shown in Appendix A. We find a companion mass of $m = 67.2 \pm 1.8 M_{\text{Jup}}$, which is consistent with previous values, but more precise. We checked the orvar results with a second fit where we model the Gaia and Hipparcos astrometry using the methodology in Xuan & Wyatt (2020). This gives $m = 67.1 \pm 2.0 M_{\text{Jup}}$, consistent with the orvar result. Furthermore, the companion mass and orbital parameters we find are also consistent with results from orbit fits that only use RV and imaging data (no Gaia-Hipparcos astrometry) from Peretti et al. (2019) and Crepp et al. (2018). We adopt the companion mass from our orvar fit for the spectral analysis in this paper.

# 3. Spectroscopic Data

## 3.1. High-resolution Spectroscopy

### 3.1.1. KPIC Observations

We observed HD 4747 B on UT 2020 September 28 with Keck/NIRSPEC. The data were collected using the first version of the KPIC FIU (Delorme et al. 2021). The FIU is located downstream of the Keck II adaptive optics system and is used to inject light from a selected object into one of the single-mode fibers connected to NIRSPEC. We obtain spectrum in the K band, which is broken up into nine echelle orders from 1.94–2.49 µm. The observing strategy is identical to that of Wang et al. (2021b). In short, we placed the companion on the fiber with the highest throughput and acquired six exposures of 600 s each, for a total integration time of 1 hr. The relative astrometry of the companion was computed using whereistheplanet.com (Wang et al. 2021a), based on data in Peretti et al. (2019). For calibration purposes, we acquired a pair of 60 s exposures of the host star before observing the companion, and a pair of 60 s exposures of a telluric standard star (HIP 6960) after the companion exposures so as to share nearly the same airmass. Using exposures on the host star, we calculated an end-to-end throughput from the top of the atmosphere to the detector of 1.8%–2.0% during the observations.
<figure><img src="image_5.png" /><figcaption>**Figure 1.** Cross-correlation functions (CCFs) in blue show detections of CO and H₂O using three spectral orders of the KPIC HRS. Gray lines represent CCFs of the background flux (from the slit) and speckle flux, whose standard deviations are used as estimates of the CCF noise. The vertical dashed lines at 15 km s⁻¹ show the expected RV of the companion from its known orbit. The strong structure in the blue CCFs outside the peaks arise because we only fit single molecules here.</figcaption></figure>

### 3.1.2. Data Reduction

To extract the spectra from the raw data, we follow the procedure outlined in Wang et al. (2021b), which the KPIC team has implemented in a public Python pipeline.<fnref n="19" /> The images for all objects were reduced in the same way. First, we removed the thermal background from the images using combined instrument background frames taken during daytime. As shown in Wang et al. (2021b), the thermal background of our data is dominated by the warm optics rather than the sky background. We also remove persistent bad pixels identified using the background frames. Then, we use data from the telluric standard star to fit the trace of each column in the four fibers and nine spectral orders, which give us the position and standard deviation of the PSF (in spatial direction) at each column. The trace positions and widths were additionally smoothed using a cubic spline in order to mitigate random noise. We adopt the trace locations and widths as the LSF positions and widths in the dispersion dimension.

For every frame, we then extracted the 1D spectra in each column of each order. To remove residual background light, we subtracted the median of pixels that are at least 5 pixels away from every pixel in each column. Finally, we used optimal extraction to sum the flux using weights defined by the 1D Gaussian LSF profiles calculated from spectra of the telluric star.

The extracted spectra have a median S/N of ≈8 per pixel element, which has a typical width of 0.2 Å, and consists of a mixture of light from the brown dwarf companion and stellar speckles. The S/N of KPIC is optimized for wavelengths around 2.3 µm, where CO has a series of strong absorption lines. For our analysis, we use three spectral orders from 2.29–2.49 µm, which contain the strongest absorption lines from the companion and have relatively few telluric absorption lines. Note that the three spectral orders have gaps in between them, so we have data over ≈0.13 µm (instead of 0.2 µm; see Figure 3).

As a preliminary analysis, we cross-correlate our KPIC spectra with single-molecule templates assuming $T_{\text{eff}} = 1400$ K and $\log(g) = 5.5$ from the Sonora model grid (Marley et al. 2021). In short, we estimate the maximum likelihood value for both the single-molecule companion flux and speckle flux in the data as a function of RV (radial velocity) shift using the

method described in Wang et al. (2021b), which is based on Ruffio (2019). We find that H₂O and CO are detected with S/Ns of 8.5 and 13.5, respectively (Figure 1). CH₄ is not detected with statistical confidence in this crude analysis, but we present evidence for a weak CH₄ detection in Section 5.4.

## 3.2. Low-resolution Spectroscopy

### 3.2.1. Gemini Planet Imager IFS

The Gemini Planet Imager (GPI) observed HD 4747 B on UT 2015 December 24 and 25, in the *K*₁ (1.90–2.19 µm, *R* = 66) and *H* (1.50–1.80 µm, *R* = 46.5) bands, respectively, and the data were published in Crepp et al. (2018). After doing some fits to the published spectrum, we found that the average flux levels of the *K*₁ and *H* bands are inconsistent, and the error bars appear to be significantly overestimated.

We therefore re-extracted the GPI spectrum using the *pyKLIP* package (Wang et al. 2015), which models a stellar PSF with Karhunen-Loève Image Processing (KLIP; also known as principal component analysis) following the framework in Soummer et al. (2012) and Pueyo (2016). We tested various model choices to minimize the residuals after stellar PSF subtraction while preserving the companion signal. A key parameter we tuned was the number of Karhunen-Loève (KL) modes. KL modes represent an orthogonal basis for patterns in the images that are used to model the stellar PSF. We chose five and 12 KL modes to subtract the stellar PSF in the *H*- and *K*₁-band data, respectively. After subtracting the stellar PSF, we first extracted the companion's relative astrometry in terms of separation and position angle, which are reported in Appendix B (Table 5). Then, we extracted the flux at the companion's determined location as a function of wavelength, which gave us the raw spectrum. Note that rather than using spectral differential imaging (SDI) to subtract the stellar PSF, we only used angular differential imaging (ADI). For a bright companion like HD 4747 B, ADI is more than sufficient to properly remove the PSF of the star given sufficient parallactic angle rotation.

To flux-calibrate the raw spectrum, we used the satellite spot flux ratios<fnref n="20" /> to find the companion-to-star flux ratio. To obtain the observed flux density of the companion, we empirically

<footnote n="19">¹⁹ https://github.com/kpicteam/kpic_pipeine</footnote>

<footnote n="20">²⁰ https://www.gemini.edu/instrumentation/gpi/capability</footnote>
<figure><img src="image_6.png" /><figcaption>**Figure 2.** Our extracted LRS of HD 4747 B from GPI and SPHERE are plotted in black circles and squares, respectively, and the published spectra for GPI is shown in gray dots from Crepp et al. (2018). Overplotted in color are spectra of a L9.5 dwarf and T0 dwarf from SPEX (Burgasser 2014), which show good agreement with our extracted spectra, demonstrating that HD 4747 B is consistent with a spectral type near the L/T transition.</figcaption></figure>

determined the flux scaling factor $R^2/d^2$ by fitting a PHOENIX model (Husser et al. 2013) of the star ($T_{eff} = 5400$, log(g) = 4.5, and $[Fe/H] = -0.5$) using the star's Two Micron All Sky Survey J-, H-, and K-band magnitudes (Cutri et al. 2003) and the Gaia G-band magnitudes (Riello et al. 2021). The zero-point fluxes and filter transmission of the photometric bands are downloaded from the Spanish Virtual Observatory (SVO) Filter Service<fnref n="21" /> and the Gaia website.<fnref n="22" /> To obtain measurement uncertainties, we injected 20 fake companions at the same separation and equally spaced position angles in the data, and repeated the same spectral extraction process. We avoided using the fake injections that were within 20° of the real companion to avoid biasing the fluxes. We inflated the uncertainties on the extracted spectra by 2.5% to account for errors in the stellar flux calibration. The value of 2.5% is estimated by comparing our empirically computed flux scaling factor with the value of $R^2/d^2$ of the star (using the radius from Peretti et al. 2019 and the Gaia parallax).

### 3.2.2. SPHERE IFS

HD 4747 B was observed on UT 2016 December 12 and 2017 September 28 with the Spectro-Polarimetric High-contrast Exoplanet Research (SPHERE; Beuzit et al. 2019). The SPHERE Integral Field Spectrograph (IFS; Claudi et al. 2008) collects data in the YH band from 0.95–1.6 µm ($R = 29$). The extracted spectrum was published in Peretti et al. (2019), but is not available. We therefore reduced the raw data using the SPHERE pipeline (Vigan 2020), and performed a similar post-processing procedure with pyKLIP as described above for the GPI spectra. The only difference is that we needed to use ADI+SDI to perform PSF subtraction for the SPHERE IFS data, which did not have enough parallactic angle rotation (only ≈0°/2). For the SPHERE IFS data, flux calibration is based on unocculted observations of the host star. We chose to use the 2017 data for our analysis since it was taken under much better

observing conditions and yields a slightly higher spectral S/N than the 2016 data, despite shorter integration times.

### 3.2.3. Results and Comparison with Previous LRS

Our newly extracted GPI and SPHERE spectra are plotted in black circles and squares, respectively, in Figure 2, and available in Appendix B (Table 4). The absolute flux scaling of our GPI spectrum agrees well with the published spectrum in gray from Crepp et al. (2018), but the uncertainties are much smaller. The shape of our new SPHERE spectrum also agrees well with that in Peretti et al. (2019). From the PSF-subtracted images, the brown dwarf companion is detected with a median S/N per wavelength bin of ≈61 and ≈26 in the GPI H and K₁ bands, and ≈20 in the SPHERE data. When comparing the extracted spectrum to that of field brown dwarfs from the SPEX library (Burgasser 2014) in Figure 2, we find that the newly extracted GPI spectrum is in better agreement compared to the previously published spectrum. As in Crepp et al. (2018) and Peretti et al. (2019), we find a spectral type near the L/T transition (the best matching spectra were from an L9.5 and T0 dwarf). The SPHERE IFS spectrum increases our wavelength coverage by a factor of ~2, which we find is important for constraining model atmosphere parameters in our fits to the LRS.

## 4. Spectral Analysis

### 4.1. Forward Modeling the KPIC High-resolution Spectrum

Here, we briefly describe the framework to forward model and fit the HRS from KPIC, which follows Wang et al. (2021b). When a companion of interest is aligned with one of the KPIC fibers, the companion light and a fraction of light from the host star's speckle field are injected into the fiber. At the projected separation of HD 4747 B (≈0.6), we find the speckles are roughly the same brightness as our companion ($K_s \approx 14.4$ from Crepp et al. 2016). Furthermore, the light is transmitted through Earth's atmosphere and modulated by the instrument optics. Thus, we build two forward models (one for the companion, one for speckles) and jointly fit them as a linear

<footnote n="21">²¹ http://svo2.cab.inta-csic.es/theory/fps/</footnote>

<footnote n="22">²² https://www.cosmos.esa.int/web/gaia/edr3-passbands</footnote>
combination. Below we detail how we generate each of the model components.

The companion spectral templates are generated with petitRADTRANS. We shift the templates in wavelength space to fit for the radial velocity. Then, we rotationally broaden the templates by a projected rotation rate $v \sin i$ using the fastRotBroad function in PyAstronomy (Czesla et al. 2019), and convolve the templates with the instrumental LSF. The effect of limb darkening is included in petitRADTRANS by integrating intensities along multiple angles between the ray and atmospheric normal.

Next, we multiply the companion model by the telluric response function, which characterizes the atmospheric transmission as a function of wavelength and includes telluric absorption lines. The telluric model is calculated by dividing the spectrum of the standard star (HIP 6960) by a PHOENIX stellar model with matching properties ($T_{\text{eff}} = 9200$ and $\log(g) = 4.0$).

To model the speckle contribution to the data, we use on-axis observations of the host star taken before the companion exposures. These observations are reduced in the same way as the companion spectra, but have a much higher S/N. Unlike the companion models, the host star observations are already modulated by telluric transmission.

The last step is to remove continuum variations. The KPIC spectra are not flux-calibrated and contain a smoothly varying continuum due to stellar speckles and wavelength-dependent atmospheric refraction. Therefore, we apply high-pass filtering with a median filter of 100 pixels (~0.002 µm) on both the data and models to subtract the continuum following Wang et al. (2021b). To determine the optimal filter size, we carried out a series of injection-recovery tests, and found that ~100 pixels is best at recovering weak molecular signals for our data set. Larger filter sizes (e.g., 200 pixels or more) do not remove the continuum sufficiently, and smaller filter sizes (50 pixels or less) tend to be overly aggressive at removing weak molecular signals.

Finally, we flux-normalize both the companion and stellar models and multiply them by different flux scaling factors, which are fitted parameters. The flux scales are in units of counts as measured by the NIRSPEC detector. After scaling, the companion and speckle models are added, and the same high-pass filter is applied on the final model before fitting it to the data.

## 4.2. Atmospheric Retrieval Setup

We implement a “retrieval” framework based on petitRADTRANS to model the data, which means that we freely retrieve the chemical abundances, vertical temperature structure, and cloud properties from the data. Previous studies have used retrievals to model HRS of self-luminous exoplanets and brown dwarfs (e.g., Burningham et al. 2017; Mollière et al. 2020), and show that it can be a powerful alternative to fitting self-consistent grid models, which solve for the abundances and temperature profiles with physical assumptions such as chemical equilibrium. The retrieval approach allows more flexibility to fit the data and can potentially provide much more detailed information about the atmospheric properties, with the caveat that it is important to check for physical plausibility of the models since retrievals need not be self-consistent.

In our main set of retrievals, we fit for the chemical abundances in terms of C/O and atmospheric metallicity [C/H]$^{23}$ along with a

quench pressure (where the chemical timescale of a certain reaction is equal to the mixing timescale) to allow for disequilibrium chemistry, the temperature profile (Section 4.2.1), the cloud structure (Section 4.2.2), and other parameters such as the radius. We denote these quenched chemistry retrievals to distinguish from free retrievals where the abundances of each gas species is fit independently. Each component of the model is described in the subsections below. We use the correlated-k and line-by-line opacity sampling methods in petitRADTRANS for the low-resolution and high-resolution retrievals, respectively. For high-resolution, we include opacities for CO, $H_2O$, $CH_4$, $NH_3$, and $CO_2$, and for low-resolution, we additionally include Na and K. This is because the alkali lines have wings, which affect the ~1 µm portion of the LRS, while their opacities are negligible over the portion of the K band covered by our HRS. We repeated our baseline HRS retrieval with Na and K included and found that the addition of these two species did not influence the results or improve the fit.

Because the native high-resolution opacities are at $R = 10^6$, much higher than the resolution of our HRS resolution ($R \approx 35,000$), we down-sampled the opacity tables by a factor of six in order to speed up the retrievals (by roughly the same factor) and reduce the corresponding computational cost. We checked that the maximum deviation in synthetic spectra obtained by using the down-sampled opacities relative to the full-resolution opacity model is <5% of the minimum HRS error bars. In addition, we repeated our fiducial HRS retrieval with the native opacities ($R = 10^6$) and found that it yielded the same results. We re-binned the correlated-k opacities to $R = 200$ for our fits to the LRS, which has a maximum resolution of 66. We also repeated our fiducial LRS retrieval at the native $R = 1000$ opacities and found the results are fully consistent.

### 4.2.1. Temperature Structure and Chemistry

We retrieve the P–T profile of the brown dwarf between $P = 10^{-4}–10^3$ bars, which sets the vertical extent of the atmosphere. We use the P–T profile parameterization from Mollière et al. (2020), which has six free parameters. The spatial coordinate is an optical depth $\tau = \delta P^{\alpha}$, where $\delta$ and $\alpha$ are the first two parameters. The atmosphere then consists of a high-altitude region (top of atmosphere to $\tau = 0.1$) fitted with three temperature points equidistant in log pressure, a middle radiative region ($\tau = 0.1$ to radiative-convective boundary) that uses the Eddington approximation with $T_0$ as the “internal temperature,” and a lower region (radiative-convective boundary to bottom of atmosphere), which is set to follow the moist adiabatic temperature gradient once the atmosphere becomes unstable to convection (Mollière et al. 2020). We ignore stellar irradiation as a source of heat because the total incident energy on HD 4747 B at periastron (≈2.7 au) is approximately four orders of magnitude less than its luminosity, which is dominated by the brown dwarf’s internal energy.

In our quenched chemistry retrievals, the C/O, [C/H], and P–T profile determine the equilibrium chemical abundances (mass fractions of molecules) as a function of pressure, by interpolating the chemical equilibrium table from Mollière et al. (2020). The opacities we include in the models are listed in Section 4.2. In petitRADTRANS, the abundances of all metals except oxygen are assumed to scale together such that $[C/H] = [\text{Si/H}] = [\text{N/H}]$, etc. Then, C/O and [C/H] are combined to set the oxygen abundance (Mollière et al. 2020).

<footnote>$^{23}$ We denote the atmospheric metallicity as [C/H] because we are only sensitive to C- and O-bearing molecules in this brown dwarf’s atmosphere.</footnote>
<table><caption>**Table 2**

Spectral Retrievals Carried Out on HD 4747 B</caption><thead><tr><th>Data/Reference</th><th>Cloud Model</th><th>C/O</th><th>[C/H]</th><th>Radius (R<sub>Jup</sub>)</th><th>log(g)</th><th>T<sub>eff</sub> (K)</th><th>B</th></tr></thead><tbody><tr><td><strong>HRS (KPIC)</strong></td><td><strong>EddySed (MgSiO<sub>3</sub>, am)</strong></td><td><strong>0.66 ± 0.04</strong></td><td><strong>-0.10<sup>+0.18</sup><sub>-0.15</sub></strong></td><td><strong>0.82<sup>+0.19</sup><sub>-0.13</sub></strong></td><td><strong>5.39<sup>+0.15</sup><sub>-0.18</sub></strong></td><td><strong>1652<sup>+128</sup><sub>-218</sub></strong></td><td>1.0</td></tr><tr><td>HRS</td><td>EddySed (MgSiO<sub>3</sub>, cd)</td><td>0.67 ± 0.04</td><td>-0.06<sup>+0.23</sup><sub>-0.18</sub></td><td>0.90 ± 0.19</td><td>5.32<sup>+0.20</sup><sub>-0.17</sub></td><td>1577<sup>+167</sup><sub>-253</sub></td><td>1.15</td></tr><tr><td>HRS</td><td>Clear</td><td>0.67<sup>+0.05</sup><sub>-0.04</sub></td><td>-0.09<sup>+0.24</sup><sub>-0.16</sub></td><td>0.87<sup>+0.19</sup><sub>-0.17</sub></td><td>5.34<sup>+0.19</sup><sub>-0.17</sub></td><td>1677<sup>+132</sup><sub>-142</sub></td><td>0.61</td></tr><tr><td>HRS</td><td>Clear (chemical equilibrium)</td><td>0.60 ± 0.02</td><td>0.73<sup>+0.40</sup><sub>-0.31</sub></td><td>0.69<sup>+0.12</sup><sub>-0.06</sub></td><td>5.27<sup>+0.20</sup><sub>-0.14</sub></td><td>1402<sup>+143</sup><sub>-110</sub></td><td>1.6 × 10<sup>-3</sup></td></tr><tr><td><strong>LRS (GPI+SPHERE)</strong></td><td><strong>EddySed (MgSiO<sub>3</sub>, am)</strong></td><td>0.55<sup>+0.06</sup><sub>-0.14</sub></td><td>0.22<sup>+0.25</sup><sub>-0.47</sub></td><td>0.70<sup>+0.05</sup><sub>-0.03</sub></td><td>5.53<sup>+0.04</sup><sub>-0.05</sub></td><td>1473<sup>+17</sup><sub>-20</sub></td><td>1.0</td></tr><tr><td>LRS</td><td>EddySed (MgSiO<sub>3</sub>, cd)</td><td>0.45<sup>+0.08</sup><sub>-0.09</sub></td><td>-0.27<sup>+0.17</sup><sub>-0.19</sub></td><td>0.77 ± 0.04</td><td>5.45<sup>+0.04</sup><sub>-0.05</sub></td><td>1443 ± 28</td><td>0.69</td></tr><tr><td>LRS</td><td>EddySed (MgSiO<sub>3</sub> + Fe, am)</td><td>0.66<sup>+0.07</sup><sub>-0.10</sub></td><td>0.21<sup>+0.18</sup><sub>-0.24</sub></td><td>0.73 ± 0.03</td><td>5.50<sup>+0.03</sup><sub>-0.04</sub></td><td>1458<sup>+21</sup><sub>-19</sub></td><td>1.54</td></tr><tr><td>LRS</td><td>EddySed (MgSiO<sub>3</sub> + Fe, cd)</td><td>0.29<sup>+0.06</sup><sub>-0.07</sub></td><td>-0.51<sup>+0.17</sup><sub>-0.19</sub></td><td>0.75 ± 0.03</td><td>5.47<sup>+0.04</sup><sub>-0.03</sub></td><td>1453<sup>+24</sup><sub>-21</sub></td><td>2.65</td></tr><tr><td>LRS</td><td>Clear</td><td>0.12<sup>+0.02</sup><sub>-0.01</sub></td><td>-1.37<sup>+0.07</sup><sub>-0.05</sub></td><td>1.10 ± 0.04</td><td>5.12 ± 0.03</td><td>1262 ± 16</td><td>7.0 × 10<sup>-26</sup></td></tr><tr><td>Peretti et al. (2019)</td><td>Cloudy retrieval<sup>a</sup></td><td>0.13<sup>+0.14</sup><sub>-0.08</sub></td><td>-1.15<sup>+0.47</sup><sub>-0.39</sub></td><td>0.85 ± 0.03<sup>b</sup></td><td>5.40 ± 0.03</td><td>1350 ± 50</td><td>...</td></tr><tr><td>Crepp et al. (2018)</td><td>Cloudy grid<sup>b</sup></td><td>...</td><td>...</td><td>...</td><td>5.2<sup>+0.5</sup><sub>-0.6</sub></td><td>1410<sup>+130</sup><sub>-140</sub></td><td>...</td></tr></tbody></table>

Notes. The rightmost column lists the Bayes factor (B) for each retrieval, with the EddySed (MgSiO<sub>3</sub>, am) model as the baseline model with B = 1 (see Section 4.3.1 for an explanation of model comparison with B). We adopt the first row (in bold) as our final results for this paper. A few key parameters and their central 68% credible interval with equal probability above and below the median are listed, along with values for common parameters from two previous studies. For our cloudy models, “am” stands for amorphous cloud particles + Mie scattering, and “cd” stands for crystalline particles + DHS model, as described in Section 4.2.2. Except for the HRS model labeled chemical equilibrium, all of our other models are quenched chemistry retrievals.

<sup>a</sup> Peretti et al. (2019) carried out cloudy retrievals on their SPHERE spectrum (1.0–1.65 µm) and archival K and L photometry with the HELIOS-R code (Lavie et al. 2017), while Crepp et al. (2018) fitted their extracted GPI spectrum (1.5–2.2 µm) to a cloudy grid model (Saumon & Marley 2008).

<sup>b</sup> Peretti et al. (2019) placed a Gaussian prior of 1.0 ± 0.1R<sub>Jup</sub> on the radius.

We use Asplund et al. (2009) as our reference for the solar metallicity in these models.

Finally, we include a quench pressure $P_{\text{quench}}$, which fixes the abundances of H<sub>2</sub>O, CO, and CH<sub>4</sub> where $P < P_{\text{quench}}$ using the equilibrium values found at $P_{\text{quench}}$ (Zahnle & Marley 2014; Mollière et al. 2020). The inclusion of $P_{\text{quench}}$ allows for the possibility of disequilibrium chemistry, which occurs where the atmospheric mixing timescale is shorter than the chemical reaction timescale. We only include a quench pressure for the net reaction between H<sub>2</sub>O, CO, and CH<sub>4</sub> because these molecules are the only ones detectable in our KPIC HRS (see Section 5.4 for the CH<sub>4</sub> detection), and chemical kinetics modeling indicates that the abundances of these three molecules are closely linked to each other by a series of reactions (e.g., Moses et al. 2013). In summary, our quenched chemistry retrievals use C/O, [C/H], and $P_{\text{quench}}$ to set the abundances of each gas species for a given P–T profile.

### 4.2.2. Clouds

Crepp et al. (2018) and Peretti et al. (2019) analyzed LRS for HD 4747 B and found evidence for a cloudy atmosphere. We summarize their results in Table 2 along with our new measurements. In this study, we consider both clear and cloudy models in order to explore the sensitivity of our retrieved abundances to the assumed cloud properties. For our cloudy model, we use the EddySed model from Ackerman & Marley (2001) as implemented in petitRADTRANS (Mollière et al. 2020). In this model, the cloud particles both absorb and scatter the outgoing photons from the atmosphere according to measured optical properties (Mollière et al. 2019). The cloud particles can be either crystalline or amorphous, and the opacities of the clouds are computed assuming either homogeneous and spherical particles, modeled with Mie theory, or irregularly shaped cloud particles, modeled with the distribution of hollow spheres (DHS; Min et al. 2005; Mollière et al. 2019).

For HD 4747 B, we consider models with two different cloud species (MgSiO<sub>3</sub> and Fe) and properties (amorphous or crystalline particles). We choose to focus on MgSiO<sub>3</sub> and Fe for several reasons. First, the condensation curves of these two species intersect the thermal profile of a $T_{\text{eff}}$ = 1400 K, $log(g)$ = 5.5 object from the Sonora atmospheric model (Marley et al. 2021) at ~10–50 bars. While the Sonora model is cloudless, it provides a rough estimate of which cloud species are relevant. Second, recent theoretical work has shown that MgSiO<sub>3</sub> is expected to be the most important cloud species for substellar objects with $T_{\text{eff}}$ > 950 K due to its low nucleation energy barriers and the relatively high elemental abundances of Mg, Si, and O (Gao et al. 2020). Finally, studies using mid-IR spectroscopy from Spitzer have found direct evidence for a MgSiO<sub>3</sub> absorption feature at ~10 µm in field brown dwarfs (Cushing et al. 2006; Luna & Morley 2021), and specifically amorphous MgSiO<sub>3</sub> (Burningham et al. 2021). Although MgSiO<sub>3</sub> and Fe clouds do not have distinct features in the near-IR, they still impact the near-IR spectrum by contributing a wavelength-dependent opacity. Our baseline model uses amorphous MgSiO<sub>3</sub> modeled with Mie theory (abbreviated MgSiO<sub>3</sub>, “am”) for the clouds. In addition, we also consider models with MgSiO<sub>3</sub>, “cd,” which assumes crystalline cloud particles modeled with DHS, as well as models with two cloud species (MgSiO<sub>3</sub> + Fe) for the LRS.

### 4.2.3. Methane Opacities

Given that HD 4747 B is located near the L/T transition for brown dwarfs, we might expect to observe methane in its atmosphere. Previous L-band studies have detected methane in field brown dwarfs with spectral types as early as mid-L, or up to $T_{\text{eff}}$ ≈ 1800 K (Noll et al. 2000; Johnston et al. 2019). In this study, we adopted the HITEMP CH<sub>4</sub> line list from Hargreaves et al. (2020), which we convert into opacities following the petitRADTRANS documentation. When cross-correlating a model generated with the HITEMP CH<sub>4</sub> opacities with a late T dwarf, we obtained a CCF S/N of ≈15, in comparison to ≈5
when cross-correlating with a model generated from the default CH₄ opacities from ExoMol (Yurchenko & Tennyson 2014) in petitRADTRANS.

### 4.2.4. Additional Fit Parameters

petitRADTRANS computes the flux density as emitted at the surface of the object. For the LRS, we scale the model by the distance and companion radius, where the radius is another free parameter, and the distance is taken from the Gaia eDR3 parallax (Brown et al. 2021). For the HRS, we also fit the companion's radial velocity and v sin i, as well as an error multiple term to account for any underestimation in the data uncertainties.

Due to imperfect starlight subtraction in the spectral extraction process, we found that our LRS likely still contains correlated noise from the wavelength-dependence of speckles, as has been noted by several previous studies on high-contrast companions (e.g., Wang et al. 2020b; De Rosa et al. 2016; Samland et al. 2017; Currie et al. 2018; Wang et al. 2021c). This is evident in the residual frames, where we can see speckles at 5%–20% of the companion intensity in the PSF-subtracted images. We therefore adopt a Gaussian process with a squared exponential kernel to empirically estimate the correlated noise in the GPI H, K, and SPHERE YJH bands when fitting models to the data. Following Wang et al. (2020b), we assume that our extracted error bars contain a fraction $f_{\text{amp}}$ of correlated noise, and $1 - f_{\text{amp}}$ of white noise, and fit for $f_{\text{amp}}$ and the scale of correlation l. This adds $2 \times 3 = 6$ additional parameters to the retrievals.

As an alternative model, we also tried fitting the LRS with error inflation terms and flux scaling factors for the SPHERE and GPI spectra along the lines of Mollière et al. (2020), but found that our results were very sensitive to our choice of prior for the flux scaling factor. We conclude that our GP model is better suited to account for correlated noise from speckles, and use it in all LRS fits presented in this work.

### 4.3. Priors

We adopt uniform or log-uniform priors for all model parameters except for the mass, for which we use a Gaussian prior of $67.2 \pm 1.8 M_{\text{Jup}}$ from the dynamical mass measurement (Section 2.2). For the parametric P-T profile parameters, we exclude profiles that contain temperature inversions, as the heat budgets of widely separated companions are dominated by their internal luminosities. For the companion's radius, we use a uniform prior between 0.6 and 1.2 $R_{\text{Jup}}$. When including a quench pressure, we use a log-uniform prior from $10^{-4}$ to $10^3$, which is the full pressure range of our models. The priors for all retrieval parameters are tabulated in Table 6 of Appendix C.

#### 4.3.1. Model Fitting with Nested Sampling

We use nested sampling as implemented by *dynesty* (Speagle 2020) to find the posterior distributions for the model parameters. Specifically, we use 200 live points and adopt the stopping criterion that the estimated contribution of the remaining prior volume to the total evidence is less than 1%. We repeated a few retrievals using 1000 live points and found the evidence remains roughly the same, implying the fits have converged when using 200 live points.

One advantage of adopting nested sampling is that we can use the Bayesian evidence from each fit to calculate the Bayes

factor B, which assesses the relative probability of model $M_2$ compared to $M_1$. We will use the Bayes factor to compare different models throughout this paper to determine whether a given $M_2$ is justified over $M_1$. In Table 2, we take a baseline model (MgSiO₃, am) to be $M_1$ and compare other models to it. Based on Jeffreys (1983), a model with 100 times lower B than the model with the highest B can be “decisively” rejected. B of $\lesssim 10$ is considered weak evidence for preferring one model over the other. We first run retrievals with only the HRS (Section 5), only the LRS (Section 6), as well as joint retrievals with both HRS and LRS (Section 7).

## 5. High-resolution Retrievals (KPIC)

### 5.1. Overview

From our HRS retrievals of HD 4747 B, we find that both clear and cloudy models yield consistent results for the atmospheric parameters (abundances, temperature structure, and quenching) and bulk properties (radius, radial velocity, and spin). A few selected parameters are plotted in Figure 4 and tabulated in Table 2. The insensitivity of the HRS retrieval results to clouds, a major finding of this paper, is discussed in Section 5.2. In Figure 3, we plot the data, a best-fit model, and residuals for the baseline HRS retrieval. We report values from this retrieval as the final results of this paper, with selected parameters shown in the first row of Table 2 and joint posterior distributions in Appendix C. We also plot the contribution from the planet and star separately at their best-fit flux levels. We compute the auto-correlation function of the residuals and find that there is no evidence for correlated noise or strong systematics. Unless otherwise specified, we quote results from the baseline EddySed cloud model (MgSiO₃, am). See Appendix C for the posterior distributions of other parameters from our baseline model.

To make sure that we are fitting the correct signal, we check the RV and flux level of the companion. From our orbital posteriors for HD 4747 B, the expected RV shift on the night of our HRS observation is $15.0 \pm 0.1 \text{ km s}^{-1}$ in the Earth's reference frame, which is a combination of the system barycenter velocity, the Earth's relative velocity with respect to HD 4747, and the companion's orbital velocity. The fitted RV of $15.0 \pm 0.4 \text{ km s}^{-1}$ agrees perfectly with this value (see Figure 4). In addition, the companion flux level in the spectral orders from 2.29 to 2.49 µm is 85 ± 10 counts, comparable to the speckle flux levels in these orders. Taking into account the difference in wavelengths and the difference in integration time (600 s for the companion, 60 s for the on-axis star), we estimate that our measured companion flux corresponds to $\Delta K_s = 8.3 \pm 0.3$ mag, which is within $3\sigma$ of the photometric $\Delta K_s = 9.05 \pm 0.14$ mag reported by Crepp et al. (2018). The agreement between these contrast values is reasonably good given the time-varying throughput of KPIC (Delorme et al. 2021), and the fact that we subtract out the continuum with high-pass filtering, effects that complicate a direct flux comparison.

Figure 4 also shows the projected spin rate $v \sin i = 13.2_{-1.5}^{+1.4}$ km s⁻¹, which is comparable to the rotation rates observed for field brown dwarfs with similar spectral types (e.g., Konopacky et al. 2012). We also plot the retrieved quench pressure $P_{\text{quench}}$ in Figure 4, which indicates that the chemical reaction timescale becomes longer than the vertical mixing timescale
<figure><img src="image_7.png" /><figcaption>**Figure 3.** The KPIC HRS used in this study are plotted in black, with error bars inflated to the best-fit value in gray. A sample full model is shown in teal (dashed), and consists of the companion model in orange (which has been RV shifted and broadened), and the stellar model in blue to model the speckle contribution. The companion model shown does not include tellurics to focus on molecular features, but tellurics are included in our fits. The residuals are shown as gray points.</figcaption></figure>

at pressures lower than $P_{\text{quench}}$. Thus, disequilibrium chemistry is clearly affecting the atmosphere (see Section 5.3 for details).

We compute $T_{\text{eff}}$ by sampling from our posterior to generate low-resolution models over a large wavelength range (0.5–30 µm) and calculating the integrated flux. We then solve for $T_{\text{eff}}$ using the Stefan–Boltzmann law. When computing $T_{\text{eff}}$, we include opacities from Na and K, which are important sources of opacity near visible wavelengths. As shown in Table 2, the retrieved radius and $T_{\text{eff}}$ from HRS have broad distributions, which reflect the relatively weak luminosity constraints from the HRS ($\log(L_{\text{bol}}/L_{\odot}) = -4.33_{-0.25}^{+0.23}$). This is because the HRS is not flux-calibrated and we remove the continuum in our fits. Comparing to values of radius and $T_{\text{eff}}$ from previous work based on LRS (Crepp et al. 2018; Peretti et al. 2019), our retrieved values from the HRS retrievals are consistent at the $1\sigma-2\sigma$ level (see Table 2). We discuss the constraints on these parameters from the LRS in Section 6.2.

We compare our retrieved [C/H], [O/H], and C/O with those of the host star (see Section 2.1) in Figure 5. Our retrieved C abundance agrees well with the host star value, while the O abundance is lower by about $1\sigma$. This results in our retrieved C/O for the companion being higher by about $2\sigma$ compared to the stellar value. Here and elsewhere in the paper, we compute the “$\sigma$ difference” between two measurements by dividing the difference in the two median values by the quadrature sum of the uncertainties from both measurements. We discuss the implications of our measured abundances for HD 4747 B in Section 8.4.

## 5.2. Why Is Our KPIC HRS Insensitive to Clouds in HD 4747 B?

Clouds represent a significant source of uncertainty in many published models of substellar atmospheres (e.g., Burningham et al. 2017; Wang et al. 2020a). However, we find that the retrieved parameters from our KPIC HRS are insensitive to the choice of cloud model for HD 4747 B. As shown in Figure 4, the posteriors for radius, RV, $v\sin i$, C/O, [C/H], and quench pressure are nearly identical across the various models. The same is true for other parameters.

Table 2 shows that the different cloud models are indistinguishable for the KPIC HRS; the clear model fits as well as the cloudy models, with $B=0.61$, which does not pass the threshold of $B=10/0.1$ to be considered statistically favored/disfavored. This indicates that the data can be fitted adequately without clouds; indeed the cloud parameters for the EddySed models span their respective prior ranges almost uniformly as shown in Appendix C. As we will discuss in Section 6, the LRS shows that the atmosphere of HD 4747 B is cloudy. This implies that cloud opacity must be minimal at the pressures probed by our HRS.

To understand this, we plot in Figure 6 the retrieved P–T profiles (black and blue lines), cloud condensation curves (dashed lines), and emission contribution functions. The left and right panels show results from the HRS and LRS retrievals, respectively. The emission contribution function for HRS shows that we are sensitive to pressures ranging from a few bars, where the continuum forms, up to $\approx 10^{-2}$ bars in the cores
<figure><img src="image_8.png" /><figcaption>**Figure 4.** Posterior distributions for five key parameters from HRS retrievals of HD 4747 B, using the EddySed cloud model (MgSiO₃, “am” and “cd” in blue and red, respectively), and the clear model in purple. The titles on each histogram show the median and 68% credible interval for the baseline retrieval (MgSiO₃, am). Regardless of the cloud model used, the results agree well between different fits for the RV, v sin i, C/O, [C/H] (discussed in Section 5.1), and quench pressure (discussed in Section 5.3).</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 5.** Retrieved C and O abundances (relative to solar) and C/O for HD 4747 B in blue. The titles on each histogram showing the median and 68% credible interval. The red points show the stellar values from Section 2.1. The [C/H] agrees well, [O/H] is consistent within 1σ between the companion and star, and C/O is consistent at the 2σ level.</figcaption></figure>

of individual lines. Note that the contribution functions use the wavelength axes on the top. In the EddySed model, the cloud base is set at the intersection of the P–T profile and a given

cloud condensation curve (dashed lines). For MgSiO₃, this corresponds to a pressure of ≈10–20 bars when using our HRS-retrieved P–T profile. As the cloud mass fraction drops
<figure><img src="image_10.png" /><figcaption>**Figure 6.** The P–T profiles from our HRS retrieval (left panel) and LRS retrieval (right panel) with the baseline cloud model (MgSiO₃, am). In each panel, black lines show the best-fit profile and blue lines are 100 random draws from the posterior. We also show a cloudless Sonora P–T profile (Marley et al. 2021) with similar bulk properties as HD 4747 B in dotted gray. The condensation curves for MgSiO₃ and Fe clouds are plotted as dashed lines. We also overplot the emission contribution function as contours, which show the fraction of flux (darker indicates higher fraction) a given pressure layer contributes to the total flux at a given wavelength (Mollière et al. 2019). Thus, these use the wavelength axes, and not the temperature axes. The HRS is sensitive to the continuum forming around a few bars and line cores, which form up to 10⁻² bars. Over the same wavelength range of 2.29–2.49 µm, the LRS shows the continuum arises from ≈1–10 bars, consistent with the HRS.</figcaption></figure>

exponentially above the cloud base in the EddySed model (controlled by $K_{zz}$ and $f_{sed}$), we find that the cloud opacity decreases to negligible levels by the time we reach pressures of a few bars where the continuum forms. For this reason, we do not consider models with Fe clouds in our HRS retrievals, since the Fe cloud base forms even deeper than that of MgSiO₃.

Therefore, our KPIC HRS is insensitive to clouds because we cover both a relatively small wavelength range (2.29–2.49 µm) and a range where molecular opacities from H₂O, CO, and CH₄ are significant. The small wavelength range means that the cloud opacity is effectively constant in wavelength. The strong molecular opacity in HRS allows us to resolve many individual absorption lines and obtain good constraints on the atmospheric composition for molecules present in this region of the spectrum. The opacity of these molecules decreases at shorter wavelength due to decreasing excitation cross sections, so the continuum shifts to higher pressures (deeper down) at shorter wavelengths. This effect is visible in the LRS contribution function, where close to 1 µm, the emission originates from roughly the same pressure as the MgSiO₃ cloud base, making the *y* and *J* bands particularly sensitive to clouds (see Figure 10).

Could the KPIC HRS be affected by clouds at lower pressures (higher altitudes) than predicted by the EddySed model? Several studies have found that including clouds at lower pressures than predicted by EddySed produces better fits to mid-IR spectra of isolated brown dwarfs (e.g., Burningham et al. 2021; Luna & Morley 2021). As shown in Figure 6, our HRS P–T profiles show a nearly isothermal region between about 0.1–1 bars, which could suggest a degeneracy with clouds (Burningham et al. 2017). To check whether the P–T parameterization affects our results, we run a retrieval with a fixed P–T, namely the self-consistent profile overplotted in gray. We find that all posteriors from this fixed P–T fit overlap within 1σ with those from our baseline retrieval. Thus, we conclude that the isothermal part of the P–T we retrieve is not biasing our conclusions. To further examine the possibility of clouds at lower pressures, we also run an opaque cloud model with infinite opacity below a retrieved pressure, and a gray

cloud model that adds a constant cloud opacity at each pressure layer. When fitting the HRS with these more flexible cloud models, we also find consistent results with the baseline model. In the second model, the gray opacity is bounded to lie below ~0.03 cm² g⁻¹, and the pressure of the infinitely opaque cloud is required to be deeper than ~1 bar. Therefore, even with these less constraining cloud parameterizations, we find that our HRS still prefers solutions with minimal cloud opacity.

### 5.3. Disequilibrium Chemistry with Deep Quenching Pressure

In our retrievals, we include a simple model for disequilibrium chemistry using the quench pressure prescription in petitRADTRANS, which is motivated by Zahnle & Marley (2014). Specifically, the abundances of CH₄, CO, and H₂O are held constant at atmospheric pressures lower than the retrieved P-quench parameter. We find that when including quenching, the goodness of fit increases drastically compared to fits with full chemical equilibrium. For example, between two clear retrievals with and without quenching, we find that $B \approx 380$, in favor of the quenched retrieval. From the Bayes factor interpretation of Benneke & Seager (2013), this represents a detection of quenching at $\approx 3.9\sigma$ significance. The quench pressure retrieved is also highly consistent between retrievals with different cloud models, with 1σ and 2σ intervals of 50–260 and 14–836 bars, respectively (Figure 4). In this section, we explore reasons why the data prefer disequilibrium chemistry in the atmosphere of HD 4747 B. The physical implications of our retrieved quench pressure, including an estimate of the vertical diffusion coefficient ($K_{zz}$), are discussed in Section 8.2.

To understand why the data prefer a deep quench pressure, we plot lines of constant log(CO/CH₄) volume-mixing ratios (VMRs) along with the best-fit P–T profile from our baseline HRS retrieval in Figure 7. We calculate CO/CH₄ from this quenched chemistry retrieval by finding the abundances of each molecule in the chemical grid, iterating over our posterior distribution of C/O, [C/H], and P–T profile. We find that $CO/CH_4 = 13.6_{-4.6}^{+5.8}$. If the atmosphere was in chemical
<figure><img src="image_11.png" /><figcaption>**Figure 7.** The best-fit P–T profile from our baseline HRS retrieval in blue, with the thicker region indicating the retrieved quench pressure (2σ interval of 14–836 bars). Lines of constant log(CO/CH₄) VMRs are shown, with black lines (gray lines) indicating the region where CO (CH₄) is more abundant. The P–T profile nearly overlaps with the log(CO/CH₄) = +1 line below ~20 bars, which is where we retrieve the quench pressure to be from the HRS.</figcaption></figure>

equilibrium, we repeat our calculation and find that we would expect CO/CH₄ = 1.35<sup>+0.21</sup><sub>-0.17</sub>, which is 10 times smaller than our retrieved value in the quenched chemistry model. Thus, the relative under-abundance of CH₄ relative to CO in our HRS leads our models to prefer a deep quench pressure. The value of CO/CH₄ also determines our retrieved the quench pressure, whose 2σ interval is indicated by the thick blue region in Figure 7. Because the P–T profile nearly overlaps the curve of CO/CH₄ = 10 at ~20 bars and deeper, a broad range of quench pressures deeper than ~20 bars are consistent with the data.

As another way to visualize the detection of disequilibrium chemistry, we plot the molecular abundances in VMR as a function of pressure in Figure 8. The solid lines show the VMRs for the HRS quenched chemistry retrieval, while dashed lines show the VMRs for the same model with quenching turned off manually. By comparing the solid and dashed lines, we see that at the pressures probed by our observations, the relative abundances of CO, CH₄, and H₂O differ by several orders of magnitude between the quenched model and expectations from chemical equilibrium.

### 5.4. Detection of Methane in the HRS

In this section, we take a closer look at the relatively weak methane absorption signal in our HRS, which leads us to prefer quenched models where the CO/CH₄ ratio is a factor of 10 higher than predicted in models assuming chemical equilibrium. We confirm the presence of detectable levels of methane in the HRS by running a pair of free retrievals, one with only H₂O and CO, and one with H₂O, CO, and CH₄. The results of these retrievals are listed in Table 3. In these free retrievals, we fit the abundances of each absorbing species independently and assume a constant abundance as a function of pressure. Although we also considered models that included NH₃ and CO₂, we only obtained upper limits on their abundances, and therefore excluded them from our fits in this section. Finally, given the insensitivity of the HRS to clouds, we carry out these tests with the clear model to save computation time.

We find that the data strongly prefer the model with CH₄, with a Bayes factor of 84 (3.4σ significance; Benneke & Seager 2013).

<figure><img src="image_12.png" /><figcaption>**Figure 8.** Solid lines: best-fit VMRs of CO, H₂O, and CH₄ from our baseline HRS retrieval with chemical quenching. Dashed lines: the corresponding VMRs when quenching is turned off (i.e., equilibrium abundances) for the same retrieval. The 1σ quench pressure is indicated by the gray region, and the blue region shows schematically the pressures where our HRS is sensitive. Within the blue region, the relative CO/CH₄ ratio can differ by orders of magnitude between the quenched abundances and the equilibrium abundances.</figcaption></figure>

<table><caption>**Table 3**
Free Retrievals Carried Out on HD 4747 B for Validating the CH₄ Detection in the KPIC HRS (Section 5.4)</caption><thead><tr><th>Molecules</th><th>log(CO)</th><th>log(H<sub>2</sub>O)</th><th>log(CH<sub>4</sub>)</th><th><i>B</i></th></tr></thead><tbody><tr><td>CO, H<sub>2</sub>O</td><td>-3.51<sup>+0.21</sup><sub>-0.17</sub></td><td>-3.77<sup>+0.19</sup><sub>-0.16</sub></td><td>...</td><td>1</td></tr><tr><td>CO, H<sub>2</sub>O, CH<sub>4</sub></td><td>-3.42<sup>+0.22</sup><sub>-0.18</sub></td><td>-3.70<sup>+0.19</sup><sub>-0.16</sub></td><td>-4.82 ± 0.23</td><td>84</td></tr></tbody></table>

**Note.** We list the log volume-mixing ratios of molecules included, and the Bayes factor between the model with and without CH₄.

As shown in Table 3, we obtain log(CH₄) = -4.82 ± 0.23 from the free retrieval, and the CH₄ posterior in Figure 9 shows no strong covariance with the abundances of either CO or H₂O. If the atmosphere was in chemical equilibrium, we would expect a CH₄ VMR that is 10 times higher than what we retrieve, according to the same calculation described in Section 5.3. We note that the abundances from the free retrieval with CH₄ also agree well with the corresponding VMRs from our quenched chemistry retrievals. This is not surprising given the deep quench pressure we retrieve, which makes the molecular abundances constant in the regions where our HRS is sensitive (see Figure 8).

We separately visualize the CH₄ detection in cross-correlation space by carrying out an analysis similar to that described in Zhang et al. (2021). First, we make a “pure CH₄ template” from the best-fit companion model with CH₄, H₂O, and CO by manually setting the abundances of H₂O and CO to zero. If the model without CH₄ is fitting poorly due to its inability to fit CH₄ lines in the data, we would expect the residuals of this model, which we denote R = (data—model without CH₄), to contain CH₄ lines. Therefore, we cross-correlate R with the pure CH₄ template, plotted as the blue CCF in Figure 9. In addition, we plot the CCF of R' = (data—model with CH₄) with the pure CH₄ template in red for comparison. The blue CCF shows a peak at 0 km s⁻¹ (solid gray line), where we expect a real signal to be since the models have been shifted by the best-fit companion RV. If the residuals were dominated by telluric CH₄ for example, the CCF peak would
<figure><img src="image_13.png" /><figcaption>**Figure 9.** Left panel: joint posterior distributions of the log(VMR) of CO, H₂O, and CH₄ from a KPIC HRS retrieval where we directly fit the molecular abundances and assume they are constant across pressure. Right panel: CCFs of a pure CH₄ template with (KPIC data—model without CH₄) in blue, and the CCF of the CH₄ template with (data—model with CH₄) in red. The CH₄ template is generated with best-fit parameters of the full model with CH₄, CO, and H₂O and manually setting opacities of CO and H₂O to zero. The gray solid line indicates the companion rest frame, and the gray dotted line is the telluric rest frame. The residuals are taken from two spectral orders (2.29–2.41 µm) with stronger CH₄ detection.</figcaption></figure>

appear at the dotted gray line (negative of the RV, or −15 km s⁻¹). Thus, even though the height of the CH₄ peak in the blue CCF is small compared to the surrounding structure, the fact that it is located at the companion RV is evidence of a real signal from CH₄.

In our CCF framework, the y-axis is the estimated flux level (in counts) of the companion signal from a least-squares minimization. As shown in Figure 9, we find a flux level of ≈50 counts for CH₄, which is an estimate of the companion flux in the residuals. Importantly, this value is consistent with the flux value found when we repeat the same CCF analysis with H₂O (i.e., comparing a model with only CO and CH₄ and the baseline model of CO, H₂O, and CH₄). For a molecule such as NH₃, which we see no evidence of in the KPIC HRS, the flux value from the CCF becomes unbounded as the least-squares routine used for computing the CCF fails to converge.

Finally, we check for crosstalk between H₂O and CH₄ by cross-correlating *R* with the pure water template and detect no CCF peaks. Furthermore, we note that in a retrieval with only CO and CH₄ (no H₂O), the retrieved CH₄ abundance is consistent with the value from the full model including CO, H₂O, and CH₄.

We therefore conclude that the data strongly favor the presence of detectable levels of methane in the HRS, with an abundance significantly lower than that predicted by equilibrium chemistry models. The detection of methane at log(CH₄) = −4.82 ± 0.23 demonstrates the ability of KPIC to retrieve species that are more than an order of magnitude lower in VMR than the dominant molecular constituents in the data, in only 1 hr of integration time.

# 6. Low-resolution Retrievals (GPI + SPHERE)

## 6.1. Overview

In this section, we present the results from our fits to the LRS and compare our retrieved parameters to those from the HRS

fits. We fit the LRS using the same models as before. These include one clear model and four different implementations of the EddySed cloud model where we vary our assumptions about the unknown cloud properties. The cloudy models consist of two MgSiO₃ retrievals with am and cd particles (explained in Section 4.2.2), and two retrievals with MgSiO₃ and Fe clouds (again, am and cd). In Figure 10, we plot the data, best-fit cloudy and clear models, the residuals, and the GP models of the residuals. The posteriors for a few key parameters from these retrievals are plotted in Figure 11 and tabulated in Table 2. See Appendix C for the posterior distributions of other parameters in the baseline model.

When comparing the clear and cloudy models in Figure 10, we see that the data shortward of ~1.2 µm is poorly fit by models without clouds. This causes the clear model to have $B \approx 7 \times 10^{-26}$; it is overwhelmingly ruled out compared to the baseline EddySed model. In addition, when we plot the models over a larger wavelength range in Figure 12, we find that the cloudy models agree with the NIRC2 L-band photometry from Crepp et al. (2016), while the clear model overpredicts the L-band flux by ≈2σ. We did not include these photometric points in our retrievals.

Figure 10 shows that the SPHERE J-band data from ≈1.2–1.35 µm is not well fit by even the cloudy model, which could either be caused by model mismatch or speckle contamination that artificially raises the flux. The GP model finds that ~60% of the SPHERE error bars and ~90% of the GPI error bars are from correlated noise, with correlation length scales of ~6 and ~2 wavelength channels, respectively. This confirms our initial intuition that the noise in the SPHERE and GPI images is likely dominated by correlated speckle noise based on visual inspection of the images. For the SPHERE data set, we estimate that the retrieved length scale is roughly equal to the number of steps that a speckle would move across the PSF for our brown dwarf’s separation; indeed, we see speckles moving across the companion PSF in the reduced images.
<figure><img src="image_14.png" /><figcaption>**Figure 10.** Top panel: best-fit cloudy model (EddySed + MgSiO₃, am) in blue and random models drawn from the posterior in light orange for an LRS retrieval of HD 4747 B. The best-fit clear model is overplotted in dashed purple, which fits visibly worse from ≈1.0 to 1.2 µm. Bottom panel: the residuals of the cloudy model shown in error bars, along with random draws of the GP models for the residuals in orange.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Figure 11.** Posterior distributions for a few key parameters from LRS retrievals of HD 4747 B, using the EddySed model with MgSiO₃ clouds (blue: amorphous; red: crystalline), and MgSiO₃ + Fe clouds (purple: amorphous, yellow: crystalline). The titles on each histogram show the median and 68% credible interval for the MgSiO₃, am model. The results disagree by as much as 3σ, especially in 2D space, and display strong covariance between C/O and [C/H]. The radius retrieved is also generally smaller than predicted by evolutionary models.</figcaption></figure>

Overall, the SPHERE spectrum is less reliable than that from GPI because only four exposures are available, compared to the ~40 exposures from GPI.

Finally, the P–T profile retrieved from our baseline LRS retrieval shows a bimodal distribution (see Figure 6). The

<figure><img src="image_16.png" /><figcaption>**Figure 12.** Random models drawn from the posterior of the baseline LRS retrieval (MgSiO₃, am), plotted over a larger wavelength range and color-coded by [C/H], the metallicity. There is a gradient in [C/H] in the L (≈3.4–4.2 µm) and M (≈4.55–4.8 µm) bands, which can be distinguished with a comparable S/N LRS in these bands. The GPI and SPHERE data are shown in black, and we also overplot the photometric data points from Crepp et al. (2016), which are not included in the retrievals but nonetheless agree with our models.</figcaption></figure>

degeneracy seen here may be related to issues with the LRS (see Section 6.3).

## 6.2. Comparison with Prior Knowledge

Because the LRS is flux-calibrated, we can check whether our retrieved radii and effective temperatures are physical and consistent with prior knowledge for this benchmark companion. Using the known age and mass of HD 4747 B ($3 \pm 2$ Gyr and $m = 67.2 \pm 1.8 M_{Jup}$), we interpolate the COND evolutionary models (Baraffe et al. 2003) to find a model-predicted radius of $0.8_{-0.03}^{+0.07} R_{Jup}$, and a predicted $T_{eff} = 1450_{-180}^{+350}$ K. Peretti et al. (2019) compared the SPHERE spectrum of HD 4747B to those of field brown dwarfs to derive a more tightly constrained $T_{eff} = 1350 \pm 50$ K (see Table 2), which we adopt in the subsequent discussion.
We calculate the effective temperatures of the models in our retrievals by integrating the flux over 0.5–30 µm. For the baseline EddySed model, we find $T_{\text{eff}} = 1473_{-20}^{+17}$ K, and a radius of $0.70_{-0.03}^{+0.05} R_{\text{Jup}}$. Compared to prior expectations, however, the radius retrieved is too small by ≈1.5σ while $T_{\text{eff}}$ is too high by ≈2σ. From substellar evolutionary models, the minimum possible radius of a brown dwarf should be ≈0.74 $R_{\text{Jup}}$, which is imposed by electron degeneracy pressure (Chabrier et al. 2009). We find that $T_{\text{eff}}$ and radius are correlated in the LRS retrievals, as shown in Figure 11, which is expected as different combinations of these two parameters can produce the same total luminosity. However, our total luminosity agrees well with the luminosity predicted by evolutionary models.

Several previous retrieval studies have also found smaller-than-expected radii for L dwarfs, which may be attributed to the presence of heterogeneous surface features, such as patchy clouds, that are not captured in current 1D retrieval frameworks (e.g., Gonzales et al. 2020; Kitzmann et al. 2020; Burningham et al. 2021). On the other hand, Gonzales et al. (2021) retrieved a radius consistent with evolutionary models for a seemingly cloudless L dwarf. Whether the radii from evolutionary models are correct is an assumption that is now being tested by a growing sample of transiting brown dwarfs from TESS (e.g., Carmichael et al. 2020).

In our retrievals with both MgSiO$_3$ and Fe clouds, we retrieve slightly larger radii that are more consistent with the evolutionary model prediction. This could indicate that a single cloud model (MgSiO$_3$) may be inadequate in attenuating the flux from the deep atmosphere. However, models with two cloud species do not improve the fit significantly ($B = 1.5–3$ compared to the baseline model with MgSiO$_3$ only). Furthermore, the MgSiO$_3$, cd model actually has the largest retrieved radius, but our data cannot distinguish between crystalline and amorphous particles. We conclude that our retrieved radius is sensitive to aspects of the cloud models that are poorly constrained by the existing data for this object.

## 6.3. LRS at Longer Wavelengths Could Improve Abundance and Cloud Constraints

While the LRS can provide tighter constraints on the cloud parameters and radius compared to the HRS, we find that many retrieved parameters, including the atmospheric abundances, are very sensitive to model choices. In Figure 11, we overplot the posterior distributions of a few parameters from our four EddySed models. The retrieved C/O and [C/H] have large uncertainties and can disagree at the 3σ level between models. The values also span a significant portion of the parameter space ($>1$ dex in metallicity), and show much stronger covariance compared to those measured from the HRS (see Figure 4). However, all cloudy models fit the LRS well, with Bayes factors within a factor of ~3 (see Table 2), so we cannot distinguish between them.

We note that Mollière et al. (2020) were able to obtain much better constraints on the composition of HR 8799e, which also has a cloudy atmosphere, using LRS data sets from 0.95–2.5 µm. Their LRS had an S/N between 4 and 11 per wavelength bin, much lower than the S/N of our LRS (between 20 and 60 per wavelength bin). Unlike Mollière et al. (2020), however, our study does not have LRS in the second half of K band (2.2–2.5 µm), which contains a strong CO bandhead as well as significant H$_2$O and CH$_4$ opacities. When we compute

the CO abundances from our LRS retrievals, we find that they are not well constrained, with $1\sigma$ intervals that are ≳3 wider than the CO constraint from HRS. In Figure 12, we plot random draws of our baseline model color-coded by metallicity out to 5 µm. As shown, the models diverge quickly in the 2.2–2.5 µm range. The fact that we miss this crucial wavelength region could explain why Mollière et al. (2020) obtained more robust constraints on atmospheric abundances and P–T profiles that agree better with self-consistent models than we do, despite using data with a lower S/N.

Figure 12 also shows a clear gradient in metallicity beyond 2.5 µm. In some of our cloudy LRS retrievals, we see a covariance between metallicity and cloud mass fraction, where lower metallicities correspond to higher cloud mass fractions as well as larger, more physically consistent radii (see Figure 11). The degeneracy between metallicity and cloud mass fraction might arise because both molecular opacities and clouds contribute opacity, and our data has insufficient wavelength coverage to probe more regions where the gas and cloud opacities are sufficiently different. From the LRS retrievals, we consistently find a factor of ~2–3 more CH$_4$ and H$_2$O than observed in the HRS, implying that the LRS retrievals could be compensating for our imperfect cloud models by increasing the gas opacities.

Using a more flexible cloud model might alleviate some of these issues. For example, Burningham et al. (2021) retrieved the 1–15 µm LRS of a field L dwarf and found the data preferred silicate clouds much higher up than the predicted cloud base locations from condensation curves. In addition, their retrieved cloud particles also have smaller sizes (submicron) than predicted by the EddySed model (a few microns). Similarly, Luna & Morley (2021) found that submicron cloud particles at lower pressures than predicted by EddySed are required to fit the mid-IR silicate feature (≈8–10 µm) of many L dwarfs. They found that the microphysical cloud model CARMA (Turco et al. 1979; Toon et al. 1988; Gao et al. 2018) allows them to fit their data much better and even place constraints on which cloud species are producing the observed absorption feature.

Both of the above-mentioned studies benefitted from data at ~10 µm that significantly helped with constraining cloud properties. Thus, to obtain better abundance measurements with LRS, it is not only important to obtain full coverage in the near-IR (which we lack), but also to acquire data in the mid-IR. JWST can obtain low- and medium-resolution spectroscopy of brown dwarfs spanning the near- to mid-IR wavelengths using the NIRSpec and MIRI instruments. Future ground-based instruments such as SCALES at Keck (Stelter et al. 2020) will also provide LRS in the mid-IR.

## 7. Joint Retrievals

In this section, we describe the results of joint retrievals to both the HRS and LRS for HD 4747 B. In practice, we set up two radiative transfer routines with petitRADTRANS using line-by-line (for HRS) and correlated-k (for LRS) opacity sampling, respectively. The HRS and LRS models share the same atmospheric parameters and priors, but each has some unique parameters (e.g., RV and v sin i for HRS, GP kernel parameters for LRS). Within one nested sampling retrieval, we add the log likelihoods from the HRS and LRS components at each step of sampling to get the total log likelihood. We
<figure><img src="image_17.png" /><figcaption>**Figure 13.** Filled areas: $1\sigma$ intervals for the CO, $H_2O$, and CH$_4$ abundances from our KPIC HRS retrieval. Hatched areas: the same for a joint retrieval (HRS + LRS). The retrieved abundances are two to three times higher in the joint retrieval, while the relative abundance ratios between species stays roughly the same (which produces a similar C/O). This highlights the fact that HRS is better at constraining relative abundances than absolute abundances.</figcaption></figure>

consider both clear and cloudy EddySed models for our joint retrievals.

Because the LRS prefer clouds, the cloudy model ($MgSiO_3$, am) is overwhelmingly preferred in our joint retrieval, with a Bayes factor in excess of $10^{34}$ compared to the clear model. From the cloudy model, we retrieve $C/O = 0.70 \pm 0.03$ and $[C/H] = 0.34 \pm 0.07$. The retrieved uncertainties on these parameters are lower than in the HRS-only retrieval (which had $C/O = 0.66 \pm 0.04$ and $[C/H] = -0.10^{+0.18}_{-0.15}$). In addition, the C/O from our joint fit is consistent with the C/O from our HRS fit. This is not surprising, because the HRS places tight constraints on the relative line depths (and hence the relative abundance ratios) of CH$_4$, $H_2O$, and CO. However, the joint fit pushes the metallicity to higher values, which corresponds to increased gas abundances as shown in Figure 13. The joint fit results translate to a $>4\sigma$ discrepancy in $[C/H]$ between HD 4747 A&B, while there is no discrepancy if we take the results from the HRS fit. This implies that the joint fit might be compensating for inadequacies in modeling clouds by increasing the gas opacities, as discussed in Section 6.3 for the LRS-only case. We ran additional joint retrievals where we varied the cloud parameters (e.g., adding Fe clouds) and found similar results.

If we compare the log likelihoods of the HRS part of the joint fit to that from the HRS-only fit, we find that the HRS is fit less well by $\approx e^{10}$ (which translates to $\approx 4\sigma$) in the joint fit, implying a trade-off between fitting the LRS and HRS. We can qualitatively compare the LRS S/N per wavelength bin to the CCF S/N of the HRS, which approximates the total constraining power of the HRS. When including all molecules in our model, we find a CCF S/N of $\approx 15$ for the HRS. For the LRS, the average S/N per wavelength bin is $\approx 20$ for the SPHERE YJH data and $\approx 60/30$ between the GPI H/K$_1$ data. This explains why the joint fits prioritize fitting the LRS at the expense of the HRS.

As discussed in Section 6, the LRS are very model-sensitive and additionally contaminated by correlated noise. For this reason, we adopt the HRS-only results as the best estimate of HD 4747 B's atmospheric properties in this paper (see first row of Table 2). We leave it to future work, preferably aided by longer wavelength coverage in LRS, to achieve a more satisfactory joint retrieval.

# 8. Discussion

## 8.1. Next Steps for High-resolution Spectroscopy

Our KPIC HRS provides a better handle on the atmospheric abundances of HD 4747 B, and is less sensitive to model choices than our LRS. In fact, our K-band HRS is essentially agnostic to clouds in the brown dwarf's atmosphere; all retrieved parameters are consistent independent of our chosen cloud model (Figure 4). As discussed in Section 5.2, this is because our HRS covers a wavelength region (2.29-2.49 µm) of high molecular opacity, and probes emission across a wide range of atmospheric pressures where cloud opacity is negligible (up to $10^{-2}$ bars in line cores). While clouds affect the continuum near 1 µm in the LRS, they have little effect on the line depths across the wavelength range of our HRS. The relative line depths are sensitive to relative molecular abundances, which directly constrain C/O. These results advocate for using HRS to measure atmospheric abundances.

In the future, it is important to explore whether these findings hold true for other substellar objects. In upcoming papers, we will present KPIC HRS retrievals of brown dwarf companions and giant planets spanning a range of effective temperatures and surface gravities. Ultimately, it would also be useful to constrain cloud properties with HRS. For transmission spectroscopy, Gandhi et al. (2020) found that their simulated near-IR HRS for warm Neptunes is more sensitive to molecular abundances than LRS for the same reasons highlighted in this study. While both clouds and metallicity affect the line depths in HRS, Gandhi et al. (2020) showed that increasing the wavelength coverage (e.g., going from 0.9-1.7 µm to 0.9-2.5 µm) helps distinguish between clouds and metallicity and provides better constraints on both. Thus, if we wish to obtain constraints on clouds and abundances at the same time, it would be important to extend our current HRS to a broader range of wavelengths. KPIC Phase II will allow us to obtain L-band data ($\approx 3.4-4.1$ µm) to complement existing K-band data (Delorme et al. 2021), and future upgrades could benefit from including the H and J bands as well.

In this study, we have assumed that the atmosphere of HD 4747 B is globally uniform. However, it would be important to examine the impact of 3D effects, including nonuniform cloud coverage. Past studies with photometry or LRS show that many brown dwarfs exhibit clear rotational variability signals (e.g., Apai et al. 2013; Biller et al. 2018; Zhou et al. 2018; Manjavacas et al. 2021; Vos et al. 2022), which appear to be caused by inhomogeneities in their cloud properties. Therefore, time-resolved observations are important to understanding clouds and 3D effects.

With HRS, we can use the time-varying line depth and shape to map the 2D brightness distributions of these objects (e.g., Crossfield et al. 2014). In this paper, we used 1 hr of KPIC data for HD 4747 B. Given our measured $v \sin i$ and assuming a radius of $0.8R_{\text{Jup}}$, we would expect a 5 or 7 hr rotation period if $i$ is equal to the orbital inclination or $i = 90^\circ$. Thus, it may be possible to sample a full rotation period within a single observing night, with the caveat that measurements of the true rotation period remain difficult for high-contrast companions (Biller et al. 2021).

## 8.2. Methane and the Presence of Disequilibrium Chemistry

Our HRS retrievals indicate that the ratio of CO/CH$_4$ (VMR) is $\approx 10$ times higher than expected by equilibrium
chemistry (see Section 5.3). To gain more physical intuition, we convert the quench pressure from our HRS retrievals to an estimate of the vertical diffusion coefficient, $K_{zz}$. To do this, we match the chemical timescale of the CO-CH$_4$ reaction from Zahnle & Marley (2014) with the mixing timescale $\tau_{\text{mix}} = L^2/K_{zz}$. While the length scale $L$ is typically taken to be the pressure scale height $H$, for lack of a better estimate, Smith (1998) showed that this assumption is not valid across several reactions in the atmospheres of Jupiter and Neptune. In fact, Smith (1998) found that $L \approx 0.1H$, which changes the inferred $K_{zz}$ by two orders of magnitude. Similarly, Ackerman & Marley (2001) also noted that the mixing length is generally shorter than the pressure scale height $H$ in stable atmospheric regions. Due to the uncertainty in $L$, we adopt $L = \alpha H$, where $\alpha$ is a scaling factor, and $H = \frac{k_B T}{\mu m g}$ ($\mu$: mean molecular weight, $g$: surface gravity, $T$: the local temperature). For each value of quench pressure from our posteriors, we compute the necessary quantities to derive a posterior for $K_{zz}$. For instance, if $\alpha = 0.1$, we find $K_{zz} = 5 \times 10^8 - 1 \times 10^{12} \text{ cm}^2 \text{ s}^{-1}$ ($1\sigma$). On the other hand, if $\alpha = 1$, we obtain $K_{zz} = 5 \times 10^{10} - 1 \times 10^{14} \text{ cm}^2 \text{ s}^{-1}$.

There have been few quantitative measurements of $K_{zz}$ for substellar companions. Miles et al. (2020) used M-band LRS to constrain the CO abundance and estimate $K_{zz}$ for seven field brown dwarfs. However, their objects have $T_{\text{eff}}$ between 250 and 750 K, much colder than HD 4747 B. In terms of objects with $T_{\text{eff}} \gtrsim 1000$ K, Barman et al. (2015) reported a detection of CH$_4$ in HR 8799b ($T_{\text{eff}} \sim 1000$ K) with Keck/OSIRIS data, which they used to estimate $K_{zz}$ between $10^6$ and $10^8 \text{ cm}^2 \text{ s}^{-1}$. However, this CH$_4$ detection was not confirmed by an independent study (Petit dit de la Roche et al. 2018), and recently Ruffio et al. (2021) concluded that future higher-resolution follow-up is needed to resolve the discrepant CH$_4$ signal strengths found by different analyses. Ruffio et al. (2021) pointed out that if the CH$_4$ abundance was overestimated by Barman et al. (2015), that would imply a larger $K_{zz}$. Using LRS, Mollière et al. (2020) found a well-constrained quench pressure for HR 8799e ($T_{\text{eff}} \sim 1100$ K) from *petitRADTRANS* retrievals, which could similarly be converted to a $K_{zz}$ constraint. In summary, our finding HD 4747 B, which is $\sim 300–400$ K hotter than HR 8799 b and c and much older (a few Gyr; Section 2.1) than most directly imaged planets, represents an important new data point because hotter objects are expected to be closer to equilibrium, making chemical disequilibrium processes harder to detect (e.g., Moses et al. 2013).

Zahnle & Marley (2014) provided an upper limit for $K_{zz}$ from mixing length theory (Gierasch & Conrath 1985) assuming full convection. For HD 4747 B, their Equation (4) translates to an upper limit of $\approx 10^9 \text{ cm}^2 \text{ s}^{-1}$. Depending on $\alpha$, our retrieved $K_{zz}$ either exceeds this upper limit by $\gtrsim 2\sigma$ (if $L=H$), or is very close to this limit (if $L=0.1H$). Together, this suggests that convection is driving the vertical mixing in HD 4747 B, and that the mixing efficiency is likely close to its predicted maximum. We check whether our inferred $K_{zz}$ makes sense by comparing them to those predicted by self-consistent atmospheric models with disequilibrium chemistry from Karalidi et al. (2021) and Mukherjee et al. (2022). For an object with properties similar to HD 4747 B, our measured CH$_4$ VMR is consistent with $K_{zz} \sim 10^8–10^{12}$ in these models (with the assumption that $L=H$). These values of $K_{zz}$ are roughly consistent with our estimate based on $P_{\text{quench}}$, and also near the upper limit from Zahnle & Marley (2014). On the modeling

front, it would be valuable to carry out 3D hydrodynamical simulations (e.g., Zhang & Showman 2018; Tan & Showman 2021) of brown dwarf interiors to independently estimate $K_{zz}$ (Tan 2022) and compare the results to that inferred by our data. Such simulations could also reveal which physical processes might cause a discrepancy between mixing length theory and our observations.

### 8.3. Dynamical versus Spectroscopic Mass Constraints

For a majority of substellar companions observed by direct imaging, there are no dynamical mass constraints. To assess whether our mass prior plays an important role in the results, we repeat our HRS and LRS retrievals with the baseline cloud model but use uniform priors in mass from 10 to 100 $M_{\text{Jup}}$ ("free-mass"). For the HRS free-mass retrieval, we find that all parameters change by less than $1\sigma$ compared to the mass-prior retrieval. The mass itself shows a broad distribution (33–76 $M_{\text{Jup}}$ at $1\sigma$) that encompasses the dynamical mass. Because our KPIC HRS is not flux-calibrated, the radius is not well constrained. In this case, we get large uncertainties in the spectroscopic mass because mass is inferred from the retrieved surface gravity, which depends on the poorly constrained radius.

Our LRS free-mass retrieval also yields posteriors for all parameters that are consistent between $1\sigma$ and $2\sigma$ with the mass-prior retrieval. Furthermore, the mass retrieved by the LRS is $59_{-8}^{+7} M_{\text{Jup}}$, which agrees within about $1\sigma$ with the dynamical mass. This provides confidence that reasonable mass constraints can be placed on substellar objects from LRS. The radius retrieved is $0.77 \pm 0.03 R_{\text{Jup}}$, consistent with evolutionary model predictions and close to the radius from the mass-prior retrieval, suggesting the two retrievals find a similar surface gravity.

### 8.4. Atmospheric Abundances of HD 4747 AB

We retrieve [C/H] and [O/H] values that are $1\sigma$ consistent with those of the host star, as discussed in Section 5.1. Both the companion and the star are mildly sub-solar in terms of their metal content. However, our retrieved C/O = $0.66 \pm 0.04$ is higher by approximately $2\sigma$ than the stellar C/O = $0.48 \pm 0.08$.

The question is whether the marginal discrepancy in C/O is from astrophysical or systematic reasons. For example, Wang et al. (2022) carried out retrieval experiments on simulated HRS (2.2–2.35 µm, $R=35,000$) and found that their formal error bars are likely underestimated due to systematic errors at the $\sim 0.15$ level in C/O. Using KPIC HRS from 2.23–2.33 µm, they found $\approx 1\sigma–1.5\sigma$ discrepancies between the [C/H] and [O/H] abundances of HR 7672 A and B, another benchmark brown dwarf system. On an earlier study of benchmark brown dwarfs, Line et al. (2015) quoted $1\sigma$ uncertainties of 0.2–0.3 in their brown dwarf C/O (much larger than our formal C/O uncertainty of 0.04), and concluded that a $2\sigma$ agreement between the stellar and companion C/O is sufficiently good given the caveats. It is also possible that the uncertainties on stellar abundances are underestimated given non-LTE effects (Line et al. 2015).

Another factor that might contribute to the $2\sigma$ discrepancy in C/O is uncertainties in the chemistry of condensates. The chemical model of *petitRADTRANS* we use accounts for the equilibrium condensation of various species and reports the
global (rather than gas phase) C and O abundances (Mollière et al. 2019). In particular, species such as MgSiO₃ and Mg₂SiO₄ contain three or four oxygen atoms per molecule, and are expected to hold a significant portion of O (Line et al. 2015). From our HRS retrievals, we find that ≈18% of O is condensed into solids such as MgSiO₃. In order to decrease the global C/O of the brown dwarf by ≈0.1 (therefore making the companion and stellar C/O agree at the 1σ level), we require a ~20% increase in the net O abundance. Keeping everything else unchanged, this means the MgSiO₃ mass fraction, which is predicted by the chemical model to be ~2 × 10⁻³ in our retrievals, needs to be doubled to ~4 × 10⁻³. From the LRS retrievals, the cloud base MgSiO₃ fraction can be as high as 10⁻². Therefore, a factor of ~2 uncertainty in the abundance of MgSiO₃ could make our C/O consistent at the 1σ level with the stellar value.

Given these caveats, we conclude that the 2σ difference between our retrieved C/O for HD 4747 B and the stellar value is not significant, and HD 4747 AB are consistent with being chemically homogeneous. Chemical homogeneity is expected by models where brown dwarf companions form via gravitational fragmentation in molecular clouds (e.g., Padoan & Nordlund 2004) or massive protostellar disks (e.g., Stamatellos et al. 2007). Simulations suggest that brown dwarfs typically form as part of unstable, high-order multiple systems, which undergo chaotic interactions that reduce the multiplicity over time (e.g., Bate et al. 2002; Thies et al. 2010; Bate 2012). With a semimajor axis of 10 au, HD 4747 B is unlikely to have been directly affected by such encounters, but its relatively high orbital eccentricity (≈0.73) could encode such a dynamically “hot” past.

## 9. Conclusions

Using high-resolution spectrum (R ~ 35,000) obtained by Keck/KPIC, we retrieve [C/H] = −0.10<sup>+0.18</sup><sub>-0.15</sub>, [O/H] = −0.18<sup>+0.18</sup><sub>-0.15</sub>, and C/O = 0.66 ± 0.04 for the benchmark brown dwarf companion HD 4747 B (formal error bars). The C and O abundances are consistent with the stellar values to ≲ 1σ, while the C/O ratio is consistent at the 2σ level, as expected for a formation scenario akin to that of binary stars. This shows that we can measure the atmospheric abundances for high-contrast substellar companions to the 20% level with KPIC and our current modeling framework, which Wang et al. (2022) also showed for another benchmark brown dwarf. We outline some other key findings from our study below.

We measure precise abundances from the KPIC HRS (2.29–2.49 µm), which are insensitive to our choice of cloud model. Our abundance measurements suggest that HD 4747 B has a CO/CH₄ ratio that is 10 times higher than predicted by equilibrium chemistry, corresponding to a quench pressure of 50–260 bars (1σ). This translates to a high vertical diffusion coefficient K<sub>zz</sub>, which depends on the assumed mixing length scale L. However, even if L is 10 times smaller than the pressure scale height, we get K<sub>zz</sub> = 5 × 10<sup>8</sup> – 1 × 10<sup>12</sup> cm<sup>2</sup> s<sup>−1</sup>, which implies a mixing strength that is at or above the upper limit predicted by mixing length theory.

The composition retrieved from our LRS (1–2.2 µm) is both sensitive to model choices and can be biased by the presence of speckles. For this reason, HRS provides a more reliable picture of the atmospheric composition in the current data sets, although the LRS could be improved with additional observation at longer wavelengths including the L and M bands. Despite these challenges, the current LRS does provide a spectroscopic mass estimate that is 1σ consistent with the dynamical mass for the brown dwarf.

Although our joint retrieval results are likely biased by the limited LRS wavelength coverage, joint analyses of LRS and HRS remain a promising avenue to constrain cloud properties and abundances simultaneously and provide a more complete picture of substellar atmospheres. When extended wavelength coverage is available, it would also be important to consider possible 3D effects, including patchy clouds. These might be constrained by obtaining multiple spectra sampling a rotation period. Additional modeling work on condensation, chemistry, and vertical mixing rates is also important to inform future observational results.

We thank the referee for helpful comments that improved the manuscript. J.X. thanks Michael Zhang for advice in computing the CH₄ opacities, Jonathan Fortney for discussions on chemical disequilibrium, and Konstantin Batygin for discussions on brown dwarf formation.

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. This research has benefitted from the SpeX Prism Library maintained by Adam Burgasser at http://www.browndwarfs.org/spexprism. Funding for KPIC has been provided by the California Institute of Technology, the Jet Propulsion Laboratory, the Heising-Simons Foundation, the Simons Foundation, and the United States National Science Foundation grant No. AST-1611623. A.V. acknowledges funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program, grant agreement Nos. 757561 (HiRISE). H.A.K. acknowledges support from the President’s and Director’s Research & Development Fund Program, which is jointly funded by the Jet Propulsion Laboratory and the California Institute of Technology under a contract with the National Aeronautics and Space Administration. The computations presented here were conducted in the Resnick High Performance Center, a facility supported by Resnick Sustainability Institute at the California Institute of Technology.

Facilities: Keck (KPIC).

Software: petitRADTRANS (Mollière et al. 2019), dynasty (Speagle 2020), PyAstronomy (https://github.com/sczesla/PyAstronomy).

## Appendix A
### Orbit Fits for HD 4747 B

Our orbit fit for the HD 4747 system is shown in Figure 14.
<figure><img src="image_18.png" /><figcaption>**Figure 14.** Results from a joint fit to host star radial velocity (top left), relative astrometry (top right), and absolute astrometry (bottom panel) for the HD 4747 system. The data together constrain the orbital parameters and mass of both the companion and host star well (Table 1).</figcaption></figure>

**Appendix B**

**Extracted Low-resolution Spectrum and GPI Astrometry**

Our extracted spectrum for HD 4747 B based on observations with GPI (Crepp et al. 2018) and SPHERE (Peretti et al. 2019) are given in Table 4. Our relative astrometry measurements based on the GPI data are listed in Table 5.
<table><caption>**Table 4**

Extracted Low-resolution Spectrum for HD 4747 B</caption><thead><tr><th>Wavelength (μm)</th><th>Flux (10<sup>-15</sup> × W m<sup>-2</sup>/μm)</th><th>Flux Error (10<sup>-15</sup> × W m<sup>-2</sup>/μm)</th></tr></thead><tbody><tr><td colspan="3">SPHERE (YJH)</td></tr><tr><td>1.008</td><td>0.726</td><td>0.167</td></tr><tr><td>1.026</td><td>0.789</td><td>0.131</td></tr><tr><td>1.044</td><td>0.767</td><td>0.088</td></tr><tr><td>1.063</td><td>0.806</td><td>0.085</td></tr><tr><td>1.081</td><td>0.839</td><td>0.086</td></tr><tr><td>1.098</td><td>0.88</td><td>0.073</td></tr><tr><td>1.116</td><td>0.848</td><td>0.071</td></tr><tr><td>1.138</td><td>0.823</td><td>0.059</td></tr><tr><td>1.159</td><td>0.907</td><td>0.056</td></tr><tr><td>1.176</td><td>1.006</td><td>0.061</td></tr><tr><td>1.194</td><td>1.114</td><td>0.054</td></tr><tr><td>1.212</td><td>1.213</td><td>0.05</td></tr><tr><td>1.23</td><td>1.33</td><td>0.049</td></tr><tr><td>1.248</td><td>1.338</td><td>0.06</td></tr><tr><td>1.267</td><td>1.378</td><td>0.061</td></tr><tr><td>1.285</td><td>1.389</td><td>0.06</td></tr><tr><td>1.302</td><td>1.33</td><td>0.063</td></tr><tr><td>1.318</td><td>1.188</td><td>0.051</td></tr><tr><td>1.33</td><td>1.106</td><td>0.061</td></tr><tr><td>1.34</td><td>0.875</td><td>0.077</td></tr><tr><td>1.358</td><td>0.506</td><td>0.086</td></tr><tr><td>1.399</td><td>0.482</td><td>0.106</td></tr><tr><td>1.418</td><td>0.596</td><td>0.078</td></tr><tr><td>1.432</td><td>0.598</td><td>0.053</td></tr><tr><td>1.447</td><td>0.668</td><td>0.042</td></tr><tr><td>1.463</td><td>0.725</td><td>0.033</td></tr><tr><td>1.479</td><td>0.741</td><td>0.034</td></tr><tr><td>1.495</td><td>0.788</td><td>0.038</td></tr><tr><td>1.511</td><td>0.912</td><td>0.036</td></tr><tr><td>1.526</td><td>0.987</td><td>0.036</td></tr><tr><td>1.54</td><td>1.059</td><td>0.038</td></tr><tr><td>1.553</td><td>1.117</td><td>0.038</td></tr><tr><td>1.568</td><td>1.179</td><td>0.04</td></tr><tr><td>1.582</td><td>1.217</td><td>0.04</td></tr><tr><td>1.599</td><td>1.266</td><td>0.044</td></tr><tr><td>1.613</td><td>1.235</td><td>0.041</td></tr><tr><td colspan="3">GPI (H)</td></tr><tr><td>1.506</td><td>0.789</td><td>0.022</td></tr><tr><td>1.51</td><td>0.798</td><td>0.022</td></tr><tr><td>1.516</td><td>0.841</td><td>0.024</td></tr><tr><td>1.522</td><td>0.89</td><td>0.025</td></tr><tr><td>1.531</td><td>0.936</td><td>0.026</td></tr><tr><td>1.539</td><td>0.983</td><td>0.028</td></tr><tr><td>1.547</td><td>1.046</td><td>0.031</td></tr><tr><td>1.554</td><td>1.107</td><td>0.031</td></tr><tr><td>1.562</td><td>1.153</td><td>0.032</td></tr><tr><td>1.572</td><td>1.183</td><td>0.033</td></tr><tr><td>1.581</td><td>1.218</td><td>0.033</td></tr><tr><td>1.589</td><td>1.233</td><td>0.034</td></tr><tr><td>1.597</td><td>1.22</td><td>0.033</td></tr><tr><td>1.605</td><td>1.216</td><td>0.032</td></tr><tr><td>1.613</td><td>1.21</td><td>0.031</td></tr><tr><td>1.621</td><td>1.201</td><td>0.032</td></tr><tr><td>1.63</td><td>1.21</td><td>0.033</td></tr><tr><td>1.638</td><td>1.208</td><td>0.032</td></tr><tr><td>1.646</td><td>1.183</td><td>0.031</td></tr><tr><td>1.654</td><td>1.181</td><td>0.032</td></tr><tr><td>1.662</td><td>1.18</td><td>0.032</td></tr><tr><td>1.67</td><td>1.172</td><td>0.033</td></tr><tr><td>1.678</td><td>1.171</td><td>0.033</td></tr><tr><td>1.686</td><td>1.153</td><td>0.032</td></tr><tr><td>1.695</td><td>1.143</td><td>0.031</td></tr></tbody></table>
<table><caption>**Table 4**
(Continued)</caption><thead><tr><th>Wavelength (μm)</th><th>Flux (10<sup>-15</sup> × W m<sup>-2</sup>/μm)</th><th>Flux Error (10<sup>-15</sup> × W m<sup>-2</sup>/μm)</th></tr></thead><tbody><tr><td>1.703</td><td>1.126</td><td>0.03</td></tr><tr><td>1.711</td><td>1.084</td><td>0.029</td></tr><tr><td>1.719</td><td>1.05</td><td>0.028</td></tr><tr><td>1.727</td><td>1.007</td><td>0.027</td></tr><tr><td>1.735</td><td>0.949</td><td>0.025</td></tr><tr><td>1.743</td><td>0.878</td><td>0.023</td></tr><tr><td>1.751</td><td>0.818</td><td>0.022</td></tr><tr><td>1.758</td><td>0.759</td><td>0.021</td></tr><tr><td>1.765</td><td>0.693</td><td>0.019</td></tr><tr><td>1.772</td><td>0.63</td><td>0.017</td></tr><tr><td>1.777</td><td>0.57</td><td>0.018</td></tr><tr><td>1.781</td><td>0.525</td><td>0.015</td></tr><tr><th colspan="3">GPI (K<sub>1</sub>)</th></tr><tr><td>1.892</td><td>0.547</td><td>0.063</td></tr><tr><td>1.898</td><td>0.627</td><td>0.054</td></tr><tr><td>1.905</td><td>0.608</td><td>0.081</td></tr><tr><td>1.907</td><td>0.671</td><td>0.032</td></tr><tr><td>1.916</td><td>0.654</td><td>0.026</td></tr><tr><td>1.924</td><td>0.624</td><td>0.02</td></tr><tr><td>1.932</td><td>0.588</td><td>0.019</td></tr><tr><td>1.941</td><td>0.567</td><td>0.019</td></tr><tr><td>1.95</td><td>0.559</td><td>0.018</td></tr><tr><td>1.96</td><td>0.561</td><td>0.019</td></tr><tr><td>1.969</td><td>0.582</td><td>0.02</td></tr><tr><td>1.977</td><td>0.582</td><td>0.018</td></tr><tr><td>1.985</td><td>0.561</td><td>0.018</td></tr><tr><td>1.993</td><td>0.564</td><td>0.018</td></tr><tr><td>2.003</td><td>0.633</td><td>0.018</td></tr><tr><td>2.016</td><td>0.694</td><td>0.021</td></tr><tr><td>2.025</td><td>0.693</td><td>0.022</td></tr><tr><td>2.033</td><td>0.692</td><td>0.022</td></tr><tr><td>2.041</td><td>0.697</td><td>0.022</td></tr><tr><td>2.049</td><td>0.727</td><td>0.022</td></tr><tr><td>2.059</td><td>0.759</td><td>0.023</td></tr><tr><td>2.069</td><td>0.762</td><td>0.025</td></tr><tr><td>2.077</td><td>0.772</td><td>0.027</td></tr><tr><td>2.086</td><td>0.744</td><td>0.022</td></tr><tr><td>2.094</td><td>0.752</td><td>0.023</td></tr><tr><td>2.103</td><td>0.764</td><td>0.028</td></tr><tr><td>2.111</td><td>0.79</td><td>0.024</td></tr><tr><td>2.12</td><td>0.79</td><td>0.028</td></tr><tr><td>2.129</td><td>0.778</td><td>0.028</td></tr><tr><td>2.138</td><td>0.787</td><td>0.024</td></tr><tr><td>2.147</td><td>0.771</td><td>0.025</td></tr><tr><td>2.155</td><td>0.749</td><td>0.022</td></tr><tr><td>2.163</td><td>0.728</td><td>0.021</td></tr><tr><td>2.17</td><td>0.695</td><td>0.022</td></tr><tr><td>2.176</td><td>0.662</td><td>0.023</td></tr><tr><td>2.181</td><td>0.615</td><td>0.022</td></tr><tr><td>2.183</td><td>0.471</td><td>0.025</td></tr></tbody></table>

<table><caption>**Table 5**
Extracted GPI Astrometry for HD 4747 B</caption><thead><tr><th>Time (BJD)</th><th>Separation (arcsec)</th><th>Position Angle (deg)</th></tr></thead><tbody><tr><td>2457380.5</td><td>0.5989 ± 0.002</td><td>183.9 ± 0.2</td></tr><tr><td>2457381.5</td><td>0.5984 ± 0.002</td><td>183.5 ± 0.2</td></tr></tbody></table>

## Appendix C
### Priors and Posteriors for Retrieval Parameters

We list the priors on our retrieved parameters in Table 6 and include joint posterior distributions of selected parameters from our baseline HRS and LRS retrievals in Figures 15 and 16.
<figure><img src="image_19.png" /><figcaption>**Figure 15.** Joint posterior distributions for the HRS retrieval of HD 4747 B. We omit the P–T profile parameters, which are better visualized by the P–T plot in Figure 6.</figcaption></figure>
<figure><img src="image_20.png" /><figcaption>**Figure 16.** Joint posterior distributions for the LRS retrieval of HD 4747 B. We omit the P–T profile parameters, which are better visualized by the P–T plot in Figure 6. The distributions for a few parameters are bimodal.</figcaption></figure>
<table><caption>**Table 6**
Priors of the HD 4747 B Retrieval</caption><thead><tr><th>Parameter</th><th>Prior</th><th>Parameter</th><th>Prior</th></tr></thead><tbody><tr><td>Mass (M<sub>Jup</sub>)</td><td>G(67.2, 1.8)</td><td>C/O</td><td>U(0.1, 1.6)</td></tr><tr><td>Radius (R<sub>Jup</sub>)</td><td>U(0.6, 1.2)</td><td>[Fe/H]</td><td>U(-1.5, 1.5)</td></tr><tr><td>T<sub>1</sub> (K)</td><td>U(0, T<sub>2</sub>)</td><td>log(P<sub>quench</sub>/bar)</td><td>U(-4, 3)</td></tr><tr><td>T<sub>2</sub> (K)</td><td>U(0, T<sub>3</sub>)</td><td>f<sub>sed</sub></td><td>U(0, 10)</td></tr><tr><td>T<sub>3</sub> (K)</td><td>U(0, T<sub>connect</sub>)<sup>(a)</sup></td><td>log(K<sub>zz</sub>/cm<sup>2</sup> s<sup>-1</sup>)</td><td>U(5, 13)</td></tr><tr><td>T<sub>int</sub> (K)</td><td>U(700, 2500)</td><td>σ<sub>g</sub></td><td>U(1.05, 3)</td></tr><tr><td>α</td><td>U(1, 2)</td><td>log(X&#773;<sub>MgSiO<sub>3</sub></sub>)<sup>(c)</sup></td><td>U(-2.3, 1)</td></tr><tr><td>log(δ)</td><td>P<sub>phot</sub> ∈ [10<sup>-3</sup>, 100]<sup>(b)</sup></td><td>log(X&#773;<sub>Fe</sub>)</td><td>U(-2.3, 1)</td></tr><tr><td colspan="4">Additional parameters for HRS</td></tr><tr><td>RV (km s<sup>-1</sup>)</td><td>U(-30, 30)</td><td>ν sin i (km s<sup>-1</sup>)</td><td>U(0, 50)</td></tr><tr><td>Error multiple</td><td>U(1, 4)</td><td>Flux scale (counts)</td><td>U(0, 200)</td></tr><tr><td colspan="4">Gaussian process parameters for LRS</td></tr><tr><td>log(f<sub>amp</sub>)</td><td>U(10<sup>-4</sup>, 1)</td><td>log(l) (μm)</td><td>U(10<sup>-3</sup>, 0.5)</td></tr><tr><td colspan="4">Mass fraction of molecules</td></tr><tr><td>log(MMR)</td><td>U(10<sup>-1</sup>, 10<sup>-7</sup>)</td><td></td><td></td></tr></tbody></table>

**Note.** *U* stands for a uniform distribution, with two numbers representing the lower and upper boundaries. *G* stands for a Gaussian distribution, with numbers representing the median and standard deviation. (a) and (b): these priors follow Mollière et al. (2020). *P*<sub>phot</sub> is the pressure where τ = 1, and *T*<sub>connect</sub> is the uppermost temperature of the “photospheric” layer, and is computed by setting τ = 0.1 in the Eddington Approximation (see Equations (1) and (2) in Mollière et al. 2020). This prior, along with those on *T*<sub>1</sub> and *T*<sub>2</sub> are used to prevent temperature inversions. (c) X&#773;<sub>MgSiO<sub>3</sub>/Fe</sub> represents the scaling factor for the cloud mass fraction, so that log(X&#773;<sub>MgSiO<sub>3</sub>/Fe</sub>) = 0 refers to a fraction equal to the equilibrium mass fraction. *f*<sub>sed</sub>, *K*<sub>zz</sub>, and *σ*<sub>g</sub> are parameters in the EddySed cloud model (Ackerman & Marley 2001). When fitting molecular abundances directly (e.g., in Section 5.4), we use the same mass fraction prior on all molecules included.

## ORCID iDs

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

Jason Wang https://orcid.org/0000-0003-0774-6502

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

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

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

Ji Wang https://orcid.org/0000-0002-4361-8885

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

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

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

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

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

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

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

J. Kent Wallace https://orcid.org/0000-0001-5299-6899

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

## References

Ackerman, A. S., & Marley, M. S. 2001, ApJ, 556, 872

Amarsi, A. M., Nissen, P. E., & Skúladóttir, Á. 2019, A&A, 630, A104

Apai, D., Radigan, J., Buenzli, E., et al. 2013, ApJ, 768, 121

Asplund, M., Grevesse, N., Sauval, A. J., & Scott, P. 2009, ARA&A, 47, 481

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

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

Bate, M. R. 2012, MNRAS, 419, 3115

Bate, M. R., Bonnell, I. A., & Bromm, V. 2002, MNRAS, 332, L65

Benneke, B., & Seager, S. 2013, ApJ, 778, 153

Beuzit, J.-L., Vigan, A., Mouillet, D., et al. 2019, A&A, 631, A155

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

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

Brandt, G. M., Michalik, D., Brandt, T. D., et al. 2021a, AJ, 162, 230

Brandt, T. D. 2021, ApJS, 254, 42

Brandt, T. D., Dupuy, T. J., & Bowler, B. P. 2019, AJ, 158, 140

Brandt, T. D., Dupuy, T. J., Li, Y., et al. 2021b, AJ, 162, 186

Brewer, J. M., Fischer, D. A., Valenti, J. A., & Piskunov, N. 2016, ApJS, 225, 32

Brown, A. G. A., Vallenari, A., Prusti, T., et al. 2021, A&A, 649, A1

Burgasser, A. J. 2014, in International Workshop on Stellar Spectral Libraries
ASI Conference Series, Vol.11 ed. H. P. Singh, P. Prugniel, & I. Vauglin
(Hyderabad: Astronomical Society of India), 7

Burningham, B., Faherty, J. K., Gonzales, E. C., et al. 2021, MNRAS, 506, 1944

Burningham, B., Marley, M. S., Line, M. R., et al. 2017, MNRAS, 470, 1177

Carmichael, T. W., Quinn, S. N., Mustill, A. J., et al. 2020, AJ, 160, 53

Chabrier, G., Baraffe, I., Leconte, J., Gallardo, J., & Barman, T. 2009, in AIP
Conf. Ser., 1094, 15th Cambridge Workshop on Cool Stars, Stellar Systems,
and the Sun, ed. E. Stempels (Melville, NY: AIP), 102

Claudi, R. U., Turatto, M., Gratton, R. G., et al. 2008, Proc. SPIE, 7014, 70143E

Crepp, J. R., Gonzales, E. J., Bechter, E. B., et al. 2016, ApJ, 831, 136

Crepp, J. R., Principe, D. A., Wolff, S., et al. 2018, ApJ, 853, 192

Crossfield, I. J. M., Biller, B., Schlieder, J. E., et al. 2014, Natur, 505, 654

Currie, T., Brandt, T. D., Uyama, T., et al. 2018, AJ, 156, 291

Cushing, M. C., Roellig, T. L., Marley, M. S., et al. 2006, ApJ, 648, 614

Cutri, R. M., Skrutskie, M. F., van Dyk, S., et al. 2003, VizieR On-line Data
Catalog: II/246

Czesla, S., Schröter, S., Schneider, C. P., et al. 2019, PyA: Python astronomy-related packages, Astrophysics Source Code Library, ascl:1906.010

De Rosa, R. J., Rameau, J., Patience, J., et al. 2016, ApJ, 824, 121

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

Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306

Gandhi, S., Brogi, M., & Webb, R. K. 2020, MNRAS, 498, 194

Gao, P., Marley, M. S., & Ackerman, A. S. 2018, ApJ, 855, 86

Gao, P., Thorngren, D. P., Lee, G. K. H., et al. 2020, NatAs, 4, 951

Gierasch, P. J., & Conrath, B. J. 1985, in Recent Advances in Planetary Meteorology, ed. G. E. Hunt (Cambridge: Cambridge Univ. Press), 121
Gonzales, E. C., Burningham, B., Faherty, J. K., et al. 2020, ApJ, 905, 46

Gonzales, E. C., Burningham, B., Faherty, J. K., et al. 2021, ApJ, 923, 19

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

Hargreaves, R. J., Gordon, I. E., Rey, M., et al. 2020, ApJS, 247, 55

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

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

Johnston, V. D., Cushing, M. C., & Noll, K. S. 2019, RNAAS, 3, 52

Karalidi, T., Marley, M., Fortney, J. J., et al. 2021, ApJ, 923, 269

Kaufer, A., Wolf, B., Andersen, J., & Pasquini, L. 1997, Msngr, 89, 1

Kitzmann, D., Heng, K., Oreshenko, M., et al. 2020, ApJ, 890, 174

Kolecki, J. R., Wang, J., Johnson, J. A., et al. 2021, AJ, 162, 125

Konopacky, Q. M., Ghez, A. M., Fabrycky, D. C., et al. 2012, ApJ, 750, 79

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

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

Line, M. R., Teske, J., Burningham, B., Fortney, J. J., & Marley, M. S. 2015, ApJ, 807, 183

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

Luna, J. L., & Morley, C. V. 2021, ApJ, 920, 146

Madhusudhan, N. 2012, ApJ, 758, 36

Madhusudhan, N., Amin, M. A., & Kennedy, G. M. 2014, ApJ, 794, L12

Mamajek, E. E., & Hillenbrand, L. A. 2008, ApJ, 687, 1264

Manjavacas, E., Karalidi, T., Vos, J. M., Biller, B. A., & Lew, B. W. P. 2021, AJ, 162, 179

Marley, M., & Robinson, T. 2015, ARA&A, 53, 279

Marley, M. S., Saumon, D., Visscher, C., et al. 2021, ApJ, 920, 85

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

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

Miles, B. E., Skemer, A. J. I., Morley, C. V., et al. 2020, AJ, 160, 63

Min, M., Hovenier, J. W., & de Koter, A. 2005, A&A, 432, 909

Mollière, P., Wardenier, J. P., van Boekel, R., et al. 2019, A&A, 627, A67

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

Mordasini, C., van Boekel, R., Mollière, P., Henning, T., & Benneke, B. 2016, ApJ, 832, 41

Moses, J. I., Madhusudhan, N., Visscher, C., & Freedman, R. S. 2013, ApJ, 763, 25

Mukherjee, S., Batalha, N. E., Fortney, J. J., & Marley, M. S. 2022, arXiv:2208.07836

Noll, K. S., Geballe, T. R., Leggett, S. K., & Marley, M. S. 2000, ApJ, 541, L75

Öberg, K. I., & Bergin, E. A. 2016, ApJ, 831, L19

Öberg, K. I., Murray-Clay, R., & Bergin, E. A. 2011, ApJL, 743, L16

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

Padoan, P., & Nordlund, Å. 2004, ApJ, 617, 559

Pelletier, S., Benneke, B., Darveau-Bernier, A., et al. 2021, AJ, 162, 73

Peretti, S., Ségransan, D., Lavie, B., et al. 2019, A&A, 631, A107

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

Piso, A.-M. A., Öberg, K. I., Birnstiel, T., & Murray-Clay, R. A. 2015, ApJ, 815, 109

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

Pueyo, L. 2016, ApJ, 824, 117

Riello, M., De Angeli, F., Evans, D. W., et al. 2021, A&A, 649, A3

Rosenthal, L. J., Fulton, B. J., Hirsch, L. A., et al. 2021, ApJS, 255, 8

Ruffio, J.-B. 2019, PhD Thesis, Stanford Univ.

Ruffio, J.-B., Konopacky, Q. M., Barman, T., et al. 2021, AJ, 162, 290

Samland, M., Mollière, P., Bonnefoy, M., et al. 2017, A&A, 603, A57

Saumon, D., & Marley, M. S. 2008, ApJ, 689, 1327

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

Smith, M. D. 1998, Icar, 132, 176

Sneden, C. 1973, ApJ, 184, 839

Soummer, R., Pueyo, L., & Larkin, J. 2012, ApJ, 755, L28

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

Stamatellos, D., Hubber, D. A., & Whitworth, A. P. 2007, MNRAS, 382, L30

Stelter, R. D., Skemer, A. J., Sallum, S., et al. 2020, Proc. SPIE, 11447, 1346

Tan, X. 2022, MNRAS, 511, 4861

Tan, X., & Showman, A. P. 2021, MNRAS, 502, 678

Thies, I., Kroupa, P., Goodwin, S. P., Stamatellos, D., & Whitworth, A. P. 2010, ApJ, 717, 577

Toon, O. B., Turco, R. P., Westphal, D., Malone, R., & Liu, M. 1988, JAtS, 45, 2123

Turco, R. P., Hamill, P., Toon, O. B., Whitten, R. C., & Kiang, C. S. 1979, JAtS, 36, 699

Vigan, A. 2020, vlt-sphere: Automatic VLT/SPHERE data reduction and analysis, Astrophysics Source Code Library, ascl:2009.002

Vos, J. M., Faherty, J. K., Gagné, J., et al. 2022, ApJ, 924, 68

Vousden, W. D., Farr, W. M., & Mandel, I. 2016, MNRAS, 455, 1919

Wang, J., Kolecki, J. R., Ruffio, J.-B., et al. 2022, AJ, 163, 189

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

Wang, J., Mawet, D., Ruane, G., Hu, R., & Benneke, B. 2017, AJ, 153, 183

Wang, J., Wang, J. J., Ma, B., et al. 2020a, AJ, 160, 150

Wang, J. J., Kulikauskas, M., & Blunt, S. 2021a, whereistheplanet: Predicting positions of directly imaged companions, Astrophysics Source Code Library, ascl:2101.003

Wang, J. J., Ruffio, J.-B., De Rosa, R. J., et al. 2015, pyKLIP: PSF Subtraction for Exoplanets and Disks, Astrophysics Source Code Library, ascl:1506.001

Wang, J. J., Ruffio, J.-B., Morris, E., et al. 2021b, AJ, 162, 148

Wang, J. J., Vigan, A., Lacour, S., et al. 2021c, AJ, 161, 148

Xuan, J. W., & Wyatt, M. C. 2020, MNRAS, 497, 2096

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

Zahnle, K. J., & Marley, M. S. 2014, ApJ, 797, 41

Zhang, X., & Showman, A. P. 2018, ApJ, 866, 1

Zhang, Y., Snellen, I. A. G., & Mollière, P. 2021, A&A, 656, A76

Zhou, Y., Apai, D., Metchev, S., et al. 2018, AJ, 155, 132