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

In this section, we summarize relevant properties of the host star. HD 4747 is a main-sequence, solar-type star located ≈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 ± 0.02), which Peretti et al. (2019) used to derive an age of 2.3 ± 1.4 Gyr from the age-log R_HK calibration of Mamajek & Hillenbrand (2008). This agrees with the gyro-chronological age estimate of 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 ± 0.05 from Peretti et al. (2019) and [Fe/H] = −0.22 ± 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 ± 0.08. The iron abundance is [Fe/H] = −0.30 ± 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 ± 0.06 and [O/H] = − 0.02 ± 0.04 for the host star. [C/H] is defined as - , 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.

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 orvara 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 (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 orvara, 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 ± 1.8M_Jup, which is consistent with previous values, but more precise. We checked the orvara 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 ± 2.0M_Jup, consistent with the orvara 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 orvara fit for the spectral analysis in this paper.

3.1.1. KPIC Observations

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. ¹⁹ 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_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.

Zoom In Zoom Out Reset image size

3.2.1. Gemini Planet Imager IFS

3.2.2. SPHERE IFS

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.

Zoom In Zoom Out Reset image size

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 (≈06), we find the speckles are roughly the same brightness as our companion (K_s ≈ 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 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 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_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.

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] ²³ 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₂O, CH₄, NH₃, and CO₂, 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⁶, much higher than the resolution of our HRS resolution (R ≈ 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⁶) 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

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).

Table 2. Spectral Retrievals Carried Out on HD 4747 B

Data/Reference	Cloud Model	C/O	[C/H]	Radius (RJup)	log(g)	Teff (K)	B
HRS (KPIC)	EddySed (MgSiO3, am)	0.66 ± 0.04					1.0
HRS	EddySed (MgSiO3, cd)	0.67 ± 0.04		0.90 ± 0.19			1.15
HRS	Clear						0.61
HRS	Clear (chemical equilibrium)	0.60 ± 0.02					1.6 × 10−3
LRS (GPI+SPHERE)	EddySed (MgSiO3, am)						1.0
LRS	EddySed (MgSiO3, cd)			0.77 ± 0.04		1443 ± 28	0.69
LRS	EddySed (MgSiO3 + Fe, am)			0.73 ± 0.03			1.54
LRS	EddySed (MgSiO3 + Fe, cd)			0.75 ± 0.03			2.65
LRS	Clear			1.10 ± 0.04	5.12 ± 0.03	1262 ± 16	7.0 × 10−26
Peretti et al. (2019)	Cloudy retrieval a			0.85 ± 0.03 b	5.40 ± 0.03	1350 ± 50	⋯
Crepp et al. (2018)	Cloudy grid a )	⋯	⋯	⋯			⋯

Notes. The rightmost column lists the Bayes factor (B) for each retrieval, with the EddySed (MgSiO₃, 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.

^a 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). ^b Peretti et al. (2019) placed a Gaussian prior of 1.0 ± 0.1R_Jup on the radius.

Download table as:  ASCIITypeset image

For HD 4747 B, we consider models with two different cloud species (MgSiO₃ and Fe) and properties (amorphous or crystalline particles). We choose to focus on MgSiO₃ and Fe for several reasons. First, the condensation curves of these two species intersect the thermal profile of a T_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₃ is expected to be the most important cloud species for substellar objects with T_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₃ absorption feature at ∼10 μm in field brown dwarfs (Cushing et al. 2006; Luna & Morley 2021), and specifically amorphous MgSiO₃ (Burningham et al. 2021). Although MgSiO₃ 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₃ modeled with Mie theory (abbreviated MgSiO₃, “am”) for the clouds. In addition, we also consider models with MgSiO₃, “cd,” which assumes crystalline cloud particles modeled with DHS, as well as models with two cloud species (MgSiO₃ + Fe) for the LRS.

4.2.3. Methane Opacities

4.2.4. Additional Fit Parameters

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 ± 1.8M_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_Jup. When including a quench pressure, we use a log-uniform prior from 10⁻⁴ to 10³, 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

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.

Zoom In Zoom Out Reset image size

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 ± 0.1 km s⁻¹ 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 ± 0.4 km s⁻¹ 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 ΔK_s = 8.3 ± 0.3 mag, which is within 3σ of the photometric ΔK_s = 9.05 ± 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.

Zoom In Zoom Out Reset image size

Figure 4 also shows the projected spin rate 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_quench in Figure 4, which indicates that the chemical reaction timescale becomes longer than the vertical mixing timescale at pressures lower than P_quench. Thus, disequilibrium chemistry is clearly affecting the atmosphere (see Section 5.3 for details).

