The Carbon-to-oxygen Ratio in Cool Brown Dwarfs and Giant Exoplanets. I. The Benchmark Late-T Dwarfs GJ 570D, HD 3651B, and Ross 458C

Measurements of the C/O ratio in brown dwarfs are lacking, in part due to past models adopting solar C/O only. We have expanded the ATMO 2020 atmosphere model grid to include nonsolar metallicities and C/O ratios in the T dwarf regime. We change the C/O ratio by altering either the carbon or oxygen elemental abundances, and we find that nonsolar abundances of these elements can be distinguished based on the shapes of the H and K bands. We compare these new models with medium-resolution (R ≈ 1700), near-infrared (0.8–2.4 μm) Gemini Near-Infrared Spectrograph (GNIRS) spectra of three benchmark late-T dwarfs, GJ 570D, HD 3651B, and Ross 458C. We find solar C/O ratios and best-fitting parameters (T eff, log(g) , Z) broadly consistent with other analyses in the literature based on low-resolution (R ∼ 100) data. The model-data discrepancies in the near-infrared spectra are consistent across all three objects. These discrepancies are alleviated when fitting the Y, J, H, and K bands individually, but the resulting best-fit parameters are inconsistent and disagree with the results from the full spectrum. By examining the model atmosphere properties we find this is due to the interplay of gravity and metallicity on H2–H2 collisionally induced absorption. We therefore conclude that there are no significant issues with the molecular opacity tables used in the models at this spectral resolution. Instead, deficiencies are more likely to lie in the model assumptions regarding the thermal structures. Finally, we find a discrepancy between the GNIRS, SpeX, and other near-infrared spectra in the literature of Ross 458C, indicating potential spectroscopic variability.


INTRODUCTION
Spectroscopy of cool brown dwarfs and giant exoplanets allows us to study the physics and chemistry of their atmospheres.Measuring spectroscopic signatures of composition has long been cited as a method of inferring the formation mechanism of planetary objects (e.g.Gaidos 2000;Lodders 2004;Kuchner & Seager 2005).Of particular importance is the C/O ratio, which in exoplanets can depart from that of their host stars due to their formation and evolution within a protoplanetary disk, thus providing an observational link between present-day atmospheric composition and formation pathway.Compared to the solar value of C/O ≈ 0.55 (Asplund et al. 2009), the C/O ratio of exoplanets can become elevated as high as C/O ≈ 1 when forming beyond the snowlines of carbonand oxygen-rich ices ( Öberg et al. 2011;Madhusudhan et al. 2011;Mollière et al. 2022).Stellar abundances of C and O establish the initial conditions for disk chemistry in planet-formation models (Bond et al. 2010), thus motivating studies of the stellar C/O ratio in the solar neighborhood.In solar-type and low-mass stars (including planet-hosting stars), the C/O ratio has been found to vary between ≈ 0.3 − 0.8 (Nissen et al. 2014;Nakajima & Sorahana 2016;Brewer & Fischer 2016;Suárez-Andrés et al. 2018).Such C/O determinations have been obtained from high resolution spectroscopy (see Nissen & Gustafsson (2018) for a review), relying on weak individual absorption lines of C and O.
In contrast, the dominant absorbers in cool brown dwarf atmospheres are water, methane and carbon monoxide, which carve out broad absorption bands across the infrared and account for almost all of the atmospheric C and O. Brown dwarfs therefore offer a unique diagnostic of the C/O ratio in the solar neighbourhood, arising from a different atmospheric chemistry than higher mass stars (Fortney 2012).Furthermore, wide-separation (> 100 au) brown dwarf companions to solar-type stars provide crucial benchmarks for such investigations thanks to the age and composition constraints derived from their host stars (e.g.Crepp et al. 2018;Rickman et al. 2020;Zhang et al. 2020).Comparing the measured C/O ratio of a brown dwarf companion to that of the host star could shed light on the companion's formation pathway, e.g.planet-like formation in a disk or star-like formation from cloud fragmentation.However, measurements of the C/O ratio in brown dwarfs, both companions and free-floating objects, are lacking in part due to past models encoding solar C/O into the model assumptions.
Over the last few years, new grids of coupled atmosphere and evolution models have been developed (Phillips et al. 2020;Marley et al. 2021).These models include the latest molecular opacities and treatments for atmospheric chemistry, benefitting from a wide range of developments since the last major set of models over a decade ago.The improvements to the chemistry are particularly important for accurately determining the C/O ratio and non-solar abundances.First, the addition of nonequilibrium chemistry that is self-consistent with the pressure-temperature profile is vital (Drummond et al. 2016), as mixing processes can drive the abundances of carbon-and oxygen bearing molecules away from their equilibrium values.Second, the addition of rainout chemistry, which models the depletion of elements due to the settling/sinking of condensates in an atmosphere (Burrows & Sharp 1999), is also important since the formation of condensates can deplete up to 20% of the oxygen atoms from the observable atmosphere (Lodders 2010), impacting the observed C/O ratio.This new generation of atmosphere models allows us to more reliably investigate the spectroscopic signatures of non-solar chemical abundances in cool brown dwarfs.
The majority of brown dwarf spectroscopy has been conducted at low spectral resolution (R < 1000) (Leggett et al. 2001;Sorahana & Yamamura 2014;Burgasser 2014;Schneider et al. 2015).Many observations have been taken with the SpeX spectrograph (Rayner et al. 2003) on the NASA Infrared Telescope Facility (IRTF), which provides R ≈ 75 − 120 near-infrared (0.8 − 2.5 µm) spectra in its prism mode.These observations have led to large compilations such as the SpeX Prism Library (Burgasser 2014), which contains thousands of spectra of ultracool M, L, and T dwarfs.The 0.8 − 2.5 µm wavelength range provided by SpeX captures the majority of the spectral energy distribution of M, L and T-type ultracool dwarfs.Mid-late T dwarfs are perhaps the most promising objects to constrain C/O ratios, with the clouds and/or thermochemical instabilities that complicate the modeling and studies of earlier spectral types having sunk below the photosphere.The numerous SpeX prism spectra available in the literature have led to many T dwarf model-data comparisons, using fully consistent one-dimensional radiative-convective atmosphere models (Oreshenko et al. 2020;Zhang et al. 2021b,c) and data-driven atmospheric retrieval frameworks (Line et al. 2015(Line et al. , 2017;;Zalesky et al. 2019;Kitzmann et al. 2020;Zalesky et al. 2022;Whiteford et al. 2023) to interpret the atmospheric properties.However, at low spectral resolution, molecular absorption lines are not resolved, so high resolutions are required to test the opacities and line lists used in the atmosphere models.A small sample of T dwarfs have been observed at medium-resolution with the Folded-Port Infrared Echelette (FIRE) spectrograph (R ≈ 6000; 0.85 − 2.5 µm) at the Magellan telescopes (Burgasser et al. 2011;Bochanski et al. 2011;Hood et al. 2023) and the Near-Infrared Integral Field spectrograph (NIFS) (R ≈ 5000; 1.0 − 2.4 µm) at the Gemini North telescope (Canty et al. 2015).T dwarfs have been observed at high resolution with IGRINS (R ≈ 45, 000; 1.46 − 2.48 µm; Tannock et al. 2022) and with NIRSPEC (R ≈ 35, 000; 2.29 − 2.49 µm; Xuan et al. 2022) and (R ≈ 20000;1.165 − 323 µm;McLean et al. 2007).However, these studies typically focused on a single object and/or had a limited wavelength range.There is thus a need for higher-resolution spectra than provided by SpeX near-infrared for a systematic sample of cool brown dwarfs.
To address this, we are conducting a spectroscopic survey of mid/late T dwarfs obtaining mediumresolution (R ≈ 1700), high signal-to-noise near-infrared spectra (0.8−2.4 µm) using the Gemini Near-Infrared Spectrograph (GNIRS).In this work, we present GNIRS spectra for the benchmark brown dwarfs GJ 570D, HD 3651B and Ross 458C, which are all wide-separation companions to higher mass stars.To analyse and interpret these spectra, we have expanded the ATMO 2020 model atmosphere grid (Phillips et al. 2020) to include non-solar abundances, through changing the bulk metallicity and C/O ratio of the model atmospheres in the mid/late T dwarf regime (T eff = 700 − 900 K).The paper is organised as follows.In Section 2.1 we outline the details of our GNIRS observations and our data reduction methods.In Section 2.2 we discuss the details of the ATMO model atmosphere grid.The results of fitting the ATMO models to the GNIRS spectra of GJ 570D, HD 3651B and Ross 458C are presented in Section 3. We then analyse and discuss the discrepancies between the models and data in Section 4, before summarizing our work in Section 5.

Gemini/GNIRS spectroscopy
We construct a sample of cool brown dwarfs from the UltracoolSheet (Best et al. 2020) which contains astrometry, photometry, spectral classifications and multiplicities for 3000+ ultracool dwarfs and imaged exoplanets from multiple sources and surveys (Dupuy & Liu 2012;Dupuy & Kraus 2013;Liu et al. 2016;Best et al. 2018Best et al. , 2021)).We select spectral types T6 and later, which correspond to T eff ≲ 1000 K (Filippazzo et al. 2015), since the clouds that complicate the modeling of earlier spectral types are believed to have sunk below the photosphere for late-T dwarfs.We then apply a cut in apparent magnitude of J < 16 to select the brightest objects, relaxing this to J < 17 for objects of interest such as benchmark companions and young moving group members.
This paper focuses on the benchmark brown dwarfs GJ 570D, HD 3651B and Ross 458C, with the wider sample being the subject of future papers.We observed these objects at the Gemini-North 8.1-m telescope with the facility near-infrared spectrograph GNIRS (Elias et al. 2006), using the short blue camera's cross-dispersed (SXD) mode with the 32 l/mm grating and 0.3" slit aligned to the parallactic angle, providing R ≈ 1700.We took four 189 s exposures of GJ 570D, eight 280 s exposures of HD 3651B, and sixteen 326 s exposures of Ross 458C, dithered in an ABBA pattern along the slit.We observed the standard stars HIP 73867 (V = 8.6 mag), HIP 6193 (V = 4.7 mag) and HIP 68868 (V = 8.6 mag) for telluric correction of GJ 570D, HD 3651B, and Ross 458C, respectively, using the same instrument setup and dither pattern.The observations are summarized in Table 1.
We reduce the Gemini/GNIRS spectra using the open-source Python package PypeIt (Prochaska et al. 2020a,b;Bochanski et al. 2009;Bernstein et al. 2015).PypeIt is a semi-automated spectroscopic data reduction pipeline that can be applied to slit-imaging spectrographs and supports longslit, multislit and cross-dispersed echelle spectra.To reduce spectra, PypeIt requires calibration flat frames for order edge tracing and pixel efficiency correction, as well as arc frames for wavelength calibration.We use 6 nighttime flats taken with a quartz-halogen lamp for the flat frames.PypeIt determines the wavelength solution in the infrared from OH sky lines, which is preferable to using the images of the Ar lamp due to flexure.Therefore the science frames are used to perform the wavelength calibration and to trace the tilt in the wavelength solution across the slit.However, for short exposures (< 60 s) the OH sky lines are too weak, and thus for the telluric standard star PypeIt uses the tilt and wavelength solution from the nearest-in-time science image.The OH sky lines provide robust wavelength solutions, with the typical RMS for each order being ≤ 0.1 pixels (6 × 10 −5 µm).PypeIt achieves Poisson limited sky-subtraction (Prochaska et al. 2020a) and then performs an optimal extraction to generate 1D science spectra.
Flux calibration and telluric correction are performed using an infrared sensitivity function.The sensitivity function is generated from the observations of the standard star, by comparing the photon flux of the standard with its absolute, previously calibrated flux.In the case of the infrared observations conducted in this work, we use the Vega spectrum to generate the sensitivity function.The telluric model is similarly derived by fitting to a grid of model atmosphere spectra, with a range of atmospheric conditions specific to different observatories, including Mauna Kea.The sensitivity function can then be applied to the science observations which converts the measured photons per second per Å to flux units (erg s −1 Å−1 cm −2 ), followed by telluric correction.
After extraction, PypeIt rebins the spectra onto a common logarithmic wavelength grid, which facilitates the coaddition of multiple exposures and the combination of overlapping echelle orders.Multiple science exposures of the same object are combined for each echelle order by taking their average, weighted by the square of the median S/N ratio.To do this, a scale factor which normalizes the spectra to a common intensity, is calculated from the frame with the largest signal.Finally, to combine the edges of overlapping echelle orders, PypeIt averages the data weighting by the inverse variance of each pixel.For more information on the algorithms used in the PypeIt data reduction pipeline, we refer the reader to Bochanski et al. (2009), Bernstein et al. (2015), and the PypeIt online documentation 1 .
The final flux-calibrated, coadded, and telluric-corrected spectra of GJ 570D, HD 3651B and Ross 458C are shown in Figures 1, 2 and 3, respectively.We flux-calibrate the spectra to the MKO H-band magnitude of all objects following Zhang et al. (2021b).

