A Uniform Retrieval Analysis of Ultra-cool Dwarfs. IV. A Statistical Census from 50 Late-T Dwarfs

The spectra of brown dwarfs are key to exploring the chemistry and physics that take place in their atmospheres. Late-T dwarf spectra are particularly diagnostic due to their relatively cloud-free atmospheres and deep molecular bands. With the use of powerful atmospheric retrieval tools applied to the spectra of these objects, direct constraints on molecular/atomic abundances, gravity, and vertical thermal profiles can be obtained enabling a broad exploration of the chemical/physical mechanisms operating in their atmospheres. We present a uniform retrieval analysis on low-resolution IRTF SpeX near-IR spectra of a sample of 50 T dwarfs, including new observations as part of a recent volume-limited survey. This analysis more than quadruples the sample of T dwarfs with retrieved temperature profiles and abundances (H$_2$O, CH$_4$, NH$_3$, K and subsequent C/O and metallicities). We are generally able to constrain effective temperatures to within 50K, volume mixing ratios for major species to within 0.25dex, atmospheric metallicities [M/H] to within 0.2, and C/O ratios to within 0.2. We compare our retrieved constraints on the thermal structure, chemistry, and gravities of these objects with predictions from self-consistent radiative-convective equilibrium models and find, in general though with substantial scatter, consistency with solar composition chemistry and thermal profiles of the neighboring stellar FGK population. Objects with notable discrepancies between the two modeling techniques and potential mechanisms for their differences, be they related to modeling approach or physically motivated, are discussed more thoroughly in the text.