We compute T_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_eff using the Stefan–Boltzmann law. When computing T_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_eff from HRS have broad distributions, which reflect the relatively weak luminosity constraints from the HRS (log(L_bol/L_⊙) = ). This is because the HRS is not flux-calibrated and we remove the continuum in our fits. Comparing to values of radius and T_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σ–2σ 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σ. This results in our retrieved C/O for the companion being higher by about 2σ compared to the stellar value. Here and elsewhere in the paper, we compute the “σ 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.

Zoom In Zoom Out Reset image size

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, , 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 ≈10⁻² bars in the cores 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 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₃.

Zoom In Zoom Out Reset image size

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^–1, 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.

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 ≈ 380, in favor of the quenched retrieval. From the Bayes factor interpretation of Benneke & Seager (2013), this represents a detection of quenching at ≈3.9σ 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₄ = . If the atmosphere was in chemical equilibrium, we repeat our calculation and find that we would expect CO/CH₄ = , 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.

Zoom In Zoom Out Reset image size

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.

Zoom In Zoom Out Reset image size

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/CH4 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.

Table 3. Free Retrievals Carried Out on HD 4747 B for Validating the CH₄ Detection in the KPIC HRS (Section 5.4)

Molecules	log(CO)	log(H2O)	log(CH4)	B
CO, H2O			⋯	1
CO, H2O, CH4			−4.82 ± 0.23	84

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

Download table as:  ASCIITypeset image

We find that the data strongly prefer the model with CH₄, with a Bayes factor of 84 (3.4σ significance; Benneke & Seager 2013). 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).

Zoom In Zoom Out Reset image size

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 (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 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.

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.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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 ≈ 7 × 10⁻²⁶; 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.

Zoom In Zoom Out Reset image size

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. 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 degeneracy seen here may be related to issues with the LRS (see Section 6.3).

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 ± 2 Gyr and m = 67.2 ± 1.8 M_Jup), we interpolate the COND evolutionary models (Baraffe et al. 2003) to find a model-predicted radius of , and a predicted 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 ± 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 K, and a radius of . Compared to prior expectations, however, the radius retrieved is too small by ≈1.5σ while T_eff is too high by ≈2σ. From substellar evolutionary models, the minimum possible radius of a brown dwarf should be ≈0.74 R_Jup, which is imposed by electron degeneracy pressure (Chabrier et al. 2009). We find that T_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₃ 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₃) 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₃ only). Furthermore, the MgSiO3, 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.

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₂O and CH₄ opacities. When we compute the CO abundances from our LRS retrievals, we find that they are not well constrained, with 1σ 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₄ and H₂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.

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⁻² 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 (≈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 and assuming a radius of 0.8R_Jup, we would expect a 5 or 7 hr rotation period if i is equal to the orbital inclination or i = 90°. 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).

Our HRS retrievals indicate that the ratio of CO/CH4 (VMR) is ≈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₄ reaction from Zahnle & Marley (2014) with the mixing timescale τ_mix = L²/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 ≈ 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 = α H, where α is a scaling factor, and (μ: 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 α = 0.1, we find K_zz = 5 × 10⁸ − 1 × 10¹² cm² s⁻¹ (1σ). On the other hand, if α = 1, we obtain K_zz = 5 × 10¹⁰ − 1 × 10¹⁴ cm² s⁻¹.

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_eff between 250 and 750 K, much colder than HD 4747 B. In terms of objects with T_eff ≳ 1000 K, Barman et al. (2015) reported a detection of CH₄ in HR 8799b (T_eff ∼ 1000 K) with Keck/OSIRIS data, which they used to estimate K_zz between 10⁶ and 10⁸ cm² s⁻¹. However, this CH₄ 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₄ signal strengths found by different analyses. Ruffio et al. (2021) pointed out that if the CH₄ 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_eff ∼ 1100 K) from petitRADTRANS retrievals, which could similarly be converted to a K_zz constraint. In summary, our finding HD 4747 B, which is ∼300–400 K hotter than HR 8799 b and e 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 ≈10⁹ cm² s⁻¹. Depending on α, our retrieved K_zz either exceeds this upper limit by ≳2σ (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₄ VMR is consistent with K_zz ∼ 10⁸–10¹² in these models (with the assumption that L = H). These values of K_zz are roughly consistent with our estimate based on P_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.

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_Jup (“free-mass”). For the HRS free-mass retrieval, we find that all parameters change by less than 1σ compared to the mass-prior retrieval. The mass itself shows a broad distribution (33–76 M_Jup at 1σ) 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σ and 2σ with the mass-prior retrieval. Furthermore, the mass retrieved by the LRS is , which agrees within about 1σ 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 ± 0.03R_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.

We retrieve [C/H] and [O/H] values that are 1σ 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 ± 0.04 is higher by approximately 2σ than the stellar C/O =0.48 ± 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 ∼0.15 level in C/O. Using KPIC HRS from 2.23–2.33 μm, they found ≈1σ–1.5σ 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σ 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σ 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σ 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 MgSiO3 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.