Atmosphere model
To analyse the GNIRS spectra we use the recently published ATMO 2020 atmosphere models (Phillips et al. 2020).Briefly, the ATMO code computes the pressure-temperature structure of a one-dimensional atmosphere by solving for radiative-convective and hydrostatic equilibrium in each model level.We adopt the same model setup described in detail in Phillips et al. (2020) and references therein, which includes descriptions of the opacity database, radiative transfer, and chemistry schemes used by the ATMO code.Notably, the ATMO models include the latest line shapes for the potassium resonance doublet (Allard et al. 2016) and self-consistent non-equilibrium chemistry due vertical mixing.The initial tranche of publically available ATMO 2020 models were generated at solar metallicity.We have developed new models with non-solar compositions, through changing the bulk metallicity and C/O ratio of the atmosphere2 .The grid spans T eff = 700 − 900 K in steps of 100 K, log(g) = 4.5 − 5.5 dex in steps of 0.5 dex, Z = −0.5 − 0.5 dex in steps of 0.5 dex, and C/O = 0.35, 0.55(solar), and 1.0.In total there are 162 models.The C/O ratio can be changed by altering either the carbon or oxygen elemental abundances.In this work, we vary either the carbon or oxygen abundance through scaling factors X C and X O , respectively, while keeping the other elemental abundances constant.The bulk metallicity is changed by scaling all elements other than H and He by a constant factor.In our models we adopt an eddy diffusion coefficient of log(K zz ) = 4, which is held constant throughout both the radiative and convective regions of the atmosphere.We note that the value of K zz is uncertain by multiple orders of magnitude, particularly in the radiative region (Mukherjee et al. 2022).We also compute some models with log(K zz ) = 6 and log(K zz ) = 8 to examine this parameter's effect on the model spectra.However in the grid-fitting analyses presented in Section 3, we use only the log(K zz ) = 4 models.We choose a value of log(K zz ) = 4 to be consistent with literature values that have been found to reasonably match the observed colors of late T and Y dwarfs (Leggett et al. 2017).
We use these models to interpret the observed GNIRS spectra.To do this, we first degrade the model spectra to the GNIRS resolution of R = 1700.This is done by convolving the model spectrum with a Gaussian kernel, and then binning the spectrum onto the same wavelength grid as the observations.We then bi-linearly interpolate between the model grid points, such that the model sub-grid steps are ∆T eff = 20 K, ∆ log(g) = 0.1 dex, ∆Z = 0.1 dex, ∆X O = 0.1 (similarly for ∆X C ), between the minimum and maxiumum value of each parameter giving 34, 607 models in total.Each model in the grid is then scaled by the geometric dilution factor α = R 2 /D 2 , where R is the radius of the object and D is the distance to the object.The dilution factor can be obtained by setting ∂χ 2 ∂α = 0, where χ 2 is the chi-squared statistic, defined as Here O i and M i are the observed and model flux densities, σ i is the uncertainties on the observed flux densities, and N is the number of data points.The maximum-likelihood value of α is then The best-fit model is then found by computing the χ 2 statistic for all the models and then choosing the minimum value.We report our results as the reduced χ 2 , defined as χ 2 ν = (1/ν)χ 2 where ν is the number of degrees-of-freedom.We note that the values are large compared to the standard expectation of a good model fit (χ 2 ν ≈ 1).This is a well-known issue in the literature, given that χ 2 values are influenced not only by the measurement uncertainties, but also the systematic deficiencies in the model assumptions.
We adopt uncertainties on the best-fit parameters of half the model grid spacing (before interpolation).This is approximately consistent with the findings of Zhang et al. (2021b), who found uncertainties typically in the range of 1/3 − 1/2 when using the Bayesian inference tool Starfish combined with the Sonora-Bobcat (Marley et al. 2021) models to conduct an analysis of low-resolution T dwarf spectra.Uncertainties of half the model grid spacing have been commonly adopted when using a traditional χ 2 based model fitting analysis (Cushing et al. 2008;Zhang et al. 2020), since the models are known to have systematic errors and thus Equation 1 will not follow a true χ 2 distribution.