Introduction
Brown dwarfs are substellar objects whose masses are intermediate between the latest M-type stars and the most massive planets (Hayashi & Nakano 1963;Shu 1977;Becklin & Zuckerman 1988;Rebolo et al. 1995;Oppenheimer et al. 1995;Saumon et al. 2006). Similar to stars, brown dwarfs form from interstellar molecular gas cloud core collapse (Uehara & Inutsuka 2000;Bate et al. 2002;Krumholz et al. 2005;Whitworth & Stamatellos 2006;Chabrier et al. 2007; Corresponding author: Joe Zalesky jazalesk@asu.edu Whitworth et al. 2007;Hennebelle 2012), but do not achieve masses high enough to sustain core H-fusion over their lifetime (Burrows et al. 2001). As the effective temperatures of brown dwarfs are much cooler than those of stars (< 2500K), molecules and condensates form in their photospheres and dominate the spectral energy distribution. It is through analyses of these spectra that we are able to infer the nature of brown dwarf atmospheres and how they evolve overtime. Hence we are able to use them as unirradiated laboratories, which provides a valuable baseline for understanding extrasolar planet atmospheres.
Fundamental properties of brown dwarfs include their mass, radius, gravity, effective temperatures (T eff ), and el-emental abundances (see review in Kirkpatrick 2005;Marley & Robinson 2015). Robust masses can be derived through dynamical means (Dupuy & Liu 2017) and radius is inferred via distance and bolometric luminosity. Determining gravity, effective temperatures, and abundances for brown dwarfs can be more challenging than it is for stars due to the lack of clear atomic lines for which classic single or multi-line spectral analyses can be performed. As such, atmospheric models which properly incorporate molecules and condensate species play a more important role in these determinations.
A common modeling approach for determining these properties is through comparisons of observations to grids of model spectra computed with self-consistent onedimensional radiative-convective equilibrium models (e.g., Allard et al. 1996;Marley et al. 1996;Burrows et al. 2001) or using more recent grid-fitting methods (Zhang et al. 2021a,b). This approach typically relies upon the use of a priori chemical/physical assumptions such as thermochemical equilibrium, molecular/atomic abundances, assumed atmospheric chemistry paradigms, and 1D radiative-convective equilibrium (Burrows et al. 2001;Marley & Robinson 2015). These assumptions are made to reduce the dimensionality of the inference problem to just a handful of parameters such as T eff , log(g), and a composition parameter like metallicity [M/H] (in some cases, alpha-element enhancement e.g. Husser et al. 2013).
A few key issues arise with this method. First, the low model dimensionality restricts any inference solely to the dimensions specified for the pre-computed grid. Second, the choice of inference tool is often not rigorous and typically does not account for grid-interpolation uncertainties (e.g., often a simple chi-square type minimizer is combined with a multi-linear-type interpolater) and can result in artificially precise constraints. Thirdly, often times the overly restrictive assumptions lead to poor model spectra fits to the data (e.g., Patience et al. (2012)), leading one to question the validity of the self-consistent modeling assumptions. Zhang et al. (2021a,b) sought to remedy the second and third issues through the use of a modernized self-consistent grid (Marley et al. (2021)) combined with the Starfish tool (Czekala et al. (2015)), which attempts to account for the finite model grid spacing, interpolation uncertainties, and datamodel misfits. However, such an approach is still restricted to a priori physical assumptions within the grid itself.
A more flexible alternative known as atmospheric retrieval, has shown success in providing robust model-data fits to low resolution spectral observations of T and L dwarfs (Line et al. 2014(Line et al. , 2015Burningham et al. 2017;Gonzales et al. 2020). Originally developed to determine temperatures and abundances from spectral soundings of the Solar System planets (Fletcher et al. 2007;Irwin et al. 2008;Greathouse et al. 2011), this technique relies on the use of a forward radiative transfer model that relaxes many of the self-consistent grid model assumptions at the expense of requiring many free parameters. The fundamental philosophy of the retrieval approach, in contrast to the grid approach, is that much of the physical/chemical mechanisms operating in low temperature atmospheres are not understood well enough to build accurate, fully self-consistent models. The aim of atmospheric retrieval is to directly determine, from the spectra the vertical temperature profiles and atmospheric composition. This approach has recently become prolific in the extrasolar planet atmosphere studies (e.g., see review by Madhusudhan (2019)). Line et al. (2015) applied the atmospheric retrieval approach to low resolution SpeX data of the benchmark late-T dwarfs, Gl 570D and HD 3651B. Late-T dwarfs were specifically chosen to mitigate the impact of uncertain cloud properties and the presentation of strong molecular absorption features from water and methane. Using a radiative transfer forward model with ∼30 free parameters combined with Markov Chain Monte Carlo (MCMC) inference (Foreman-Mackey et al. 2013), they were able to obtain bound constraints on the molecular mixing ratios for [H 2 O, CH 4 , NH 3 and Na+K], the vertical thermal profiles (temperature vs. pressure), gravity, and photometric radii (given the parallactic distances). The key findings were (1) that ammonia abundance could be constrained from low-resolution nearinfrared spectra alone, a surprise given the lack of obvious visual spectral features (typical of longer wavelengths Saumon et al. 2006 or higher resolutions (Canty et al. 2015)); (2) the retrieved molecular (and alkali) abundances were consistent with self-consistent chemical predictions (albeit constant-with-altitude mixing ratios were assumed in the retrieval); (3) derived metallicities and carbon-to-oxygen ratios were consistent with their host star abundances; and (4) the vertical thermal temperature profiles agreed with radiativeconvective equilibrium. Taken together, these findings lend support that the retrieval paradigm can be used as a complimentary tool to grid-modeling for inferring fundamental brown dwarf atmospheric properties.
Having validated the retrieval methodology against late-T benchmark systems, Line et al. (2017) performed a systematic retrieval analysis on the spectra of 11 late-T dwarfs (T7-T8, spanning 600 -800K) available in the SpeX prism library (Burgasser 2014). This uniform analysis found that (1) the large number of free parameters required in retrievals, compared to self-consistent grid models (27 vs. 4 parameters) is justified owing to their much better fits; (2) the T7/T8 atmospheres are cloud free (upper limits on the cloud optical depth of unity were obtained); (3) the temperature profiles for all objects were again consistent with radiative-convective equilibrium; (4) the retrieved gravities, radii, and inferred effective temperatures agreed with evolution model predictions; (5) abundances for ammonia, methane, and water were found to be constant with effective temperature but a strong decreasing trend in the alkali abundances was observed to occur with decreasing effective temperature; and (6) the late-T dwarf ensemble had somewhat lower metallicities and higher carbon-to-oxygen ratios than the local FGK stellar population. These findings provided a first look at the carbon-tooxygen ratio and metallicities of a sample brown dwarfs as well as the first direct determination of the possible influence of alkali rainout on their abundances.
Building upon Line et al. (2015Line et al. ( , 2017, Zalesky et al. (2019) extended the uniform retrieval analysis into the cooler Y-dwarfs. This sampled comprised eleven Y dwarfs and three T dwarfs as observed with the Hubble Space Telescope's Wide Field Camera 3 (WFC3) (Schneider et al. 2015). The conclusions were similar to the late-T results Line et al. (2017), finding that (1) the retrieved temperature profiles for most objects were consistent with radiativeconvective equilibrium predictions; (2) water and methane abundances were consistent with equilibrium chemistry; (3) the ammonia abundances showed an upward rise with decreasing temperature, with scatter consistent with both equilibrium and quenched abundances; (4) constraints and/or upper limits on the alkali abundances consistent with rain-out predictions; (5) findings of very high gravities, pushing the limits of the evolution models; and (6) elemental abundance ratios broadly inline with those from the local FGK star population. They also compared the results of a self-consistent grid fit to the retrieval results finding that (1) the retrievals fit better and the large number of free parameters were justified, and (2) constraints on common parameters (T eff , log(g), metallicity, and radius) were inconsistent at 1σ.
Additional works by Burningham et al. (2017) and Gonzales et al. (2020) focused on determining the fundamental properties of L-dwarfs using similar retrieval methods. Ldwarfs are more complicated to characterize due to the presence of condensate clouds and additional higher temperature species (hydrides, oxides), and the reduced vertical grasp on the thermal structures. Overall these cloudy investigations showed that the retrieved temperature structures could be degenerate with the presence of clouds, but that plausible abundances of the hydrides/oxides could be retrieved, opening up the possibility of abundance determinations at higher temperatures.
In this work, we extend the late-T analysis in Line et al. (2017) to a broader sample of 50 T7-T9-dwarfs which were contained in a volume limited survey (Best et al. (2020)), including new IRTF SpeX spectroscopy (Zhang et al. 2021a,b). Our work here represents the largest uniform retrieval analysis on a nearly complete sample of late-T dwarfs, with a factor of 5 larger sample over the analysis in Line et al. (2017). As in the past papers of this series, we focus on late-T dwarfs as these have been shown to be largely free of influence from clouds and they present deep methane and water absorption features, enabling simultaneous carbon and oxygen constraints.
We follow closely the methodologies and analysis as in Line et al. (2017) and Zalesky et al. (2019). In Section 3, we briefly discuss the source of our data and review our retrieval methodology. Section 4 presents our retrieved constraints for our thermal structures (4.1), composition (4.2), and evolutionary parameters (4.3). Finally we briefly summarize our findings in Section 5.

Data
Our entire dataset is derived from late-T dwarf IRTF SpeX prism observations (0.95 -2.5 µm, R∼120) largely described in (Zhang et al. 2021a,b, see Table 3). There are a total of 55 objects with SpeX spectra with 39 available within the SpeX Prism Library, and 16 from our own observational campaigns. Briefly, the sample is comprised of objects from the volume-limited (d <25 pc) survey in Best et al. (2020) with our SpeX objects out to 20 pc. 54 of the 55 have well measured parallaxes and all spectra are flux calibrated using MKO H-band photometry. Our analysis, for uniformity purposes, also includes the 11 objects from Line et al. (2017). Due to low signal-to-noise, our final analysis excludes 5 objects, for a total of 50 in our analysis. This subset is high- T9   T8   T7   T6   T5   T4  T3 T2  T1 T0  L9  L8 L7 L6 L5 L4 Figure 1. A color-magnitude diagram of a population of field brown dwarfs (Best et al. 2018(Best et al. , 2020, and our sample of brown dwarfs from this study (T5-T9). This collective ensemble shows L and T dwarfs, in descending order from the top of the diagram to the bottom. References-References for both spectral types and parallax distances in this table are (1)  lighted in a color-magnitude diagram shown in Figure 1 with photometry and parallax measurements enumerated in Table  1. As in our past works we only analyze every 3rd spectral point as SpeX samples the line-spread shape with ∼2.5 pixels.

Methods
Our retrieval framework follows closely that described in Line et al. (2017) and Zalesky et al. (2019). The model computes the thermal emission spectrum of a brown dwarf and includes 31 free-parameters to describe the atmosphere: constant-with-altitude volume mixing ratios (VMR's) for: H 2 O, CH 4 , NH 3 , CO, CO 2 , H 2 S, Na, and K (8 in total), gravity, radius-to-distance scale factor (R/D) 2 , 15 independent temperature-pressure (TP) profile points implemented within a Gaussian-Process-like smoothing framework, and a simple cloud parameterization (Burrows et al. 2006); summarized in Table 2.
Compared to our previous work we have made several adjustments to this model. We have separated the alkali constraints into distinct Na and K free parameters (as opposed to Na+K fixed at a solar ratio), removed H 2 S due to low predicted abundances at these temperatures, and updated the alkali metal wing profiles (Allard et al. 2016). All other opacity sources are identical to Line et al. (2017) and Zalesky et al. (2019). Model cross-section sampling is used at a constant R=10000, a resolution proven sufficient for interpreting SpeX brown dwarf observations (Line et al. 2015). Parameter prior ranges are the same as those of Line et al. (2015) which are restated in Table 3. Fits were performed as in the mentioned previous studies using the emcee package (Foreman-Mackey et al. 2013) with 224 walkers run out to 60,000 iterations. To assure convergence several test cases were ran to between 120K and 1 million iterations with no significant differences in model posteriors. emcee, like all MCMC-based methods, requires an initial guess to initiate the walkers. The initial guess is based on a "Gaussian-ball" about a lose by-eye fit to the spectra of each object. The final solutions were found insensitive to the initial guess as was the case in previous publications.
To reduce computational runtime for such a large number of objects, our radiative transfer core (Appendix A in Lacis & Oinas (1991)) was rewritten to make use of graphics processing units (GPUs). This was done using the Anaconda Numba guvectorize framework on NVIDIA Tesla V100 GPUs. Forward model times improved by a factor of ∼100 (0.01s or so per model, at an R=10,000 over the 0.95 -2.5µm wavelength range). Given the limited memory of the GPUs (32 GB), we could only run up to 16 simultaneous CPU threads at a time. The overall computational improvement between this work and that of Line et al. (2017) is about a factor of 10 (about 6 hours to hit 60,000 iterations). GPU and CPU routines were tested to produce identical model spectra with no impact on any science results.

Results & Discussion
As it is largely uninformative to show the full posteriors and discuss each individual object (available in the attached Zenodo link 1 ), here we simply summarize the key results broken down by thermal structure, composition (molecular abundance trends/chemistry and elemental abundances), and evolutionary parameters. Figures 2 and 3 provide a snapshot of the fits and temperature-pressure profiles, respectively, of a sub-sample of representative objects. The best fits and temperature profiles for all objects in our sample are provided in the Appendix ( Figures A2 and A3). As has been discussed in previous papers in this series, the visual quality of the model fits in Figures 2 and A2 are substantially better than those by more traditional grid-based models, with residuals following only the uncertainty in the data itself with little other deviation. errorbar inflation exponent γ, β TP-profile smoothing hyperparameters (eq. 5, Line et al. (2015)) κP 0 , P0, α Cloud opacity profile parameters (eq. 1, , cloud base opacity, cloud base pressure, cloud fractional scale height) Inverse Gamma(Γ(γ; α, β)), α=1, β= 5 ×10 −5 log(κ C ), log(P 0 ), α (-12, 0), (2.3, 2.8), [0,10) Table 4 summarizes the nominal constraints for the key properties of individual objects. As in our past works, the effective temperatures are derived by integrating over an ensemble of best fits for each object, extrapolated out to between 0.7 -20 µm, and radius is derived from the retrieved (R/D) 2 scale factor and the measured distances (from Table 1). The elemental C and O abundances are derived from the retrieved molecular abundances (more details below). When comparing the retrieved quantities to those predicted from self-consistent grids (given an effective temperature and gravity and assuming solar abundance chemistry and cloud free), we refer to those produced by the ScCHIMERA model

Thermal Structures
Late-T type brown dwarfs are ideal for temperature profile retrievals due to the lack of obscuring, optically thick condensate clouds that would otherwise limit the range of pressures that could be probed via spectroscopy. The high degree of flexibility of the thermal profile parameterization allows us to compare thermal profiles derived directly from observational spectra against assumptions often made in more selfconsistent models. Since the thermal structure is a reflection of the total energy balance, identifying similarities and differences between retrieved profiles and those predicted from 1D radiative-convective-thermochemical equilibrium (RCE) for such a large number of objects enables us to assess and quantify potential sources of heating outside of the standard RCE assumptions. Figure 3 summarizes (see Appendix, Figure A3 for all objects in the sample) the retrieved thermal structures for nine representative objects (red) compared with 1D RCE predicted TP profiles (blue). Uncertainties for the RCE temperature profiles are derived by propagating uncertainties in the retrieved log(g), T eff , metallicity and C/O. Overlaid for context are condensation curves for several notable condensate species (Morley et al. 2012(Morley et al. , 2014. As in previous studies we reiterate that our SpeX observations largely probe only the region between 0.1 and 10bars with the rest of the profile's constraints being largely dependent upon the smoothing prior; any structure outside this range should be interpreted with caution. We find that the uncertainties on the retrieved profiles are similar to previous studies, with 1σ errors on the order of 100K with larger deviations being owed to larger uncertainties on the spectra themselves (Line et al. 2015Zalesky et al. 2019). While some objects appear to be in agree-ment with the assumption of 1D RCE to 2σ (see WISEPC J0223 in Figure 3, middle left) 18 objects are discrepant with the 1D RCE predicted profiles. There is a systematic offset whereby the retrieved profiles for a given T eff and log(g) are hotter (at a given pressure) than the corresponding 1D RCE predicted profiles. In addition to this systematic offset, 9 objects display an interesting "kink" in the profiles near the 10bar pressure level.
We first tested if the overall offset was physical or a simple result of over or under-constraining either model. To make clear, the error bounds for the grid model are calculated by taking the constraints from the retrieval model and propagating them into error bounds for the grid's thermal profile using several key parameters: T eff , log(g), metallicity, C/O, as well as a reasonable range of K zz values for late-T dwarfs. We tested the robustness of these constraints by artificially inflating the errors of the retrieval results by 1σ to see if the resulting thermal profiles would then overlap. Inflating the errors in this way removed the evidence of any kink structure near 10bars. As this pressure range is near the edge of the contribution functions that our spectra probe, the presence of such a feature should be interpreted with caution.
The most likely reason for the retrieval profiles being systematically hotter than the 1D RCE profiles is owed to how T eff is computed in the retrieval model coupled with the gridbased models being forced onto RCE. As explained in previous papers, a sample of several thousand spectra are taken from the posterior and propagated out to roughly 20µm. This flux is then integrated and T eff is computed from L Bol . However, its important to note that a majority of the output flux for these T dwarfs lie at longer wavelengths outside the range of our SpeX dataset. Since we cannot detect molecules such as CO or CO 2 in our SpeX data, these major absorptive species at longer wavelengths may not be properly included in these models which stretch out to 20µm. Therefore the result may be an over-estimate for the flux and thereby T eff . This estimate for T eff is then used as an input into the 1D RCE model which, while at similar T eff , may be limited by being forced onto a strict 1D RCE prescription for the PT profile. Another complicating factor is that, as shown in Figure 8 and discussed in Section 4.3.1, some objects can be offset by several hundred kelvin depending upon the choice for priors on log(g).
If we assume the offset is physically motivated, then it might be possible that these late-T dwarfs may be impacted by the presence of iron-silicate condensate clouds. To test this we calculated the cloud column optical depth, τ cld by integrating the random draws of our retrieved opacity profile from our cloud model. Figure A1 shows the result for two representative objects (all targets available in the above Zenodo link), where the most likely values of τ cld are far below unity. This suggests that optically thick clouds do not  Figure 2. Subset of the data (black round points with errors) and resulting fit (binned best fit in blue, residuals about the dashed line below in red). In blue is a higher resolution (R∼25,000) representation of the best fit spectrum. The effective temperature and gravity constraints are given in the upper right hand corner of each panel. strongly affect our dataset. Even targets such as ROSS 458C, which has been seen in previous works to better fit grid-based models that incorporate these clouds (Burgasser et al. 2010;Burningham et al. 2011;Morley et al. 2012), still strongly favors a cloud-free retrieval result.
Additionally, these targets may have additional heating mechanisms in their atmospheres that could be driving the system away from pure radiative-convective equilibrium. However it seems more likely that differences in the physical and chemical assumptions between the 1D RCE models and retrieval can bias the thermal profile using the same T eff .

Molecular Abundances
A key goal of this investigation is to develop a large, diagnostic dataset of retrieved molecular abundances across the late-T spectral class to serve as a baseline unirradiated comparison for directly imaged and transiting exoplanets. For these objects, species such as H 2 O, CH 4 , NH 3 , and K are the most readily accessible over the SpeX wavelength range. The change in the abundances of these gases with T eff are diagnostic of the driving physical and chemical processes in brown dwarf atmospheres (Burrows & Sharp 1999;Lodders & Fegley 2002;Sharp & Burrows 2007). Figure 4 shows the behavior of the retrieved molecular abundances with temperature combined with the retrieved Y-dwarf molecular abundances from Zalesky et al. (2019) (cyan). Overlaid for context are predictions from the 1D RCE models over a range of atmospheric metallicity and C/O ratios. Unless labeled otherwise each curve assumes a solar composition ([M/H]=0, C/O=0.55). To compute effective atmospheric abundances from the 1D RCE models, gas volume-mixing-ratios (VMRs) are averaged over the 1 -30bar pressure levels based on the range of pressure levels probed by the SpeX spectra. The exact retrieved abundances for all 50 objects are enumerated in the Appendix under Table 4. We find that the precision of our constraints is slightly improved relative to Zalesky et al. (2019) as expected from . TP-profiles of representative T dwarfs of our collection of 50 objects. The thermal profiles indicated here are of the retrieved objects (in red), and of grid models (in blue, see text) that they closely compared to based on their similarities in a mix of T eff , log(g), C/O, and metallicity parameters. For each of the thermals profiles, both 1 and 2σ confidence intervals are shown. Equilibrium condensation curves of Mg2SiO4, Fe, MaSiO3, Na2S, and K are shown by the cyan, blue, black, yellow, and magenta dashed lines, respectively.
the increase in quality of data. For all species we find that our abundances are constrained to within ∼0.25dex at 1σ. For both H 2 O (top left) and CH 4 (lower left) there is a trend of decreasing abundance with higher T eff , consistent with what was found in Zalesky et al. (2019). We are able to confirm that this trend extends beyond two late-T dwarfs and continues with increasing temperature though there is considerable spread. For H 2 O, this spread is consistent with what is expected for these objects between a metallicity of 0.1 to 10 times solar and various C/O ratios as shown. As there is no known physical mechanism for why these two populations should be disparate in their water content, and given the error bounds and sampling size of the Y-dwarfs, it seems likely that such a trend is a result of low-sampling from this population. Future studies which could retrieve higher precision (∼0.25dex) water abundances for brown dwarfs between 300-700K would strongly help either confirm or deny such a trend with T eff .
For CH 4 , most objects are within the predicted ranges of from the grid model. However, there are interestingly several objects with T eff between 900-1000K which deviate from this assumption. This would require either anonymously high C/O ratios (some well-above 1.0), high metallicities, or both when compared to the curves from the grid model. These objects must be interpreted with caution for two reasons. First is that only 4 objects lie above the expectations from the grid model at only 1σ. Second, we emphasize such curves do not represent all possible combinations of parameter space. Despite this, there may be physical motivation for finding objects with anomalously high CH 4 abundances near 900K. It has been well studied that L and T dwarfs have both theoretical and observational evidence for the presence of CO  2019) are not included because that work considered Na+K as one alkali species, whereas herein we are treating them as individual species. The four sub-plots show how the molecular and alkali abundances of each species varies with T eff . All the modelled molecular abundance trends are defined by log(g) = 5.0 and solar composition (M/H = 0.0, C/O = 0.5), and we assume equilibrium chemistry unless otherwise specified. We observe that the overall trends of water, methane, and ammonia are consistent with solar chemistry curve models. in their photospheres driven by disequilibrium convective mixing (Noll et al. 1997;Oppenheimer et al. 1998;Saumon et al. 2003;Golimowski et al. 2004;Burgasser et al. 2006;Geballe et al. 2009;Sorahana & Yamamura 2012;). More recently, Miles et al. (2020) found that the CO constraints for several brown dwarfs between 250-700K suggest that at lower temperatures (∼250K) there is evidence for stronger mixing (higher CO, less CH 4 ) and at higher temperatures (∼700K) there is evidence for reduced mixing strength (less CO, more CH 4 ). It is possible that we could be seeing objects that follow such a trend up to higher temperatures, however given the low number of objects in our sample and the precision of our constraints, this is largely inconclusive. There is a wealth of available early T and L spectral class objects to study to confirm such a trend, but have the added complication of obscuring condensate clouds which is beyond the scope of this work.
We find that the retrieved NH 3 abundances are largely consistent with grid model expectations, with the exception of 9 objects which show additional NH 3 near 900K. This result is surprising given that the onset of ammonia is the defining feature of the Y-dwarf spectral class. While again this could be a result of small number sampling, it is possible that this could be evidence for enhanced mixing strength at higher temperatures resulting in an excess of NH 3 compared to expectations from equilibrium. While Miles et al. (2020) found no evidence for enhanced mixing at these temperatures, disequilibrium abundances of NH 3 have been detected in the past (Saumon et al. 2007;Geballe et al. 2009;Leggett et al. 2010Leggett et al. , 2015. As was found in Part II and III, there is a strong systematic trend of decreasing potassium between the late-T to early-Y spectral types as a result of equilibrium rainout chemistry removing aluminium and silicate reserves near 1300K, leading to a delayed depletion of Na+K closer to T eff ∼700K where they condense into KCl and Na 2 S. Figure 4 shows the the trend with our larger sample which continues to indicate a steadly decreasing trend with decreasing T eff , with potassium abundances falling roughly 1 dex between -6.5 and -7.5 log(VMR) consistent with what was found in Parts II and III. We do not show the results from the Y dwarfs here as in that work, Na and K were treated as a single gas with most objects only providing upper limits.
While we are not able to fully explain the observed abundance trends with standard 1D RCE assumptions, it may indicate that some underlying physical or chemical process is not being taken into account in these models. Regardless, this sample of retrieved abundances provides a unique dataset to further explore the chemistry of sub-stellar atmospheres.  In addition to the molecular composition of a brown dwarf atmosphere, the total elemental abundance inventory, and ratios between those abundances, can assist in placing brown dwarfs in a larger context with their stellar and planetarymass cousins. Stars and field brown dwarfs are expected to form via fragmentation and core-collapse of a molecular cloud. Thus one would expect their elemental inventories to be highly similar. In contrast, planets that are formed in disks might have a completely different elemental composition due to a wide variety of processing mechanisms, differential icelines and thus elemental ratios within the disk, and/or migration of the planet within the disk (e.g. Öberg et al. 2011;). These mechanisms can change the elemental composition to range from anywhere from stellar composition (Fortney et al. 2013;Mordasini et al. 2016), to vastly different metal enrichment, and/or to have completely different C/O ratios, sometimes several factors greater or lesser than the host star (Madhusudhan et al. 2014;Eistrup et al. 2016). Identifying at what masses this dramatic compositional shift may occur for a large sample of brown dwarfs can better place brown dwarfs as a whole in a larger astronomical context.  Figure 6. log g vs. T eff for our collection of early Y dwarfs and late-T dwarfs. Here we notice that we almost have a handful of T dwarfs with low surface gravity values while Y dwarfs tend to have higher gravity values. Based on evolutionary models from (Marley et al. (2020), in prep) (shown by blue dashed curves), T dwarfs having the lowest log g values would be relatively young (∼100 Myr) while Y dwarfs having the highest log g would be relatively old (∼12 Gyr). Most of the objects above the 10 Gyr curve are considered atypical although their high gravities can be justified by the logic that under-estimations of gravities that are done in grid models (Schneider et al. 2015;Zalesky et al. 2019).

Metallicity and Elemental Ratios
As noted in previous studies in this series, the bulk of the metals found in brown dwarf atmospheres with T eff ≤1000K are largely comprised of C and O contained in the water and methane content with a smaller amount of N stored in NH 3 . The metallicity is computed as, where the metallicity of the T dwarf is taken to be the summation of the number of elemental species included in the retrieval model. The C/O ratio is taken as, (2) Figure 5 shows the [M/H] compared against our computed C/O for both this study and those objects in Zalesky et al. (2019). For context a representative sample of local FGK stars are overlaid (dark red) (Hinkel et al. 2014). As in Parts II and III, we note that we additionally weight our water abundance by a factor of 1.3 to account of the loss of oxygen to condensate rainout species such as enstatite and forsterite (Fegley & Lodders 1994). We find that, in agreement with previous papers in this series, the broad-scale metallicities and C/O ratios for our brown dwarfs largely match both in value and in distribution that of the local FGK stellar population. The new population (largely green-orange color) seems to fill in the gap between our two previous late-T and early-Y populations. Additionally, we performed a two-dimensional, two-sample, Kolmogorov-Smirnov test and found a p-value of 0.18 between the brown dwarf and stellar populations (Press et al. 1992). This suggests that our measurements display no statistically significant difference from the local FGK population. Though there is considerable scatter, there is a slight positive correlation between metallicity and T eff .

T eff & Gravity
A brown dwarf's effective temperature and gravity are highly diagnostic of both its age and evolutionary history (Burrows et al. 2001;Baraffe et al. 2003;Saumon & Marley 2008). Given that our sample should be comprised of largely field-age brown dwarfs, deviations in the these parameters from evolutionary models serves as both a sanity-check for our retrieval model, but may also serve as a sign of youth if a particularly low surface gravity is found. Figure 6 shows our sample of late-T objects (orange) against results from Zalesky et al. (2019) (cyan). Overplotted are evolutionary cooling curves for several different ages (blue-dashed, Marley et al. 2021, in prep). A large and extensive discussion for why the population from Zalesky et al. (2019) have anomalously high gravities is given in that paper.
Here we focus on our sample, which for a roughly forty objects agree with being near 1-3Gyr at 1σ. However, seven T dwarfs appear to have lower surface gravity values compared to the rest of the objects by several sigma with ages of near 0.1Gyr. Additionally several objects are within 1σ of being consistent with a 10Gyr age. Both of these ages are unusual given our sample and necessitate a sanity check of our retrieval model.
The seven youngest objects are WISE J1254, PSO-J224, ROSS-458C, WISE-J1257, WISE-J1322, WISE-J1959, and WISE J0325. All the aforementioned T dwarfs, with the exception of ROSS-458C, have relatively low water and methane abundances. Furthermore, all seven of these T dwarfs have very low ammonia abundances. This observation of brown dwarfs exhibiting low gravities and metallicites is highly atypical (Helling & Casewell 2014;Yuan et al. 2011;Helling & Casewell 2014). To test the credibility of the surface gravity constraint, we conducted an experiment on both the lowest gravity object (WISE 0325) and on one of our two highest gravity objects (2MASS 0729) by varying their gravity priors.

Metallicity -Surface Gravity Relationship
In this test we took a representative object that had both a relatively high and low log(g) constraint. The low-gravity object, WISE 0325, was used to perform a retrieval with a "full" log(g) prior range of 0 ≤ log(g) ≤ 6 compared to a retrieval done with a "constrained" log(g) prior range of 4.3 ≤ log(g) ≤ 6. The high-gravity object, 2MASS 0729, was set to 0 ≤ log(g) ≤ 4.3 while the full range was maintained at 0 ≤ log(g) ≤ 6. In the constrained cases, the gravity ranges were deliberately restricted to lower and higher gravity ranges to see how the retrieval results changed if forced to be in agreement with the system's likely true age. Figure 7 shows stair-pairs plots for both object's molecular abundance and gravity constraints in both the constrained (red) and unconstrained (blue) test cases. The constraint on log(g) shows a strong several-sigma offset between the two cases and importantly that the posterior in the constrained case is highly non-Gaussian as it is pushing up against the limit of the prior. This indicates that, despite being "guided" to a lower/higher gravity in-line with field-age expectations, the spectrum itself has features which favor a higher/lower gravity than the priors allow.
Alterations in gravity prior ranges also impact the posteriors of other species like water, ammonia, and methane as the parts of the spectra where gravity would be most constraining are the same regions where the absorption of these gases are most prominent. Potassium notably appears to be insignificantly affected by changes in gravity. This is likely because potassium absorption features only affect a small portion of the spectrum between 0.95-2.5 µm, mostly in the Y and J bandpasses. This high correlation between gravity and metal- Figure 7. Corner plot posteriors of 2MASS-J072290002 (highest gravity object) and WISE-J032547.72 (lowest gravity object). This is based on a test of constraining priori ranges of gravity for the two objects. Gravity range restrictions affect gravity's correlation with the other gas species. Limited logg priors' results for surface gravity lean more toward full log(g) priors' results. licity also impacts the constraints on the thermal profiles as well. as shown in Figure 8.
Though a large majority of our objects are within reasonable expectations for field age brown dwarfs, the age estimates from several of these objects are highly nonphysical being either less than 0.1Gyr or older than the age of the Universe in the most extreme cases. These results show that retrievals, 1D RCE models, and comparison/sanity checks be-tween both are needed to fully understand the atmospheres of brown dwarfs. One area this is possible is with different types of benchmark systems, where either metallicity or age is constrained. In our sample we have three systems which are in a binary (ROSS 458B,HD 3651B) or quaternary (Gliese 570D), all of which are within 1σ for their stellar companion's age and metallicity, giving confidence to the reliability of our retrieval results despite the many caveats of using this technique.

Summary and Conclusions
Below we list the key conclusions of our analysis: 1. For most of our sample, we find that the retrieved thermal structures are in agreement with with predictions from self-consistent grid based approaches (Figure A3). However, 18 objects have retrieved thermal profiles that are much hotter, at roughly 2σ, than the grid models. Though possibly physical in origin, it is more likely owed to the different physical and chemical assumptions in each model. Longer wavelength information, namely beyond 2.5 microns where most of the flux is, would significantly help in testing if this is the true cause of the discrepancy.
2. In addition to the thermal profiles, the highlight of this work is the retrieval of abundance constraints for the major atmospheric molecular species [H 2 O, CH 4 , NH 3 , K] with a precision of roughly 0.25dex at 1σ (Figure 4, Table 4). For a majority of our sample, we find that H 2 O, CH 4 , and NH 3 are consistent with expectations from equilibrium chemistry. There are however, a small subset of objects which do not fit these trends. There are four objects with elevated CH 4 and roughly ten with elevated NH 3 constraints. Though it may be a simple result of small number sampling, there is literature showing that objects near this temperature range may begin to show signs of elevated CH 4 and NH 3 due to non-equilibrium vertical mixing. Developing a larger sample of earlier T dwarfs (T eff > 900K) requires accurately taking into consideration a realistic cloud model into our retrieval framework, but would show if this trend continues into higher temperatures. Finally we note that our sample helps confirm previous findings that the depletion of Na and K near 500-600K owed to the chemical rainout of species like enstatite and forsterite is present.
3. Using the retrieved molecular abundances we are able to infer the atmospheric metallicity and C/O ratio. We find that the distribution of our 50 targets is broadly consistent with measurements of the local stellar FGK population. Though there is substantial scatter, there also is a slight positive correlation between atmospheric metallicity and effective temperature.
4. We compare our constraints on both T eff and log(g) to expectations from evolutionary models and previous results of Y dwarfs (Figure 6). For a majority of objects, we find them to be consistent with expectations of typical field age brown dwarfs. However we also find several outliers that have unphysically old and young ages. We preformed several tests in which we altered the priors on our retrieval model to "guide" these objects to more physically realistic ages. Despite these efforts, we found that the resulting posteriors were highly non-Gaussian and were clustered near the limit of the new prior ranges. This indicates that there is some as yet undetermined feature in the spectra which favors these anomalous ages. Given that these age estimate are likely incorrect, this result highlights the importance of complimenting such retrieval studies with more traditional grid models.
This uniform atmospheric retrieval analysis of 50 T dwarfs provides a large and unique dataset of retrieved thermal profiles, gravities, and atmospheric abundances which offer insight into the physical and chemical mechanisms at work in their atmospheres. This dataset now can serve as a baseline for comparisons between other stars and planets to place brown dwarfs in a holistic astronomical context.   Table 4 continued  Figure A2. A continuation of Figure 2 showing the SpeX spectra (black), best-fit spectra (blue) and residuals (red) for all objects in our sample.  Figure A3. A continuation of Figure 3 showing the retrieved (red) and grid-model (blue) PT profiles, overlaid with relevant condensate curves, for all objects in our sample.