RESULTS
First we discuss the results of our atmosphere modeling and explore the spectral signatures of non-solar metallicities and C/O ratios in Section 3.1.We then present the results of fitting these atmosphere models to our GNIRS spectra of GJ 570D, HD 3651B and Ross 458C in Sections 3.2, 3.3 and 3.4, respectively.

Non-solar ATMO models
Our new non-solar ATMO models allow us to explore the effect that changing the C/O ratio of the atmosphere has on the near-infrared spectrum.The dominant near-infrared absorbers H 2 O, CH 4 and CO vary as a function of C/O ratio, as shown in Figure 4. CH 4 and CO are the dominant carbon-bearing species at low (< 1300 K) and high (> 1300 K) temperature, respectively.Oxygen is predominantly contained within H 2 O at low temperature, and in H 2 O and CO at higher temperatures.It can be clearly seen that changing the C/O ratio, the abundances of these key molecules varies depending on whether the abundance of C or O is changed.
When varying the oxygen abundance of the atmosphere, at lower temperatures CH 4 is independent of C/O, and H 2 O and CO both decrease with increasing C/O as the number of oxygen atoms decreases.At higher temperatures the abundance of CO is limited by the fixed carbon abundance, and is independent of C/O for C/O < 1.For C/O > 1, the decreasing oxygen abundance limits the CO abundance.The abundance of H 2 O increases as C/O decreases due to the increasing number of oxygen atoms.
When varying the carbon abundance of the atmosphere, at low temperatures H 2 O is independent of C/O and CH 4 increases with increasing C/O as the number of available carbon atoms increases.At higher temperatures, the abundance of CO increases with increasing C/O for C/O < 1 as the number of available carbon atoms increases, but become independent of C/O for C/O > 1 due to the fixed oxygen abundance.The H 2 O abundance also varies when varying the carbon abundance at high temperatures, since it depends on the available oxygen atoms and thus the abundance of CO.For C/O > 1, when varying both carbon and oxygen abundances, trace species of C 2 H 2 and HCN become more abundant than H 2 O, due to the increasing number of free carbon atoms.These trace species have spectral features in the mid-infrared, at ∼ 3, ∼ 7.5 and ∼ 14 µm.
Figure 5 shows the contribution function and abundance-weighted absorption cross-sections of the main atomic and molecular opacity sources in a typical T dwarf atmosphere from the ATMO 2020 model grid (T eff = 800 K, log(g) = 4.5, assuming chemical equilibrium with solar elemental abundances).Combined with the variation of the molecular abundances with C/O ratio shown in Figure 4, we can interpret the impact the model grid parameters have on the observable near-infrared emission, shown in Figures 6 through 9. Decreasing the effective temperature causes the absorption bands of H 2 O and CH 4 to become broader and deeper, impacting all four of the Y , J, H and K bands.Changing the surface gravity impacts the K band, due to increased H 2 − H 2 collisionally induced absorption, which peaks at this wavelength range (see Figure 5).We note that while H 2 − H 2 CIA peaks in the K-band, it can also be an important absorber in the Y , J and H bands.The Y and J band can probe pressures greater than 20 bar.At these pressures H 2 − H 2 CIA is the second most prominent absorber in the Y -band behind the broad red wings of the K I resonance doublet, which is often thought as being the sole opacity source in this wavelength range.In the J-band, H 2 − H 2 CIA is the second most important opacity source behind H 2 O at high pressures.While the H-band does not probe as deep into the atmosphere as the Y or J bands, H 2 − H 2 CIA can be the most dominant absorber in the H-band at high pressures.H 2 − H 2 CIA therefore plays a role in shaping T-dwarf spectra across the near-infrared, not just in the K band.It is particularly important in high-gravity atmospheres in which the peak of the contribution function shifts to higher pressures.
Changes in the metallicity of the atmosphere also impact the H 2 − H 2 collisionally induced absorption.Decreasing the metallicity leads to a relative increase in the abundance of H and He atoms, and a decrease in the abundances of important near-infrared molecular absorbers.This leads to a less opaque atmosphere, shifting the photosphere to higher pressures, and enhancing the H 2 − H 2 collisionally induced absorption.This effect primarily impacts the K band where the H 2 − H 2 collisionally induced absorption peaks, but also the peak fluxes in all the bands.We note that the degeneracy between metallicity and surface gravity for T dwarf atmospheres has been quantified by Zhang et al. (2021b).
Similarly to metallicity, decreasing the oxygen abundance of the atmosphere decreases the abundance of H 2 O, increasing the relative abundance of H atoms and enhancing the H 2 − H 2 collisionally induced absorption in the Y and K band.Additionally, changing the oxygen abundance also manifests as changes in flux in the red side of the H band as seen in Figure 8, due to the increasing CH 4 abundance as the oxygen abundance is decreased (see Figure 4, top panel).Altering the carbon abundance impacts either side of the H band and also the blue side of the K band.The methane abundance increases when increasing the carbon abundance leading to differences in the red side of the H band.The H 2 O abundance also decreases with increasing carbon abundance leading to differences in the blue side of the H and K bands.While these models indicate that non-solar carbon and oxygen abundances could be distinguished from the shapes of the H and K band, we found in our model fitting analysis that there are negligible differences in the spectra of the best-fitting models to GJ 570D when changing the C/O ratio through either the carbon or oxygen abundances (Section 3.2).

GJ 570D
The best-fit model spectrum to the GNIRS spectrum of GJ 570D is shown in the top panel of Figure 10.We find the best-fit parameters to be T eff = 820 K, log(g) = 4.6, Z = 0.0 and C/O = 0.48, using models with varying oxygen abundance.The best-fit parameters are well contained within the model grid parameter space, as shown by the χ 2 surfaces in Figure 11.The minimum χ 2 values are on the order of 10 5 with 3927 degrees of freedom, and almost all the probability contained within the best-fit model.We also compare the best-fit models when fitting grids in which the C/O ratio is changed by altering either the carbon or oxygen elemental abundances, as shown in Figure 12.We find negligible differences in the best-fitting model parameters, except for gravity and C/O ratio where there are small differences which do not impact the spectra meaningfully.Therefore, in the grid-fitting analyses we use only models which alter the oxygen abundance of the model atmosphere.
We also fit the ATMO models to the Y , J, H and K bands of GJ 570D individually, as shown in Figure 13.The individual-band fits reproduce the flux and shape of each band signficantly better than the full-spectrum model fit.However, the best-fit parameters are inconsistent and disagree with the results from the full-spectrum.The best-fit parameters of the full-spectrum and individual band fits are summarised in Table 5.
The age and metallicity of the GJ 570 system are constrained through observations of the host star GJ 570A (see Table 2 of Zhang et al. 2021b, and references therein).High resolution spectroscopy of the host star indicates a metallicity of Z = −0.05− 0.05.The host star C/O ratio has been found to be 0.65-0.97from a stellar abundance analysis (Line et al. 2015).Our best-fit metallicity is consistent with that of the host star.However our best-fit C/O ratio of 0.48 is lower compared to that of the host star (C/O = 0.65 − 0.97).The age of GJ 570A is in the range 1.4 − 5.2 Gyr, determined from stellar isochrone fitting, gyrochronology and stellar activity.Comparing our best-fit T eff , log(g) and R to isochrones from the ATMO 2020 evolutionary models in Figure 14, we can see a disagreement with the values expected for an object with an age range of 1.4 − 5.2 Gyr.Both the gravity and radius are underestimated compared to the evolutionary tracks.
Most previous grid fitting studies of GJ 570D used an IRTF/SpeX prism spectrum, providing 1.0 − 2.5 µm wavelength coverage and R ≈ 150 − 400.Fitting analyses have been done using forward modeling (Oreshenko et al. 2020;Zhang et al. 2021b) and retrieval frameworks (Line et al. 2015(Line et al. , 2017;;Kitzmann et al. 2020;Zalesky et al. 2022;Whiteford et al. 2023), and the results from these studies compared to this work are shown in Table 2.We find a comparable effective temperature and higher gravity compared to Zhang et al. (2021b), and consistent parameters compared to Oreshenko et al. (2020).Notably we find an approximate solar metallicity, whereas the forward modelling work of Zhang et al. (2021b) and the retrieval studies of Line et al. (2015Line et al. ( , 2017)); Kitzmann et al. (2020) all find sub-solar metallicities (see Table 2).Furthermore, our best-fit C/O ratio is approximately solar, whereas the retrieval studies find a super-solar C/O ratio (see Table 2).
These differences could be due to differences in the modelling framework, or due to the higher resolution and signal-to-noise of the GNIRS spectrum.To test this we also fit the ATMO models to the IRTF/SpeX spectrum of GJ 570D (Burgasser et al. 2004).We use the same methods as described in Section 2 to degrade the model spectra.However to implement the wavelength-dependent resolution of the SpeX instrument (R ∼ 150−400, 0.8−2.5 µm, slit width 0.5") we vary the width of the Gaussian kernel used to smooth the spectra as described in Zhang et al. (2021b).The best-fitting model is shown in Figure 15, and we can see that we obtain the same gravity and C/O ratio as the best-fitting model for the higher resolution GNIRS data.We obtain a cooler effective temperature of T eff = 800 K and a sub-solar metallicity of Z = −0.1.The higher resolution GNIRS data appears to better constrain the metallicity of this object, given that the metallicity is constrained to appoximately solar from observations of the host star.We also see that the same discrepancies between model and data exist in the best-fit to the SpeX observations, namely an underprediction of the peak flux in the H band, underprediction of the flux in the red side of the Y band, and an overprediction of flux on the blue side of the J band.

HD 3651B
The best-fit model spectrum to the GNIRS spectrum of HD 3651B is shown in the middle panel of Figure 10.We find the best-fit parameters to be T eff = 740 K, log(g) = 4.3, Z = −0.1 dex and C/O = 0.44, using models with varying oxygen abundance.The best-fit parameters are well contained within the model grid parameter space, as shown by the χ 2 surfaces in Figure 16.The minimum χ 2 values are on the order of 10 5 with 3923 degrees of freedom.
In addition to fitting the full wavelength range, we also fit the Y, J, H and K bands individually as shown in Figure 17.Similar to GJ 570D, fitting the individual bands alleviates the model-data discrepancies seen in the full wavelength range best-fit model, and better reproduces the flux and shape of each band.However, the best-fit parameters are again inconsistent between the individual bands and disagree with the results from the full-spectrum.The best-fit parameters and wavelength ranges of the individual band fits are summarized in Table 6.
The age and metallicity of HD 3651B are constrained through observations of the K0V host star HD 3651A (Liu et al. 2007;Zhang et al. 2021b).High resolution spectroscopy of HD 3651A indicates a metallicity of Z = 0.1−0.2dex, and the host star C/O ratio has been found to be 0.51−0.73from a stellar abundance analysis (Line et al. 2015).Both our best-fit metallicity and C/O ratio are subsolar and below that of the host star, though consistent within our adopted uncertainties of half the model grid spacing.The age of the HD 3651 system is in the range 4.5 − 8.3 Gyr, determined from stellar isochrone fitting, gyrochronology and stellar activity (Zhang et al. 2021b).We compare our best-fit T eff , log(g) and R to that expected for HD 3651B from ATMO 2020 model isochrones.Similar to GJ 570D, the gravity is underpredicted compared to that expected for a 4.5 − 8.3 Gyr object.However unlike GJ 570D, the radius is overestimated.
HD 3651B has been the subject of numerous previous grid fitting studies using its IRTF/SpeX prism spectrum.Burgasser (2007) used a semi-empirical method comparing measured spectral indices to spectral models that were calibrated on GJ 570D to determine the physical properties.Leggett et al. (2007) also obtained a GNIRS spectrum of HD 3651B for comparisons with atmosphere models, and Liu et al. (2007) fit the SpeX spectrum of HD 3651B to several model grids.Most recently, fitting analyses on the SpeX spectrum have been performed using forward modelling (Zhang et al. 2021b) and retrieval frameworks (Line et al. 2015(Line et al. , 2017;;Zalesky et al. 2022).The results from these studies compared to this work are shown in Table 3.We find a lower T eff compared to the Sonora model fits of Zhang et al. (2021b), but a similar surface gravity.The retrieval analyses obtain a larger surface gravity (see Table 4) more in line with that expected from evolutionary models for this object (see Figure 14), and effective temperatures comparable to that found in this work.
We also fit the ATMO models to the IRTF/SpeX spectrum of HD 3651B (Burgasser 2007), and the best-fitting model is shown in Figure 15.We find a larger T eff of 800 K compared to the GNIRS model fit (T eff = 740 K), more in line with the T eff obtained in Zhang et al. (2021b) and older fitting analyses (Burgasser 2007;Leggett et al. 2007;Liu et al. 2011).The gravity, however, is comparable to the GNIRS model fit, and remains low compared to evolutionary models.The metallicity increases by 0.1 dex, but is still low compared to the metallicity of the host star.The C/O ratio does not differ significantly.Interestingly, the radius of the SpeX model fit is lower compared to the GNIRS model fit, and is consistent with that predicted by the evolutionary models for this object (see Figure 14).This occurs due to the higher effective temperature of the SpeX model fit allowing a lower radius to be used when fitting the data.

Ross 458C
The best-fit model spectrum to the GNIRS spectrum of Ross 458C is shown in Figure 10.We find the best-fit parameters to be T eff = 740 K, log(g) = 3.9, Z = 0.1 dex and C/O = 0.48, using models with varying oxygen abundance.The best-fit parameters are well contained within the model grid parameter space.This is shown in the χ 2 surfaces in Figure 18.The minimum χ 2 values are on the order of 10 5 with 3907 degrees of freedom.
Similarly to GJ 570D and HD 3651B, we also fit the ATMO models to the Y , J, H and K bands of Ross 458C individually, as shown in Figure 19.As with the previous two objects, the individual-band fits reproduce the flux and shape of each band significantly better than the full spectrum model fit.The best-fit parameters are again inconsistent and disagree with the results from the full-spectrum.The best-fit parameters of the full spectrum and individual-band fits are summarised in Table 7.
The age and metallicity of the Ross 458 system are constrained through observations of the host binary Ross 458AB.Photometry and medium-resolution (R ≈ 2000) spectroscopy of Ross 458A indicate a super-solar metallicity in the range Z = 0.17 − 0.33 dex (Burgasser et al. 2010;Gaidos & Mann 2014).Our best-fit metallicity of Ross 458C is super-solar, but slightly lower than that observed for the host star at Z = +0.1 dex.The age of Ross 458A is in the range 0.15 − 0.8 Gyr, which has been determined from gyrochronology, stellar activity, and spectroscopic youth indicators (Burgasser et al. 2010;Zhang et al. 2021b).We compare our best-fit T eff , log(g) and R to that expected for Ross 458C from ATMO 2020 model isochrones.Similarly to both GJ 570D and HD 3651B, the gravity of Ross 458C is underpredicted compared to that expected for a 0.15 − 0.8 Gyr object.Additionally the radius of Ross 458C is also lower compared to the evolutionary models.
Ross 458C has been studied extensively in the literature using near-infrared spectra from multiple different spectrographs.Burgasser et al. (2010) and Morley et al. (2012) used a low-resolution (R ∼ 250 − 350) Magellan/FIRE spectrum of Ross 458C and cloudy atmosphere models to derive its physical parameters.Burningham et al. (2011) obtained a Subaru/IRCS spectrum and, combined with photometry and age-estimates of the system, estimated the bolometric luminosity and physical parameters of Ross 458C.Most recently, fitting analyses of the IRTF/SpeX prism spectrum of Ross 458C have been performed using forward modelling (Zhang et al. 2021b) and retrieval frameworks (Zalesky et al. 2022).Gaarn et al. (2023) also performed a retrieval analysis on the Magellan/FIRE spectrum of Ross 458C.The results from these studies compared to this work are shown in Table 7.We find a comparable effective temperature and surface gravity to the retrieval studies, and note that the higher gravity retrieved by Gaarn et al. ( 2023) is more consistent with that expected from evolutionary models.We obtain a comparable gravity (log(g) = 3.9) but lower effective temperature (T eff = 740 K) compared to the Sonora forward modelling work of Zhang et al. (2021b).Furthermore, we get a lower metallicity (Z = 0.1 dex) than the Sonora models obtained on the SpeX spectrum.
We also fit our ATMO models to the IRTF/SpeX spectrum of Ross 458C (Zhang et al. 2021b).The best-fitting model is shown in Figure 15, and we obtain notably different best-fitting parameters compared to our higher-resolution GNIRS data.This is due to a large discrepancy between the observed GNIRS and SpeX spectra of Ross 458C, with the H and K bands in the SpeX spectrum having a larger flux in the normalised spectra.The near-infrared colors of the SpeX spectrum are therefore redder than the GNIRS spectrum for Ross 458C, with J − H = −0.18 and J − K = 0.31 for the SpeX spectrum, and J − H = −0.38 and J − K = 0.19 for the GNIRS spectrum.Near-infrared spectra of Ross 458C have been obtained using IRTF/SpeX in July 2015 (Zhang et al. 2021b), Magellan/FIRE in April 2010 (Burgasser et al. 2010) and Subaru IRCS in May and December 2009 (Burningham et al. 2011).Comparing these spectra to our Gemini/GNIRS spectrum in Figure 20, we can see differences in the J − H and J − K colors.While we see variability in these colors, we do not see any variability in the shape of the individual bands.It should be noted that these spectral differences could be caused by imperfect order stitching in the GNIRS and IRCS data.However the SpeX and FIRE spectra also show slight differences in the J − H and J − K colors, and are not impacted by order stitching since these spectra were obtained in prism mode, which produces continuous wavelength coverage.
Ross 458C has been previously observed to be variable by Manjavacas et al. (2019) who obtained time-resolved spectroscopy using the Hubble Space Telescope Wide Field Camera 3.Over their ∼ 10 h of observing time they found Ross 458C to have peak-to-peak rotational modulations of 2.63 ± 0.02% over the entire 1.1 − 1.64 µm wavelength range, with an estimated 6.75 ± 1.58 h rotation period.This has been attributed to heterogeneous sulfide clouds in the atmosphere.The peak-to-peak variability amplitude is smaller than seen in our comparisons of the spectra of Ross 458C in Figure 20.It should also be noted that Metchev et al. (2015) did not detect this object to be variable in the Spitzer IRAC channels 1 (3.6µm) and 2 (4.5 µm), observing variability amplitudes of < 1.37% and < 0.72% over time baselines of 14 and 7 h respectively.

DISCUSSION
There are three regions of the near-infrared spectrum that are not fitted satisfactorily, and these regions are consistent for the best fits of all three objects.The best-fit models underpredict the flux in the red side (≈ 1.07 − 1.10 µm) of the Y band, overpredict and underpredict the flux on the blue (≈ 1.22 − 1.26 µm) and red (≈ 1.27 − 1.31 µm) side of the J band, respectively, and underpredict the peak flux (≈ 1.54 − 1.60 µm) of the H band.This is seen in the spectral-fitting residuals shown in Figure 21.These discrepancies are alleviated when fitting to the Y , J, H and K bands individually (Figures 13,17 and 19).These individual-band fits reproduce the flux and shape of each band significantly better than the full-spectrum fit, indicating that deficiencies lie in the temperature profile and/or the corresponding chemical abundances.For example, the issues could lie in our assumption of radiative-convective equilibrium, or in using a constant K zz in the convective and radiative regions of the atmosphere (e.g.Mukherjee et al. 2022).Furthermore, potential uncertainties in the order stitching of the GNIRS data by PypeIt could introduce uncertainties to the baseline of the spectral shape.We now discuss each of the Y , J, H and K bands in turn, beginning with the H band.

H band
The peak flux in the H band is underestimated in the best-fit models of GJ 570D, HD 3651B and Ross 458C.The individual spectral band fits better reproduce the flux in the H band, and the reduced χ 2 values are lower in the individual spectral band fits than the full-spectrum fits (see Tables 5, 6 and 7).By comparing the best-fit parameters of the individual H band fits to those of the full-spectrum fits, we can see a clear trend.The best-fit models move to higher gravity and higher metallicity in order to better reproduce the data.To examine why this is the case, we compare the properties at the photosphere (pressure, temperature, abundances) in the H band, for models with log(g) = 4.5, Z = 0.0 and models with log(g) = 5.5, Z = 0.5.
We first begin by examining the contribution function at a given wavelength.The pressure and temperature along with the H 2 O and CH 4 mole fractions for the 1.58 µm contribution function are shown in Figure 22.It should be noted that these are the properties for the contribution function at a single wavelength, and there will be further variation of these properties if a larger wavelength range is selected (for example, the entirety of the H band).In Figure 22 we can see that in highgravity, high-metallicity atmospheres, the contribution function moves to higher pressures and stays at approximately the same temperature.This is the result of two competing effects.Increasing the gravity of a model atmosphere shifts the contribution to higher pressures due to a decrease in the pressure scale height and the optical depth at a given pressure level.At a given pressure level, the temperature decreases with increasing gravity.Conversely, increasing the metallicity of a model atmosphere increases the abundances of metal-bearing species such as H 2 O and CH 4 , increasing the opacity of the atmosphere and shifting the contribution to lower pressures and higher temperatures.In the case shown in Figure 22, the contribution function shifts to higher pressures and remains at approximately the same temperature when increasing gravity and metallicity.Figure 22 also shows that the abundances of important opacity sources H 2 O and CH 4 increase at the location of the contribution function, primarily due to the increase in metallicity.
Figure 23 shows how the H-band model spectra change when increasing the gravity and metallicity of the atmosphere.At log(g) = 4.5 and Z = 0.0 (i.e.parameters representative of the best-fit model to the full spectra of the objects), increasing either the gravity or metallicity does little to increase the peak H-band flux and reproduce the observations.However, it can be seen at high gravities and metallicities, increasing these parameters increases the peak flux in the H band.This shows that the effects gravity/metallicity have on a spectrum is dependent on the metallicity/gravity of the atmosphere.To understand why this is the case, Figure 24 shows the opacities shaping the H band spectrum at the location in the model atmosphere of the peak of the contribution function.The blue side of the H band is shaped by H 2 O opacity, and the red side by CH 4 opacity.The peak flux is controlled by these two opacities sources, along with H 2 − H 2 collisionally induced absorption, the opacity of which depends on gravity and metallicity.Increasing gravity increases H 2 − H 2 CIA due to the higher atmospheric pressures.Increasing metallicity however, decreases the H 2 − H 2 CIA due to more hydrogen being sequestered into molecules such as H 2 O, CH 4 and NH 3 .The abundances of these molecules increase due to the increased availability of metals, increasing the opacity of these species and decreasing H 2 − H 2 CIA.The interplay/competition of the effects of gravity and metallicity on H 2 − H 2 CIA is what drives the individual band fits to converge on high gravity, high metallicity models.These parameters allow the model to control the peak H band flux through the level of H 2 − H 2 CIA.
While these model fits reveal that the peak H-band flux can be affected by the level of H 2 − H 2 CIA, in cooler atmospheres the abundance of NH 3 has also been shown to affect the flux emitted through the H band.The abundance of NH 3 increases in cooler atmospheres, leading to increased opacity in the H band.However, non-equilibrium chemistry models are required to quench the abundance of NH 3 through vertical mixing, increasing the flux in the H-band compared to chemical equilibrium models that is required to reproduce observations of late-T and Y dwarfs (Tremblin et al. 2015).Since we include non-equilibrium chemistry due to vertical mixing in the ATMO models, the continued under-prediction of the flux in the H band could indicate that the abundance of NH 3 needs to be further reduced, implying inadequacies in the chemistry model.To explore this, we removed NH 3 opacity from the calculation of the emission spectra for models close to the best-fitting parameters of GJ 570D and Ross 458C.We note that this treatment is not self-consistent since we do not recompute the T-P profile.The result is shown in Figure 25, in which we can see a H-band brightening in the models without NH 3 opacity which becomes larger when decreasing the effective temperature.However, while removing the NH 3 opacity does brighten the H band, it does not fully account for the under-prediction of flux seen in our best-fits of GJ 570D and Ross 458C.

Y band
The best-fit models for all objects underestimate the flux in the red-side of the Y -band (≈ 1.07 − 1.10 µm).The best-fit model of Ross 458C also overestimates the flux in the blue side of the Y -band.These discrepancies can be seen in the scaled residuals shown in Figure 21.The individual band fits rectify these issues for all objects, and we can see a similar trend as H band, in that the best individual band fit models move to higher gravity and metallicity in order to better reproduce the data (see Tables 5, 6 and 7).
Figure 26 shows how the Y-band model spectra change when increasing the gravity and metallicity of the atmosphere.At log(g) = 4.5 and Z = 0.0, increasing the metallicity does not increase the peak Y-band flux, or alter the shape of the Y-band significantly.However, increasing the metallicity at high gravity (log(g) = 5.5), increases the peak flux of the Y band.Increasing the gravity, regardless of the metallicity, leads to an increase in the flux passing through the Y band, with a larger increase seen at solar metallicity.Figure 27 shows the opacities shaping the Y band spectrum at the location of the peak of the contribution function.The peak flux of the Y band is primarily controlled by the red wing of the potassium resonance doublet, with H 2 O opacity defining the red side of the band.At lower gravities, FeH contributes to the Y band opacity, however increasing gravity lowers its abundance due to the formation of Fe-bearing condensates leading to the rainout of this element.Increasing gravity also leads to an increase in the opacity of potassium and H 2 − H 2 CIA, which becomes the second most important absorber in the Y -band.This is due to the increase in pressure at the peak of the contribution function, suppressing the flux in the Y band.Increasing metallicity, further increases the potassium opacity along with H 2 O opacity, but decreases the H 2 − H 2 CIA due to more hydrogen being contained within metal-bearing molecules.The combination of increasing the potassium opacity by increasing the gravity, and subsequently decreasing the H 2 − H 2 CIA by increasing the metallicity, reshapes the Y band such that the individual band fits converge on a higher-gravity, higher-metallicity model.

J band
The full-spectrum fits of all three benchmark objects do not correctly reproduce the shape of the Jband.They suffer from the same deficiencies; an overprediction of flux on the blue (≈ 1.22−1.26µm) side of the band and an underprediction of flux on the red (≈ 1.27 − 1.31 µm) side of the band.This is seen in the scaled spectral-fitting residuals in Figure 21.When performing individual spectral band fits, the shape of the J-band is better reproduced by the best-fit models, as shown in Figures 13, 17 and 19.In these individual band fits, the best-fit models have different parameters for all three objects (Tables 5, 6 and 7).The best-fit model for GJ 570D shifts to lower T eff , higher gravity, and a super-solar metallicity and C/O ratio.The best-fit parameters for HD 3651B shift to higher gravity (edge of the model grid), a more super-solar metallicity, and C/O ratio while remaining at the same T eff as the full spectrum fit.Finally, the best-fit parameters for Ross 458C follow a different trend, shifting to a lower gravity, sub-solar metallicity and C/O ratio while remaining at the same effective temperature.We can therefore see that the individual J-band fits show a similar trend to the H and Y bands, in that the best-fit parameters shift to high gravity and metallicity (except for Ross 458C).
Figure 28 shows how the J-band model spectra change when increasing gravity and metallicity in the atmosphere.At log(g) = 4.5 and Z = 0.0, increasing the metallicity does not impact the J band.However increasing the metallicity at high gravity (log(g) = 5.5) increases the peak J-band flux and alters the shape of the band.Increasing gravity at solar and super-solar metallicities does not significantly impact the J-band, apart from the depth of the K I resonance doublet at ≈ 1.25 µm which is strongest at log(g) = 4.5, Z = +0.5. Figure 29 shows the opacities shaping the J-band spectrum at the location of the peak of the contribution function.The J band is primarily shaped by H 2 O absorption, with potassium resonance doublets (≈ 1.175 µm and ≈ 1.250 µm) and H 2 − H 2 CIA contributing to the opacity.At high gravities, the strength of H 2 − H 2 CIA increases as the contribution function shifts to higher pressures, and the K I resonance doublets decrease in strength as this element begins to rainout of the atmosphere due to condensation.Increasing metallicity strengthens the K resonance doublets, and decreases H 2 − H 2 CIA.The individual J-band fits of GJ 570D and HD 3651B therefore shift to higher gravity, higher metallicity atmospheres to better reproduce the shape of the J band through these opacity sources.We also note that decreasing the effective temperature (as in the case of the GJ 570D best-fit parameters) also alters the shape of the J-band due to an increase in strength of the water absorption bands, (Figure 7).

K band
The shape and flux level in the K band is well reproduced by the best-fit models for all objects.The peak flux of the K-band is slightly underestimated in the best-fit model of HD 3651B, and the best-fit model of Ross 458C slightly underestimates the flux on the blue side of the K-band.The individual spectral band fits visually do a better job at reproducing the depth of the water absorption features in the blue side of the K band for Ross 458C and the peak flux for HD 3651B.

SUMMARY AND CONCLUSIONS
In this paper we have expanded the ATMO 2020 model atmospheres (Phillips et al. 2020) to include non-solar abundances in the mid/late-T dwarf regime (T eff = 700 − 900 K).This is done through altering the bulk metallicity of the atmosphere and the C/O ratio by separately changing the carbon and oxygen elemental abundances.We use these new models to assess the spectral signatures of non-solar abundances and C/O ratios at the resolution of the Gemini/GNIRS instrument in the 0.8 − 2.5 µm wavelength range (Figures 6-9).We find that altering the oxygen abundance has effects in the Y and K bands that are degenerate with metallicity, but can be distinguished by changes in flux on the red side of the H band. Altering the carbon abundance impacts the shape of the H band and the flux on the blue side of the K band.Our results indicate that non-solar abundances of oxygen and carbon could be distinguished from one another at this resolution based on the shapes of the H and K bands.
We compare our models with new Gemini/GNIRS R ≈ 1700 spectra of three benchmark T dwarfs, GJ 570D, HD 3651B and Ross 458C.Our model fitting analysis finds negligible differences in the best-fitting model parameters when changing the C/O ratio through either the carbon or oxygen abundances (Figure 12), and we therefore use only models which alter the oxygen abundance of the model atmosphere.We find approximately solar C/O ratios for all objects and best-fitting parameters (T eff , log(g), Z) broadly consistent with other analyses in the literature (Tables 2, 3 and 4, Figure 10).
To examine the benefits that R ≈ 1700 spectra provide over the more widely available low-resolution (R ∼ 150) IRTF/SpeX spectra, we also fit the ATMO models to the SpeX spectra of GJ 570D, HD 3651B and Ross 458C (Figure 15).For GJ 570D, we obtain the same gravity and C/O ratio as the best-fitting model to the higher-resolution GNIRS data.However, the best-fit model to the SpeX data gives a ≈ 50 K cooler effective temperature and a somewhat sub-solar metallicity (−0.1 dex).The higher-resolution GNIRS data appears to better establish the metallicity of this object, given that the metallicity is constrained to be approximately solar from observations of the host star.Fitting the SpeX spectrum of HD 3651B we obtain the same gravity but a larger effective temperature compared to the best-fit of the GNIRS data, more similar to previous grid fitting analyses.The best-fit metallicity increases by 0.1 dex, but is still low compared to the metallicity of the host star.Finally, we also fit the ATMO models to the SpeX spectrum of Ross 458C.However this analysis is complicated by the fact that the GNIRS and SpeX spectra differ substantially, particularly in the J − H and J − K colors.Further investigation reveals that the near-infrared (0.8 − 2.5 µm) spectra of Ross 458C differs across different epochs taken with multiple spectrographs, indicating that this object could be potentially variable (Figure 20).
There are three regions of the near-infrared spectrum that are not fitted satisfactorily, and these regions are consistent in the best fits of all objects.The best-fit models under-predict the flux on the red side (≈ 1.07 − 1.10 µm) of the Y band, overpredict and under predict the flux on the blue (≈ 1.22 − 1.26 µm) and red (≈ 1.27 − 1.31 µm) side of the J band respectively, and under-predict the peak flux (≈ 1.54 − 1.60 µm) of the H band.These discrepancies can be seen in the spectral residuals shown in Figure 21.To investigate these discrepancies, we also fit the ATMO models to the Y , J, H and K bands individually (Figures 13,17 and 19).These individual-band fits reproduce the flux and shape of each band significantly better than the full-spectrum fit.However, the resulting best-fit parameters are not consistent with each other and disagree with the results from the full spectrum.This leads us to conclude that the opacities used in the atmosphere model are sufficiently complete and accurate at this spectral resolution, and deficiencies are instead more likely to lie in the temperature profile and/or the corresponding chemical abundances.
While the individual band best-fit parameters are inconsistent with the results from the full spectrum, we find a consistent shift to high-gravity, high-metallicity atmospheres.This is particularly apparent in the individual H-band fits, which show a clear trend by moving to high gravity and high metallicity to better reproduce the H-band spectra of all objects.By examining the model atmosphere properties at the H-band photosphere, we find that this is due to the interplay/competition of the effects of gravity and metallicity on H 2 − H 2 CIA (Figures 23 and 24).These parameters allow the model to control the peak H-band flux through the level of H 2 − H 2 CIA and drive the individual-band fits to converge on high-gravity, high-metallicity models.This same process occurs in the individual-band fits of the Y and J bands, for which the best-fit parameters again shift to high gravity and high metallicity (except for the J-band fit of Ross 458C), with H 2 − H 2 CIA playing a key role (see Figures 24,27 and 29).Such high-gravity, high-metallicity atmospheres are unrealistic for all objects given the age and metallicity determined from their host stars, and thus it is clear that these discrepancies are due to incorrect or missing physics in the atmosphere models.
To explore the underprediction of the H-band further, we removed the NH 3 opacity from the calculation of the model emission spectra, to simulate a more extreme quenching of NH 3 due to vertical mixing than currently modeled by our non-equilibrium chemistry scheme.The quenching of NH 3 has been shown to be particularly important in reproducing the H-band flux in observations of late-T and Y dwarfs (Tremblin et al. 2015).We find that while removing NH 3 opacity does brighten the H band, it does not fully account for the under-prediction of flux seen in our best-fits of GJ 570D and Ross 458C (Figure 25).A full comparison of the impact the choice of non-equilibrium chemistry scheme (e.g.relaxation scheme [Tsai et al. 2018], kinetics network [Venot et al. 2012], reduced kinetics network [Venot et al. 2019]) has on the predicted H-band peak flux in cool brown dwarfs is needed to better understand the discrepancies seen between model and data.
The wide-seperation companions studied in this work provide crucial benchmarks thanks to the age and composition constraints derived from their host stars.In future works, we will compare our new non-solar ATMO models to a wider sample of high-quality, medium-resolution Gemini/GNIRS spectra of cool T type brown dwarfs.This wider sample will include free-floating brown dwarfs in the solar neighbourhood and confirmed members of nearby young (< 100 Myr) moving groups (Zhang et al. 2021a).This will produce measurements of the fundamental properties and C/O ratios for cool brown dwarfs, and allow us to compare the different formation pathways among companions, young moving group members and field brown dwarfs.This research was funded in part by the Gordon and Betty Moore Foundation through grant GBMF8550 and by the National Science Foundation through grant NSF AST 1838016.Z. Z. acknowledges support from the NASA Hubble Fellowship grant HST-HF2-51522.001-A awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under contract NAS5-26555.This work has benefitted from The UltracoolSheet at http://bit.ly/UltracoolSheet, maintained by Will Best, Trent Dupuy, Michael Liu, Aniket Sanghi, Rob Siverd, and Zhoujian Zhang, and developed from compilations by Dupuy & Liu (2012); Dupuy & Kraus (2013); Deacon et al. (2014); Liu et al. (2016); Best et al. (2018Best et al. ( , 2021)); Schneider et al. (2023); Sanghi et al. (2023).This work was enabled by observations made from the Gemini North telescope, located within the Maunakea Science Reserve and adjacent to the summit of Maunakea.We are grateful for the privilege of observing the Universe from a place that is unique in both its astronomical quality and its cultural significance.
Software: Astropy (Astropy Collaboration et al. 2013Collaboration et al. , 2018Collaboration et al. , 2022)), IPython (Pérez & Granger 2007), Numpy (Harris et al. 2020), Scipy (Virtanen et al. 2020), Matplotlib (Hunter 2007).b Luminosity analysis refers to using spectroscopic and photometric observations, along with atmosphere models, to estimate the bolometric luminosity of an object.Combined with age estimates, the bolometric luminosity can be used to estimate a range of T eff and log(g) from evolutionary tracks, as described in Saumon et al. (2006).
Table 3. Parameters derived for HD 3651B from spectroscopic analysis

b
A semi-empirical analysis comparing measured spectral indices to those from model spectra that have been calibrated on GJ 570D, as described in Burgasser et al. (2006) and Burgasser (2007).
Tucson refers to models from the Tucson group (Burrows et al. 2003;Burrows et al. 2006).
d SM refers to atmosphere models from the Saumon & Marley group (Saumon et al. 2006;Saumon & Marley 2008;Saumon et al. 2012;Morley et al. 2012).b Luminosity analysis refers to using spectroscopic and photometric observations, along with atmosphere models, to estimate the bolometric luminosity of an object.Combined with age estimates, the bolometric luminosity can be used to estimate a range of T eff and log(g) from evolutionary tracks, as described in Burningham et al. (2011).
c SM refers to atmosphere models from the Saumon & Marley group (Saumon et al. 2006;Saumon & Marley 2008;Saumon et al. 2012;Morley et al. 2012).Y band  J band         Wavelength ( m) Figure 12.The best-fit ATMO spectra to the GNIRS spectrum of GJ 570D when varying carbon (green) and oxygen (orange) in the model atmosphere.The uncertainties on the GNIRS spectrum are plotted as the blue line.Note that these model fits are limited to log(g) ≥ 4.5, hence why the best-fit parameters when varying oxygen differ to those presented elsewhere in the paper (e.g. Figure 10).

a
These values are fixed in the model fitting analysis.

a
These values are fixed in the model fitting analysis.

a
These models assume solar abundances.
Table 6.Best-fit parameters derived for HD 3651B full-spectrum and individual-band fits Fitted wavelength range (µm) χ 2 Table 7. Best-fit parameters derived for Ross 458C full-spectrum and individual-band fits Fitted wavelength range (µm)

Figure 5 .
Figure 5. Top: Contour plot of the normalised contribution function, illustrating the pressure levels of the model atmosphere contributing to the emission from the top of the atmosphere at a given wavelength.The model has T eff = 800 K, log(g) = 4.5, Z = 0.0 and C/O = 0.55.Middle & Bottom: Abundance weighted absorption cross-sections of the main atomic and molecular opacity sources in a typical T dwarf atmosphere, with the same parameters as the top panel.The two panels show the cross-sections at two different atmospheric levels.

Figure 6 .
Figure 6.Model ATMO Y -band emission spectra varying parameters around a T eff = 800 K, log(g) = 5.0, Z = 0.0 dex and C/O = 0.55 (solar) model atmosphere.The oxygen and carbon abundances are changed through the scaling factors X O and X C respectively.Changes in the flux when increasing and decreasing the model grid parameters by the values indicated on the plot, are plotted as blue and yellow shaded regions respectively.The spectra have been degraded to the GNIRS resolution R = 1700, and the models have been offset for clarity.

Figure 7 .
Figure 7. Same as Figure 6 but for the J band.

Figure 8 .
Figure 8. Same as Figure 6 but for the H band.

Figure 9 .
Figure 9. Same as Figure 6 but for the K band.

Figure 10 .Figure 11 .
Figure10.The best-fit ATMO spectra (orange) to the GNIRS spectra (black ) of GJ 570D (top), HD 3651B (middle) and Ross 458C (bottom) .The uncertainties on the GNIRS spectra are plotted as the blue line at the bottom of each plot.

Figure 17 .
Figure17.Y -, J-, H-, and K-band (from top to bottom) spectra of HD 3651B including the best-fitting model to the full wavelength range (orange, same as Figure10) and the best-fitting model to the individual wavelength ranges shown in each panel.

Figure 18 .
Figure 18.Ross 458C χ 2 ν surfaces for each parameter in the model grid.There are ν = 3907 degrees of freedom.The location of the best-fit value (χ 2 ν = 11.7) is indicated by a red cross in each panel.
Figure19.Y -, J-, H-, and K-band (from top to bottom) spectra of Ross 458C including the best-fitting model to the full wavelength range (orange, same model as Figure10) and the best-fitting model to the individual wavelength ranges shown in each panel.

Figure 21 .
Figure21.Top: Spectral-fitting residuals for the 3 benchmark brown dwarfs analysed in this work.Each residual has been normalised by the object's median flux at the peak of the J-band.The residuals have been smoothed using a median filter with a window size of 10 pixels.Bottom: The GNIRS spectrum of GJ 570D as a reference.

Figure 23 .Figure 24 .
Figure23.Model H-band emission spectra for a T eff = 800 K atmosphere, and varying gravity and metallicity as indicated in the legend.Changing gravity/metallicity has a larger effect on the peak H-band flux at high metallicity/gravity.

Figure 27 .Figure 29 .
Figure 27.Same as Figure 24 but showing absorption cross-sections at the location of the peak of the contribution function at 1.09 µm.
Y S/N J S/N H S/N

Table 2 .
Parameters derived for GJ 570D from spectroscopic analysis

Table 4 .
Parameters derived for Ross 458C from spectroscopic analysis T eff (K)

Table 5 .
Best-fit parameters derived for GJ 570D full-spectrum and individualband fits Mole fractions of the dominant infrared opacity sources in cool brown dwarf atmospheres as a function of C/O ratio for two different temperatures.The two panels approximately represent the temperatures at which the Y and J band (∼ 2000 K) and H and K band (∼ 1000 K) flux peaks originate.Solid and dashed lines correspond to varying the oxygen and carbon abundances respectively.
(Phillips et al. 2020) our best-fit gravity, effective temperature and radius to the ATMO 2020 evolutionary tracks(Phillips et al. 2020).The best-fit values to the GNIRS spectrum of each object are shown as the colored stars with error bars.The best-fit values from spectroscopic analyses in the literature listed in Tables 2, 3 and 4 are shown as smaller color circles with error bars.Note that some literature sources do not provide best-fit radii.The corresponding colored shaded regions are the range of estimated ages of each benchmark system, with the dotted lines as model isochrones.
(Zhang et al. 2021b04and K-band (from top to bottom) spectra of GJ 570D including the best-fitting model to the full wavelength range (orange, same as Figure 10 model) and the best-fitting model to the individual wavelength ranges shown in each panel.Figure15.The best-fit ATMO spectrum (orange) to the IRTF/SpeX spectra (black ) of GJ 570D (top)(Burgasser et al. 2004), HD 3651B (middle)(Burgasser 2007), and Ross 458C (bottom)(Zhang et al. 2021b).The uncertainties on the SpeX spectra are plotted as the blue line in all panels.