Retrieval of the d/sdL7+T7.5p Binary SDSS J1416+1348AB

We present the distance-calibrated spectral energy distribution (SED) of the d/sdL7 SDSS J14162408+1348263A (J1416A) and an updated SED for SDSS J14162408+1348263B (J1416B). We also present the first retrieval analysis of J1416A using the Brewster retrieval code base and the second retrieval of J1416B. We find that the primary is best fit by a nongray cloud opacity with a power-law wavelength dependence but is indistinguishable between the type of cloud parameterization. J1416B is best fit by a cloud-free model, consistent with the results from Line et al. Most fundamental parameters derived via SEDs and retrievals are consistent within 1σ for both J1416A and J1416B. The exceptions include the radius of J1416A, where the retrieved radius is smaller than the evolutionary model-based radius from the SED for the deck cloud model, and the bolometric luminosity, which is consistent within 2.5σ for both cloud models. The pair’s metallicity and carbon-to-oxygen ratio point toward formation and evolution as a system. By comparing the retrieved alkali abundances while using two opacity models, we are able to evaluate how the opacities behave for the L and T dwarf. Lastly, we find that relatively small changes in composition can drive major observable differences for lower-temperature objects.


Introduction
Brown dwarfs are a class of astronomical objects that straddle the mass boundary between stars and planets with masses 75 M Jup Chabrier & Baraffe 1997) and effective temperatures of 250-3000 K, corresponding to late-type M, L, T, or Y spectral types (Burgasser et al. 2002;Kirkpatrick 2005;Cushing et al. 2011). Due to electron degeneracy, they never reach a core temperature high enough for stable hydrogen burning but instead contract and cool through their lifetimes, progressing through spectral classifications as they age.
Field-age brown dwarfs anchor the spectral type scheme; however, low-gravity, low-metallicity, and color outliers expand the standard scheme. Low-metallicity sources, known as subdwarfs, have unusually blue near-infrared (NIR) J−K colors (Burgasser et al. 2003) compared with equivalent field sources. Spectral features distinguishing them from field dwarfs include enhanced metal-hydride absorption bands (e.g., FeH), weak or absent metal oxides (TiO, CO, VO), and enhanced collisionally induced H 2 absorption (Burgasser et al. 2003 and references therein). Subdwarfs also exhibit substantial radial velocities, high proper motions, and inclined, eccentric, and sometimes retrograde Galactic orbits indicating membership in the Galactic halo (Burgasser et al. 2008;Dahn et al. 2008;Cushing et al. 2009). To date, as classified by Zhang et al. (2017Zhang et al. ( , 2018bZhang et al. ( , 2018aZhang et al. ( , 2019, there are approximately 66 L subdwarfs and 41 T subdwarfs, although most T subdwarfs are not classified as such in previous literature (see Table 3 of Zhang et al. 2019). To be identified as a T subdwarf in Zhang et al. (2019), T dwarfs need to have a suppressed K-band spectrum.
Presently there is only one subdwarf L+T system, SDSS J14162408+1348263AB (hereafter J1416AB), and it is ideally suited for low-metallicity bd-bd binary atmospheric characterization via retrievals. In this paper, we determine and examine fundamental parameters and atmospheric features of J1416AB via two methods: (1) by coupling the empirical bolometric luminosity, from the distance-calibrated spectral energy distribution (SED), with evolutionary models and (2) by atmospheric retrievals where we explore similarities and differences between the pair to determine their formation and evolution and to understand their individual atmospheric structure.
In Section 2 we present literature data on J1416AB. Section 3 presents the data used for creating distance-calibrated SEDs and the retrievals, as well as the resultant fundamental parameters derived from creating the SED. Section 4 describes our retrieval framework and setup for J1416AB. Retrieval results for J1416A and J1416B are discussed in Sections 5 and 6, respectively. Fundamental parameters derived from SED and retrieval methods are compared to the literature and evolutionary models in Section 7. Lastly, Section 8 brings together the individual retrievals of J1416AB to discuss the alkali abundance, metallicity, and carbon-to-oxygen (C/O) ratios derived and what we can interpret for the system as a whole.

Literature Data on SDSS J1416AB
At the time of discovery, J1416AB was one of the few known widely separated L+T systems, thus allowing for the properties of both to be examined in tandem. This system is a benchmark, as features of the primary indicate an old age for the system. Here we present the literature data for the independent discoveries of the L and T dwarfs.
Many studies aimed to determine the fundamental properties of J1416A by fitting its spectrum to self-consistent grid models (Burningham et al. 2010;Bowler et al. 2010;Cushing et al. 2010;Schmidt et al. 2010). Its atmosphere was determined to be relatively dust-free (Burningham et al. 2010), and like other blue L dwarfs it possibly had a thin or patchy cloud deck with large grains that could cause the observed blue NIR colors. Additionally, J1416A might have an older age and higher surface gravity (Bowler et al. 2010;Schmidt et al. 2010). Cushing et al. (2010) found evidence for vertical mixing in the atmosphere due to the lack of CH 4 absorption at 3.3 μm. Temperature estimates of J1416A vary from 1500 to 2200 K (Burningham et al. 2010;Bowler et al. 2010;Cushing et al. 2010;Schmidt et al. 2010), while the literature agrees on a surface gravity of 5.5 dex (Burningham et al. 2010;Bowler et al. 2010;Cushing et al. 2010) and a weakly or unconstrained age (Burningham et al. 2010;Bowler et al. 2010;Schmidt et al. 2010).
J1416A was examined for variability in Khandrika et al. (2013), Metchev et al. (2015), and Miles-Páez et al. (2017). Khandrika et al. (2013) found marginal evidence of variability detected in one night of their observations using Gemini camera J and K′ bands on the Shane telescope. Metchev et al. (2015) monitored J1416A using Spitzer ch1 (14 hours) and ch2 (7 hours) as part of their Weather on other Worlds survey to look for variability attributed to patchy clouds, finding no evidence for variability. In Miles-Páez et al. (2017), variability correlated to activity was tested using the Gemini Multi-object Spectrograph (GMOS-N) with the R831-G5302 grating, but no evidence for variability was found.
Values for J1416A from the literature and those determined in this work are listed in Table 1. All literature values are also listed in Table 9 for comparison in Section 7.

Literature Data on SDSS J1416B
ULAS J141623.94+1348836.30 (hereafter J1416B) was discovered by Burningham et al. (2010) through a cross match of SDSS and UKIRT, finding a separation of 9″ between the A and B component. J1416B was also independently discovered by Scholz (2010) with a projected separation of 75 au, which we have updated (now 83.7 au) using the Gaia DR2 parallax and the angular separation from Burningham et al. (2010). Like J1416A, J1416B has unusual features of a late-T dwarf. In particular, the CH 4 − J-early peculiarity (where the CH 4 -J index on the red side of the J-band peak suggests an earlier spectral type than the H 2 O-J index on the blue side of the Jband peak), the very blue H−K color, and the extremely red H−[4.5] color lead to its classification as a T7.5p (Burningham et al. 2010). At the time of its discovery, J1416B was both the bluest H−K and reddest H−[4.5] T dwarf. The CH 4 −J-early peculiarity of 1416B pointed toward either low metallicity or high surface gravity (Burningham et al. 2010). It was noted that J1416B forms a sequence with other lowmetallicity and high-gravity T dwarfs, and because of the extremely red H−[4.5] color it could not be ruled out as a binary itself (Burningham et al. 2010). Burgasser et al. (2010a) classified J1416B as a T7.5 but noted strong water and methane bands, a possible detection of ammonia between 1 and 1.3 μm, and a broadened Y-band peak and suppressed K band indicative of high gravity or low metallicity in its spectrum. Kirkpatrick et al. (2016) regarded J1416B as an sdT7.5 in relation to J1416A, andZhang et al. (2017) also classified it as sdT7.5 via their subdwarf metallicity classification scheme. Presently, J1416B has three NIR prism spectra (IRCS, SpeX, and FIRE) from Burningham et al. (2010) and Burgasser et al. (2010aBurgasser et al. ( , 2010b, respectively. Fundamental parameters of J1416B were determined through comparison with grid models (Burgasser et al. 2010a(Burgasser et al. , 2010b, SED fitting (Filippazzo et al. 2015), and atmospheric retrieval . Burgasser et al. (2010a) found that J1416B was well matched to the archetype blue T dwarf 2MASS J09393548−2448279. They determined a T eff =650±50 K, log g=5.2±0.4, [Fe/H]−0.3, and K zz =10 4 using the Saumon & Marley (2008) models and used Baraffe evolutionary models to find an age range of 2-10 Gyr, mass between 22 and 47 M Jup , and radius of 0.83R Jup . Both cloudless and cloudy models were fit to the spectrum of J1416B in Burgasser et al. (2010b), with cloudy   Metchev et al. (2015) for variability with no evidence found in Spitzer ch1 and ch2. All literature values are also listed in Table 10 for comparison in Section 7.

Data Used and Results from Generating the SED
The fundamental parameters for J1416AB were determined using the technique of Filippazzo et al. (2015), where we create a distance-calibrated SED using the spectra, photometry, and parallax. 13 The SED of J1416A uses the SpeX short-crossdispersed (SXD) and long-cross-dispersed (LXD) spectrum from Cushing et al. (2010), while J1416B uses the SpeX prism spectrum from Burgasser et al. (2010a). The photometry and Gaia parallax used for both sources are listed in Table 1. Table 2 lists the spectra used in the SEDs and the retrieval models, which differ for J1416A due to the current time constraints on data resolution for our retrieval model.
To generate the SED of J1416A, the SpeX SXD and LXD spectra were stitched, linearly interpolating to fill gaps in the data, into a composite spectrum and then scaled to the absolute magnitudes of the observed photometry. For J1416B we scale the SpeX prism spectrum to the absolute magnitudes of observed and synthetic (those calculated based on empirical relations) photometry. Synthetic photometry for J1416B is included because if we linearly interpolated between W2 and W3, without including the synthetic MIR IRAC Ch3 and Ch4 photometry calibrated based on field dwarfs, we would likely overestimate the mid-infrared flux compared to most T dwarfs, causing a noticeable change in the T eff . As there are no known low-metallicity T dwarfs with IRAC Ch3 or Ch4 MIR photometry, we cannot place a level of error on their difference from field T dwarfs. The SEDs of J1416A and J1416B are shown in Figure 1, with the synthetic magnitudes used for J1416B plotted as transparent squares in Figure 1 The bolometric luminosity (L bol ) was determined by integrating under the distance-calibrated SED from 0 to 1000 μm, using a distance of 9.3±0.3 pc based on the Gaia Collaboration et al. (2018) parallax measurement. The effective temperature (T eff ) was calculated using the Stefan-Boltzmann law with the resultant inferred radius from the cloudless Saumon & Marley (2008) low-metallicity (−0.3 dex) evolutionary model. The low-metallicity models were chosen for the assumed radius due to the literature spectral type classification of sd for both components. Additionally, as done in Filippazzo et al. (2015), the Chabrier et al. (2000), Baraffe et al. (2003), and cloud-free Saumon & Marley (2008) evolutionary models were also used to determine the radius. The final radius range was set as the maximum and minimum from all model predictions as done in Filippazzo et al. (2015). An age range of 0.5-10 Gyr for the system was chosen to conservatively encompass possible field and subdwarf ages. Additional details on the SED generation can be found in Filippazzo et al. (2015). Fundamental parameters derived for J1416A and J1416B using this approach are listed in Table 1 and are compared to the literature in Section 7 (also see Tables 9 and 10).

The Brewster Retrieval Framework
Our retrievals use the Brewster framework (Burningham et al. 2017) with a modified setup from the one in Burningham et al. (2017) in order to optimize for low-metallicity atmospheres. A summary of the Brewster framework with our modifications is discussed below. We differ from Burningham et al. (2017) with a higher resolution for opacity sampling, using a second method (thermochemical equilibrium with rainout) for determining gas abundances, and expanded temperature and mass priors. A more detailed description of Brewster can be found in Burningham et al. (2017).   13 SEDkit is available on GitHub athttps://github.com/hover2pi/SEDkit. The Eileen branch was used for this work.

The Forward Model
The forward model in Brewster uses the two-stream radiative transfer technique of Toon et al. (1989), including scattering, as first introduced by McKay et al. (1989) and subsequently used by, for example, Marley et al. (1996), Saumon & Marley (2008), and Morley et al. (2012). We use a 64 pressure layer (65 levels) atmosphere with geometric mean pressures between = -P log 4 and 2.3 bars in 0.1 dex spaced intervals. The temperature in each layer is characterized by the three exponential functions as done following the Madhusudhan & Seager (2009) parameterization, splitting the atmosphere in three zones where the pressure and temperature are related by < < < < > = a a --P P P P e P P P P e P P T T : Zone 1 , : Zone 2 , : Zone 3 , where P 0 and T 0 are the pressure and temperature at the top of the atmosphere and the atmosphere becomes isothermal at pressure P 3 with temperature T 3 . Since P 0 is fixed in our model and continuity at the zonal boundaries requires fixing two parameters, we consider six free parameters: α 1 , α 2 , P 1 , P 2 , P 3 , and T 3 . A thermal inversion can occur when P 2 >P1; however, this is ruled out by setting P 2 =P 1 , thus further simplifying this to five free parameters.

Gas Opacities
Layer optical depths due to absorbing gases are calculated using opacities sampled at a resolving power R=10,000 taken from Freedman et al. (2008Freedman et al. ( , 2014. Line wing profiles based on the unified line shape theory (Allard et al. 2007a(Allard et al. , 2007b) are used to account for the broadening of the D resonance doublets of Na I (∼0.59 μm) and K I (∼0.77 μm) in brown dwarf spectra. Tabulated line profiles (Allard N., private communication) are calculated for the Na I and K I D1 and D2 lines broadened by collisions with H 2 and He, for temperatures in the 500-3000 K range and perturber (H 2 or He) densities up to 10 20 cm −3 with two collisional geometries considered for broadening by H 2 . Within 20 cm −1 of the line center there is a Lorentzian profile with a width calculated from the same theory. While there are updated versions of these opacities (Allard et al. 2016(Allard et al. , 2019Phillips et al. 2020), we did not have access to them for this work. We also use the Na and K alkali opacities from Burrows & Volobuyev (2003) to be consistent with Line et al. (2017) in the J1416B retrievals.
Across our temperature-pressure regime, the line opacities are tabulated in 0.5 dex steps for pressure and in steps ranging from 20 to 500 K as we move from 75 to 4000 K in temperature where we then linearly interpolate this to our working pressure grid. We include free-free continuum opacities for H − and -H 2 and bound-free continuum opacity for H − , which are influenced by the H − metallicity and determined from the thermochemical equilibrium grid (see Section 4.3). Continuum opacities for H 2 -H 2 and H 2 -He collisionally induced absorption, using cross sections from Richard et al. (2012) and Saumon et al. (2012), are included, as well as Rayleigh scattering due to H 2 , He, and CH 4 , but we neglect the remaining gases. Neutral H gas fraction abundance was determined from the thermochemical equilibrium grid. The atmosphere is assumed to be dominated by H 2 and He, with proportions of (0.84H 2 + 0.16He) based on solar abundances.
After including the retrieved gases, neutral H, H − , and electrons H 2 and He are assumed to make up the remainder of the gas in a layer. The former is drawn from the thermochemical equilibrium grids discussed later in this section.

Determining Gas Abundances
As done in Burningham et al. (2017), we use the uniformwith-altitude mixing ratios method for absorbing gases and retrieve these directly, also known as "free" retrievals, for all of our retrieval models. While simple, the uniform-with-altitude mixing method cannot capture important variations in gas abundance with altitude for some species (i.e., see metal oxides and metal hydrides of J1416A and the alkalies for J1416B), which can vary by several orders of magnitude in the photosphere and are expected to make a large contribution to the flux we observe. Freely retrieving abundances that vary with altitude would be preferred; however, the resultant large number of parameters to solve for in this approach is computationally difficult. To address this issue we use a second method, the chemical equilibrium method, which instead retrieves [Fe/H] and C/O. Gas fractions in each layer of this method are pulled from tables of thermochemical equilibrium abundances as a function of T, P, [Fe/H], and C/O ratio along with the thermal profile of a given state vector. The thermochemical equilibrium grids we use were calculated using the NASA Gibbs minimization CEA code (McBride & Gordon 1994), based on previous thermochemical models (Fegley & Lodders 1994Lodders 1999Lodders , 2002Lodders & Fegley 2002Visscher et al. 2006Visscher et al. , 2010Lodders 2010;Visscher 2012;Moses et al. 2012Moses et al. , 2013 and recently utilized to explore gas and condensate chemistry over a range of conditions in substellar atmospheres (Morley et al. 2012(Morley et al. , 2013Skemer et al. 2016;Kataria et al. 2016;Wakeford et al. 2017). The chemical grids in this work determine equilibrium abundances of atmospheric species over pressures ranging from 1 microbar to 300 bars, temperatures between 300 and 4000 K, metallicities in the range −1.0<[Fe/ H]<+2.0, and C/O abundance ratios of 0.25-2.5 times the solar abundance.

Cloud Model
The cloud model follows that of Burningham et al. (2017), with options for a "deck" or "slab" cloud parameterization. Both clouds are defined similarly, where the cloud's opacity is distributed among layers in pressure space, with the optical depth either gray or as a power law (τ=τ 0 λ α , where τ 0 is the optical depth at 1 μm).
The deck cloud is parameterized by (1) a cloud top pressure, P top , the point at which the cloud passes τ=1 (looking down); (2) the decay height, D P log , over which the optical depth falls to lower pressures as t log log ( ( )) ( ); and (3) the cloud particle single-scattering albedo. The deck cloud becomes optically thick at P top . At P>P top , the optical depth increases following the decay function until it reaches Δτ layer =100. With this decay function, the deck cloud can quickly become opaque with increasing pressure, and therefore we obtain essentially no atmospheric information from deep below the cloud top. Because of this, it is important to note that the pressure-temperature (PT) profile (and spread) below the deck is an extension of the gradient (and spread) at the cloud top pressure.
Unlike the deck cloud, it is possible to see the bottom of the slab cloud and thus include an additional parameter for determining the total optical depth at 1 μm (τ cloud ), bringing the total number of parameters for the slab cloud to 4. The optical depth is distributed through the slab cloud extent as dτ/ dP∝P (looking down), reaching its total value at the bottom (highest pressure) of the slab. In principle, the slab can have any optical depth; however, we restrict our prior as 0.0 τ cloud 100.0. Because it is possible to see to the bottom of the slab cloud, a physical extent in log-pressure (D P log ) is determined, instead of the decay scale, as was done for the deck cloud.
If the deck or slab cloud is nongray, an additional parameter for the power (α) in the optical depth is included.

Retrieval Model
As described in Burningham et al. (2017), we use EMCEE (Foreman-Mackey et al. 2013) to sample posterior probabilities. Table 3 shows our priors for both J1416A and J1416B. We differ from the Burningham et al. (2017) setup by extending the thermal profile temperature up to 6000 K for both J1416A and J1416B and extending the mass prior up to 100M Jup for J1416A (up to only 80 M Jup for J1416B) to expand the surface gravity in an effort to encompass likely ranges for subdwarfs. In our retrievals of J1416A and J1416B we use their distance-calibrated SpeX prism spectra (output from generating our SED) trimmed to the 1.0-2.5 μm region and set the distance to 10 pc. This spectrum calibration differs from Burningham et al. (2017), where they calibrated the spectrum to the 2MASS J-band photometry and used the true distance in their initialization.
For each retrieval of J1416A and J1416B we initialize 16 walkers per parameter in a tight Gaussian for the gases, surface gravity, wavelength shift between the model and data (Δλ), and scale factor where R≈1.0 R Jup . Gases are centered around the approximate solar composition equilibrium chemistry values for gas volume mixing ratios, while the surface gravity is initiated centered around the SED-derived value. The tolerance parameter has a flat distribution across the entire prior range. For cloud parameters, the cloud top pressure and power law are initialized as tight Gaussians, while the optical depth, albedo, and cloud thickness are flat across the entire prior range. As in Burningham et al. (2017), we use the fiveparameter thermal profile, as we do not expect a temperature inversion for either of these objects, and use the Saumon & Marley (2008) T eff =1700 K log g=5.0 model to initialize α 1 , α 2 , P 1 , P 2 , P 3 , and T 3 for both J1416A and J1416B. Differences in the individual setups between J1416A and J1416B are discussed in the following subsection.

J1416A
To explore the atmosphere of J1416A, we retrieved for the following gases: H 2 O, CO, CO 2 , CH 4 , TiO, VO, CrH, FeH, K, and Na. As done in Burningham et al. (2017) and Line et al. (2015), we tie K and Na together as a single element in the state vector assuming a solar ratio taken from Asplund et al. (2009). Additionally, we include the H − bound-free and free-free continuum opacities to account for the possibility of the profile going above 3000 K in the photosphere. As stated above, the log g mass prior ranges from 0 to 100 M Jup . The multiple cloud parameterizations are tested building up from the cloudless to the four-parameter power-law slab cloud model. We also test both the uniform-with-altitude mixing ratios and chemical equilibrium methods for determining the gas abundances.

J1416B
The retrieval setup and initialization for 1416B is similar to J1416A with the following exceptions: (1) as J1416B is much cooler, we retrieve only H 2 O, CH 4 , CO 2 , NH 3 , K, and Na (where Na and K are tied together), and (2) we do not include the H − bound-free and free-free continuum opacities as the profile is cooler than the L dwarf and does not warrant them. As the T dwarf should be less massive, the log g mass prior ranges from 0 to 80 M Jup . We also differ from the retrieval setup of that in Line et al. (2017) by (1) excluding CO 2 and H 2 S in our gas list as Line et al. (2017) could only derive upper limits and (2) testing both the Allard and Burrows alkali opacities.

Model Selection
A variety of parameters were tested in our retrievals of J1416A and J1416B, with some aspects remaining constant throughout (the gases included in each model), while others differed. The aspects that were allowed to differ in our retrievals include cloud parameterization, gas abundance method, and alkali opacities. To compare all of our retrievals, model selection was assessed using the Bayesian information criterion (BIC) where the lowest BIC is preferred. We use the following intervals from Kass & Raftery (1995) for selecting between two models, with evidence against the higher BIC as 1. 0<ΔBIC<2: no preference worth mentioning 2. 2<ΔBIC<6: positive 3. 6<ΔBIC<10: strong Notes. a Gravity prior upper limit only to 80M Jup for J1416B. b For the deck cloud, this is the pressure where τ cloud =1; for a slab cloud, this is the top of the slab. c Decay height for deck cloud above the τ cloud =1.0 level. d Thickness and τ cloud only retrieved for slab cloud.

ΔBIC>10: very strong
A variety of cloud assumptions are explored in our retrievals by building up from the least complex cloud-free model to the most complex slab cloud model. Prior to moving from the cloud-free to cloudy models, we tested the impact of assuming different metallicities when determining the neutral H, H − , and electron abundances used for the continuum opacity calculations as both targets are expected to be low metallicity. We found using low-metallicity ([M/H]=0.3) ion fractions to be indistinguishable from the solar-metallicity ion fractions and thus proceeded using the solar ion abundances for the cloudy models.
Once the "winning" model was determined, we tested two additional methods for calculating gas abundances: (1) the thermochemical equilibrium assumption and (2) alternate opacities based on the Burrows and Allard line-broadening treatments (Burrows for J1416A and Allard for J1416B) in the uniform-with-altitude assumption. We examined both Allard and Burrows alkali opacities as there is no agreement in the literature as to which is the preferred choice in retrievals or grid models (Saumon & Marley 2008;Todorov et al. 2016;Burningham et al. 2017;Line et al. 2017;Gravity Collaboration et al. 2020, Marley et al. 2020. As done in Burningham et al. (2017), we started with the Allard opacities for the L dwarf, and as done in Line et al. (2017), we started with the Burrows opacities for the T dwarf. By testing the alternative line profile treatments, we aim to establish the impact of this choice on the derived alkali abundances for J1416AB.

Retrieval Model of J1416A
The ΔBIC for all tested models for J1416A are shown in Table 4. The best-fitting model is parameterized as a power-law deck cloud. However, this model is indistinguishable from the power-law slab cloud (ΔBIC=1.40), meaning both models provide similarly good fits to the spectroscopic features observed in J1416A. The retrieval results of the power-law deck and power-law slab cloud models are discussed in Sections 5.1 and 5.2, respectively. The "winning" deck and slab cloud models were also indistinguishable when using the Burrows alkali opacities instead. Section 8.1 provides further discussion of the preferred choice of alkali opacities for comparing J1416A to J1416B.  deeper. However, one should note that our deep PT profile (below photosphere) is an extrapolation of the shape at lower pressures as there is little contribution to the observed flux. At pressures lower than 1 bar (higher up in the atmosphere), the median PT profile is hotter than the Sonora  models, which was also seen for two L dwarfs in Burningham et al. (2017). The median deck cloud, shown in the center of Figure 2(a), becomes optically thick deeper than ∼10 bars with the cloud top location in pressure space quite tightly constrained to log = -+ P 1.14 0.21 0.18 bars. However, the extent of the cloud (gradient region) where the optical depth falls to τ=1/2 (dashed black line) is poorly constrained. Figure 2(b) shows the contribution function for this model along with the τ=1 gas and cloud contributions. The contribution function in a layer is defined as where B(λ, T(P)) is the Planck function, zero is the pressure at the top of the atmosphere, P 1 is the pressure at the top of the layer, and P 2 is the pressure at the bottom of the layer. The majority of the flux contributing to the observed spectrum of J1416A comes from the approximately 1 to 18 bar region, corresponding to the photosphere. The observed flux in the Y band is dominated by the gas at shorter wavelengths (1.11 μm), while the cloud opacity dominates from ∼1.06 to 1.11 μm. The J band is shaped by the gas opacity, with the cloud opacity sitting just below the τ=1 gas line, potentially contributing minor amounts of opacity. In the H and K bands, the gas opacity dominates our observed flux as it becomes optically thick well before (higher up) the cloud contribution. The lack of the cloud's contribution to the J band is a possible factor for J1416A's observed unusually blue J−K color of 1.03±0.03. In Figure 3 we compare the contribution functions for J1416A with the two L dwarfs in Burningham et al. (2017), 2MASS J05002100+033050 (hereafter J0500 +0330), and 2MASSW J2224438−015852 (hereafter J2224 −0158). We can see that the median τ cloud =1 level is reached at shallower pressures for both comparison targets and lies above the median τ gas =1 level in most of the Y and the entire J band for J0500+0330 and the entire Y and J bands in the case of J2224−0158. This points toward seeing deeper into the atmosphere at the J band of J1416A potentially due to its lower metallicity, than the field source J0500+0330 and the red L dwarf J2224−0158, as the possible cause of the observed blue J−K color.  Table 5 for ease of reading. An extrapolated value for L bol is not shown in Figure 4 as L bol showed no interesting correlations with any parameter. Our retrieved gas abundances are compared to values expected from chemical equilibrium grids in Section 5.1.4.

Retrieved Gas Abundances and Derived Properties
The derived radius and mass are calculated from the retrieved scaling factor (R 2 /D 2 ) and log g values, along with the parallax measurement. To derive the T eff , we use the radius and integrate the flux in the resultant forward model spectrum across 0.6-20 μm. Our retrieval-derived T eff is ∼200 K hotter than our semiempirical K). This is due to the retrieval-based radius being 0.12R Jup smaller than the model radius from the SED method. Our retrieved gravity and extrapolated mass agree within 1σ to the gravity and mass we derive from evolutionary models when generating the SED (retrieval: log = To derive the C/O ratio we exclude all carbon-and oxygenbearing molecules that are not constrained for both cloud models of J1416A, thus assuming all of the carbon exists in CO where f H 2 is the H 2 fraction, N H is the number of neutral hydrogen atoms, N element is the number atoms for the element of interest (C, O, V, Cr, Fe, and Na+K), n atom is the number of atoms of that element in a molecule (e.g., two for oxygen in Figure 4. J1416A power-law deck cloud posterior probability distributions for the retrieved parameters and extrapolated parameters. One-dimensional histograms of the marginalized posteriors are shown along the diagonals, with 2D histograms showing the correlations between the parameters. The dashed lines in the 1D histograms represent the 16th, 50th, and 84th percentiles, with the 68% confidence interval as the width between the 16th and 84th percentiles. Parameter values listed above are shown as the median± 1σ. Gas abundances are displayed as log 10 (X) values, where X is the gas. T eff , radius, mass, C/O ratio, and [M/H] are not directly retrieved parameters but are calculated using the retrieved R 2 /D 2 and log g values along with the predicted spectrum. Our derived C/O ratio is absolute, where solar C/O is 0.55, while our [M/H] is relative to solar. Values for CO 2 , CH 4 , and TiO are not constrained and thus only provide upper limits.
CO 2 ), f gases is the total gas fraction containing only the constrained gases, and N tot is the total number of gas molecules. ). We note that for both C/O and [M/H] it does not matter if we include or exclude VO, which is done when comparing to J1416B; the ratios agree within 1σof each other. This C/O ratio does not account for oxygen lost to silicate formation, which we address in further detail in Section 8.2. By examining the overplotted condensation curves on the PT profile ( Figure 2) to identify the possible cloud deck species, we find that no condensation curves intersect the profile at the cloud top location. Burningham et al. (2017) found iron or corundum as likely cloud compositions for their two L dwarfs as these condensation curves intersected the PT profile at the top of the deck cloud. Thus for J1416A, iron or corundum (Al 2 O 3 ) could be possible deck cloud candidates; however, the cloud optical depth continues to increase beneath the phaseequilibrium condensation point on our thermal profile. This could be due to cloud opacity deriving from the condensation of other species at deeper layers or opacity arising from a process such as virga: when condensed material (rain) falls through the atmosphere before vaporizing again.

Cloud Properties
Interestingly, we find a slight positive correlation between the retrieval-derived radius of J1416A and the α parameter. With a more negative α, the cloud has a lower optical depth at longer wavelengths, allowing for flux to escape from hotter, brighter layers. The retrieval compensates for this to provide a good fit by reducing the scale factor (R 2 /D 2 ), resulting in a smaller radius estimate.

Retrieved Spectrum and Composition
Figure 6(a) compares the observed SpeX prism data and Sonora model spectra, which are cloud-free and consistent with the retrieved PT profile (see Figure 2). Figure 6(b) compares the retrieved forward model spectrum for the deck cloud model to the observed SpeX prism data. To compare the Sonora spectra to our retrieved forward model spectrum, the Sonora models were scaled to the median retrieved scale factor. Even though J1416A is best fit with a power-law deck cloud, the fits to the cloudless Sonora models are not very far off. This is  likely due to the deck cloud affecting only a small portion of J1416A's spectrum; thus these models can do a fair job at fitting the observed data. When comparing the observed spectrum of J1416A to the Sonora model spectra, we find the 1900 K solar-metallicity model provides the best fit overall but struggles to fit features in the J band and the H-band plateau. The J band is best fit by the 1800 K solar model, while the peak of the H band is best fit by the 1900 K low-metallicity model, and while the 1900 K solar model fits some of the K-band pseudo continuum, it is a poor match to the CO feature.
In Figure 6(b), the retrieval spectrum fits the overall shape of the observed spectrum quite well but has difficulties fitting the Na I doublet, K I doublets, and the FeH feature between the K I doublets in the J band. Issues in fitting the Na I and K I doublets are likely due to how the pressure broadening is treated in the opacity models for these lines. With pressure broadening from the 0.77 μm K I doublet impacting the slope in the J band through about 1.1 μm, the retrieved spectrum is likely unable to fit both the broad slope of the J band in this region as well as the narrow K I and Na I doublet features. We find that the Allard alkali opacities provide a better fit to J1416A than the Burrows alkali opacities, discussed in further detail in Section 8.1. In the H band, the retrieval does a much better job of fitting the FeH band to the data. This is likely driven by the H-band feature being broader than the J-band FeH feature and thus has a larger impact on the goodness of fit. The FeH fitting issue is an example of a problem introduced by the assumption of uniform-with-altitude mixing ratios, as the Jand H-band features are at different pressures and should have different abundances at those pressure layers.
Figure 6(c) shows the retrieved abundances for constrained gases compared to the solar-metallicity and solar C/O thermochemical equilibrium model values from the grid introduced in Section 4.3. Here we see the Na+K and H 2 O abundances are less than expected from models, pointing toward a subsolar metallicity for J1416A. The median retrieved CO abundance is also less than the solar model value but is just Figure 6. (a) Retrieved forward model spectra for the deck cloud model of J1416A. The maximum-likelihood spectrum is shown in dark green, the median spectrum in yellow, and 500 random draws from the final 2000 samples of the EMCEE chain in red. The SpeX prism data are shown in black. For comparison the cloud-free Sonora grid model solar-metallicity spectra for log g=5.0 and T eff =1600 K, 1700 K, and 1800 K (solid teal, blue, and purple), as well as [M/H]=−0.5 for log g=5.0 and T eff =1800 K and 1900 K (dotted teal, blue, and purple), are shown. These T eff values bracket the range of the SED-derived and retrieval-derived T eff . (b) Retrieved uniform-with-altitude mixing abundances for constrained gases compared to solar-metallicity and C/O model abundances. The approximate location of the photosphere is shown in gray.
within the 1σ confidence interval. The photosphere is shown as a gray strip to guide reasonable abundance ranges for metal oxides and metal hydrides. As these are not close to uniformwith-altitude, it is difficult to compare our retrieved values to the models. We do find that our abundances for TiO, VO, CrH, and FeH all fall within the very wide range of possible model abundances in the photosphere. Examining our FeH abundance, we see that the retrieved value is less than the maximum abundance of ≈−6 that is possible in the photosphere. This maximum abundance corresponds to those deeper into the atmosphere where we see the J-band FeH feature. With our lower than expected abundance, this points toward Fe being condensed in the atmosphere and agrees with Fe as our predicted cloud species.
Interestingly, we find that the uniform-with-altitude model is preferred over the thermochemical equilibrium model. At these temperatures, J1416A is expected to be in thermochemical equilibrium as the thermochemical timescale should be faster than the mixing timescale (Visscher et al. 2006, Section 5.1). Based on Figure 6(b), the alkalies are likely to be driving this preference as their abundance is the only one that is discrepant with the thermochemical grid abundance. Therefore, the uniform-with-altitude method is able to capture this discrepancy while still allowing for the other gas abundances to be in agreement with the thermochemical grid.

It's a Different Cloud, Which Is Indistinguishable: The
Power-law Slab Cloud Tells the Same Story As listed in Table 4, the power-law slab cloud is indistinguishable from the power-law deck cloud model and thus should tell a similar story about the atmosphere of J1416A. Here we present the retrieval results of the power-law slab cloud retrieval.

PT Profile and Contribution Function
Figure 7(a) shows the retrieved PT profile, slab cloud location, and total optical depth of the cloud. For this model, we find the bulk of the flux roughly between 1 and 18 bars like the deck cloud. The median retrieved profile in this region agrees within the 1σ confidence interval with the Sonora solar metallicity, log g=5.0, 1900 K, and 1700 K models and the [M/H]=−0.5, log g=5.0, 1900 K model. Compared to the 1700 K/5.0/solar and 1900 K/5.0/−0.5 models, the retrieved profile is slightly hotter at the same pressure, while it is slightly cooler than the 1900 K/5.0/solar model at the same pressure. At higher pressures, deeper in the atmosphere, the retrieved profile has a similar slope to that of the Sonora models, while at pressures lower than the photosphere the retrieved profile is more isothermal than the models. This is similar to the behavior of the power-law deck cloud profile compared to the Sonora models. The location in pressure space of the slab cloud, as well as its vertical height, are both poorly constrained due to the cloud being primarily optically thin with a total median optical depth across the cloud thickness of τ=1.08 at 1 μm, with a λ −1.27 drop-off to longer wavelengths. Figure 7(b) shows the contribution function for this model, which shows the opacity from the slab cloud having a small effect on the overall flux emitted. The optically thick portion of the slab cloud is only between ∼1 and 1.06 μm, whereafter it becomes optically thin and no longer significantly contributes to the observed flux. In the optically thick Y-band region we see that even though the cloud contributes ∼1% of the total flux in this region, the optical depth of τ median =1.08 is primarily from here. Unlike the deck cloud, we see that the slab contributes to the flux at higher altitudes; however, this only contributes ∼1%-10% of the total flux observed. The lack of cloud opacity in the J band contributes to the unusually blue J−K color in the same way as the power-law deck cloud model. With the cloud only affecting part of the Y band, the flux from the J band likely coming from a deeper pressure layer than that of field L dwarfs causing the bluer J−K color (see the comparison in Section 5.1.1).  Table 6 for ease of readability. Comparisons to the chemical equilibrium grid values of the gases are discussed in Section 5.2.4. The The contribution function associated with this cloud model, with the median cloud (magenta) and gas (aqua) at an optical depth of τ=1 overplotted. majority of the gas abundances, T eff , radius, mass, C/O, and [M/H] values for the slab model agree with those from the deck cloud model. The exception is the Na+K abundance, which differs from the deck cloud abundance by 1.4σ. This key difference in alkali abundance will be discussed in more detail in Section 8.1, when we compare the alkali abundances between the retrievals for J1416A and J1416B.

Cloud Properties
Retrieved cloud properties for the total optical depth, the pressure level for the base of the cloud (log P base ), the height of the cloud, the single scattering albedo, and the wavelength exponent α that describes how "nongray" the cloud is for the slab model are shown in Figure 9. The cloud base, height, and albedo are unconstrained for this model. The power α is more Figure 8. J1416A power-law slab cloud posterior probability distributions for the retrieved parameters and extrapolated parameters. One-dimensional histograms of the marginalized posteriors are shown along the diagonals, with 2D histograms showing the correlations between the parameters. The dashed lines in the 1D histograms represent the 16th, 50th, and 84th percentiles, with the 68% confidence interval as the width between the 16th and 84th percentiles. Parameter values listed above are shown as the median± 1σ. Gas abundances are displayed as log 10 (X) values, where X is the gas. T eff , radius, mass, C/O ratio, and [M/H] are not directly retrieved parameters but are calculated using the retrieved R 2 /D 2 and log g values along with the predicted spectrum. Our derived C/O ratio is absolute, where solar C/O is 0.55, while our [M/H] is relative to solar. Values for CO 2 and TiO are not constrained and thus only provide upper limits.
tightly constrained than in the deck cloud model and agrees within 1σ. The slab cloud also has a negative power, corresponding to a reddening cloud with submicron-sized particles likely described by a Hansen distribution. As the slab cloud is higher in the atmosphere, multiple condensates are stable at its location. As the slab and deck cloud models are indistinguishable, distinguishing between the condensates is critical to atmospheric understanding and will be the subject of future work. Like the deck cloud, the slab cloud also has a positive correlation between the radius and α, causing a smaller opacity at longer wavelengths, thus allowing for a smaller radius.

Retrieved Spectrum and Composition
The forward model spectrum for the slab cloud model is shown in Figure 10(a) compared to the observed SpeX spectrum and various temperature and metallicity Sonora models that bracket the retrieved T eff . For the slab cloud forward maximum-likelihood model spectrum, we find it is best fit by the 1700 K solar-metallicity model in the J band, while the 1800 K solar metallicity or 1800 K [M/H]=0.5 models fit better in the H and K bands. In Figure 10(b), we compare the retrieved spectrum and the observed spectrum. The spectrum from the slab cloud model is quite similar to that of the deck cloud, fitting both the FeH feature and the 1.25 μm K I doublet in the J band poorly, for similar reasons as discussed in Section 5.1.4. Figure 10(c) compares the retrieved gas abundances for the constrained gases to the solarmetallicity values expected from the thermochemical equilibrium model values from the grid introduced in Section 4.3. Unlike the deck cloud model, the retrieved CO abundance is below the solar model expected values. All of our retrieved gas fractions for this model are consistent with the deck cloud model, with the exception of the Na+K abundance. These low abundances of H 2 O, CO, and the tied Na+K again confirm the low-metallicity atmosphere that we derive.

Retrieved Model of 1416B
We initially used the Burrows alkali opacities as done in Line et al. (2017) for J1416B, which produced the best-fit model. However, we find that the Allard alkali opacities give consistent abundances between J1416A and J1416B, and thus we effectively treat the cloud-free Allard alkali model as the best model for J1416B. Thus in this section, we present the results of the second best fitting model (our winning model, ΔBIC=10), the cloud-free, uniform-with-altitude mixing ratio, Allard alkali opacity model for J1416B. The ΔBIC for all tested models for J1416B are listed in Table 7, and the cloud-free Burrows alkali opacity model results are shown in Section A.2. Detailed examination of our choice of alkali line models is discussed in Section 8.1. We will compare our J1416B results to retrieval results from Line et al. (2017) throughout this section. Intercomparison of retrieved and extrapolated parameters between J1416B and J1416A, as well as comparisons to the literature, will be discussed in Section 8. Figure 11(a) compares our retrieved median profile to Sonora 500 K, 600 K, and 700 K solar and [M/H]=−0.5 models and the median retrieved profile from Line et al. (2017). We see that our retrieved profile has a similar slope and is consistent within 1σ across the entire profile with all models except the solar 700 K and low-metallicity 500 K Sonora models. Compared to the median profile from Line et al. (2017), we find our profile consistent within 1σ; however, the shape of our profile differs from Line et al. (2017) at pressures  Notes.Molecular abundances are fractions listed as log values. For unconstrained gases, 1σ confidence is used to determine the upper limit. a Additional comparatives are listed in Table 1. b Atmospheric values. below ∼−0.5 bar. Many of the retrieved T dwarf profiles in Line et al. (2017) were more isothermal than the models, and they suggested it could be due to additional heating; however, temperature constraints are unreliable in this region of the profile. Figure 11(b) shows the contribution function for this model with the photosphere ranging from about 1 to 100 bars.

Retrieved Gas Abundances and Derived Properties
Posterior probability distributions for gases, surface gravity, T eff , radius, mass, C/O, and [M/H] are shown in Figure 12 with their values along with the derived L bol listed in Table 8. 29 ), while our radius, surface gravity, and metallicity agree within 1σ. Comparing our retrieved gas abundances to those of Line et al. (2017), we find all the gases we have in common are consistent except for the Na+K alkali abundance. Line et al. (2017) used the Burrows alkali opacities, while we use the Allard opacities in this model. When comparing with our model that used the Burrows opacities we find the Na+K abundance is consistent with Line et al. (2017). Similar to Line et al. (2017), we detect ammonia with our constraints equally as tight.
Our retrieved abundances yield a C/O ratio of = C O -+ 0.52 0.07 0.09 . To consider the effect of oxygen sequestration by Figure 10. (a) Retrieved forward model spectra for the slab cloud model of J1416A. The maximum-likelihood spectrum is shown in dark green, the median spectrum in yellow, and 500 random draws from the final 2000 samples of the EMCEE chain in red. The SpeX prism data are shown in black. For comparison, the Sonora grid model solar-metallicity spectra for log g=5.0 and T eff =1600 K, 1700 K, and 1800 K (solid teal, blue, and purple), as well as [M/H]=−0.5 for log g=5.0 and T eff =1800 K and 1900 K (dotted teal, blue, and purple), are shown. These T eff values bracket the range of the SED-derived and retrieval-derived T eff . (b) Retrieved uniform-with-altitude mixing abundances for constrained gases compared to solar-metallicity and C/O model abundances. The approximate location of the photosphere is shown in gray.  Line et al. (2017) value, and is subsolar relative to the solar C/O ratio of C/O=0.55 (Asplund et al. 2009). It should be noted that the correction used from Line et al. (2017) (where 3.28 oxygen atoms are removed per silicon atom) is under the assumption of uniform metallicity variations in elemental abundance ratios (e.g., Si/H∼M/H; see Visscher et al. 2010), as variations in the abundances of rock-forming elements (such as Mg and Si) will affect the proportion of oxygen removed by silicate condensation. However, as J1416B is subsolar, corrections to the C/O ratio may differ as subdwarf atmospheres have weak or absent metal oxides. If there is a relative depletion or lack of rock-forming elements, less oxygen would be sequestered, yielding a smaller correction in the C/O ratio. In Figure 13 of Nissen et al. (2014), the authors show that as metallicity ([Fe/ H]) decreases, the C/O is expected to decrease for thin-disk stars. Using our uncorrected metallicity, we find that our C/O ratio lies within the scatter of their expected metallicity prediction. The Line et al.  Figure 13(a) compares our retrieved median and maximumlikelihood spectra to the SpeX prism J1416B data and the bestfitting Sonora solar and [M/H]=−0.5 grid model spectra. We find that our retrieval spectrum fits quite well, with the exception of the Y-band peak being slightly below the data. In comparison to the Sonora model spectra, we find none of the models fit the Y-band peak, the 600 K solar-metallicity model does a good job fitting the J-band peak but is unable to fit the slope on either side quite well, and the H-and K-band features are best fit by the 700 K low-metallicity model. We find our retrieved gas abundances for H 2 O, CH 4 , and NH 3 are subsolar in the photosphere, while the alkalies are broadly consistent with the solar value in Figure 13(b). These values are consistent with those in Line et al. (2017). We find that relatively small changes in composition can drive major observable differences in the spectrum, particularly at lower temperatures. Therefore, with a slightly subsolar metallicity for J1416B, its spectrum differs quite drastically from field T dwarfs. Table 9 compares our SED-and retrieval-based fundamental parameters to the literature. Additionally, we list new UVW values using Gaia Collaboration et al. (2018) proper motions and parallax along with the radial velocity from Schmidt et al. (2010). Our empirical L bol is 2.5σ and 1.5σ discrepant from our deck and slab retrieval-based bolometric luminosities, respectively; however, all three measurements have very small uncertainties. The largest discrepancy between our SED and retrieval-derived parameters are the T eff and radius, with our T eff for the deck cloud at minimum 81 K hotter and the slab 50 K hotter than the semiempirical T eff of 1694 K. This is due to our small retrieved radius of R deck =0.7±0.04, = R slab -+ 0.77 0.06 0.10 , which is about 20% smaller than the evolutionary model radius from the SED method (see Section 7.3 for further discussion). Compared to the literature, our retrieval-based T eff is hotter than all, except the model-based T eff from Bowler et al. (2010) (which also calculates a L bol , but using an atmospheric spectra model), while the retrieval-based masses are consistent with our SED method value and Bowler et al. (2010). The log g we derive agrees between the SED and the retrieval methods. As this work is the first to derive a metallicity for J1416A, we find that the metallicity is consistent between both cloud models.

J1416B Fundamental Parameter Comparison
Table 10 lists our SED and retrieval method fundamental parameters compared to the literature. Comparing L bol , we find that both our SED and retrieval method values agree with Filippazzo et al. (2015) within 1σ. Our semiempirical and retrieval-based T eff radius, mass, and log g are consistent with one another and the literature within 1σ, with the exception of As the SED-based parameters T eff , mass, and radius are drawn from different evolutionary models (see Section 3 for evolutionary models used), these are plotted for comparison to the retrieval-inferred values and not to the evolutionary models themselves.
A comparison of radius versus L bol is shown in Figure 14(a), with the retrieval shown in black and the SED in pink. It is Figure 12. J1416B cloud-free posterior probability distributions for the retrieved and derived parameters using the Allard alkalis. One-dimensional histograms of the marginalized posteriors are shown along the diagonals, with 2D histograms showing the correlations between the parameters. The dashed lines in the 1D histograms represent the 16th, 50th, and 84th percentiles, with the 68% confidence interval as the width between the 16th and 84th percentiles. Parameter values listed above are shown as the median± 1σ. Gas abundances are displayed as log 10 (X) values, where X is the gas. T eff , radius, mass, C/O ratio, and [M/H] are not directly retrieved parameters but are calculated using the retrieved R 2 /D 2 and log g values along with the predicted spectrum. Our derived C/O ratio is absolute, where solar C/O is 0.55, while our [M/H] is relative to solar. T eff , radius, mass, C/O ratio, and [M/H] are not directly retrieved parameters but are calculated using the retrieved R 2 /D 2 and log g values along with the predicted spectrum. CO abundance is not constrained and thus only provides an upper limit. quite clear that the derived retrieval radius for the deck cloud model of J1416A is smaller than predicted by the evolutionary models, while the slab cloud is consistent with the lowmetallicity 6-10 Gyr and the solar-metallicity 3-6 Gyr models. While the radius of the deck cloud model may appear to be unphysically small, Sorahana et al. (2013) estimated the radii of brown dwarfs from the scale factor, similar to our method, using AKARI spectra and found that most of their mid-to late-L dwarfs had radii smaller than predicted from evolutionary models. We find both the deck and slab cloud radii for J1416A fall within the radius range of 0.64-0.81 R Jup for the mid-to late-L dwarfs in Sorahana et al. (2013) with a T eff between 1500 and 2000 K. The problem of unphysically small radii is an ongoing problem for atmospheric retrievals (e.g., Zalesky et al. 2019) and has been seen as an issue for the directly imaged exoplanets as well. We caution the reader against using the retrieved radii for J1416A.
Compared to our SED method radius, we see that it is only consistent with J1416A's slab cloud model radius. For J1416B, the retrieval radius is consistent with the Sonora evolutionary models and the SED method radius. As seen with J1416A, J1416B's SED method radius is larger than the retrievalderived radius. J1416A and J1416B's empirical L bol from the SED are fainter than the retrieval-derived L bol , which is inferred from integration under the retrieved forward model spectrum. The retrieval-derived radius for J1416B constrains the age to be >6 Gyr.
Comparison of the retrieved and evolutionary model-based (from the SED method) surface gravity versus L bol compared to the Sonora Bobcat evolutionary models is shown in Figure 14(b). The surface gravity for both retrieval models of J1416A is consistent with the SED value, and the same is seen between J1416B's retrieved and SED surface gravity. Here we see that both J1416A's slab and deck retrieval, as well as the SED, log g gives an age range of 1-10 Gyr. For J1416B, we find the retrieved log g produces an age range of 1-10 Gyr, which is broader than the range given from the radius. Figure 14(c) compares the log g versus T eff , where here we also compare J1416A to literature results from model values in Cushing et al. (2010) and the retrieval results for J1416B from Line et al. (2017). While the log g for J1416A and J1416B are consistent across the SED, retrievals, and the literature values plotted, the T eff measurements vary over a wider range, particularly for J1416A. When comparing mass versus L bol in Figure 14(d), we find that the retrieval places J1416A with a very young age of likely less than 1 Gyr, which is strikingly different from the very old age estimate from the radius. This age disagreement is likely due to the mass being tied to the radius, and with a larger radius the derived mass range would be higher.   Table 11 lists all the parameters we will discuss when comparing between J1416A and J1416B, particularly those of interest when determining whether the system formed and evolved together. We list the retrieved alkali abundances, C/O, and [M/H] determined when using both the Allard and  . d For the retrieval, the distance-calibrated spectrum from the SED was used; thus it was set to a distance of 10 pc. Distance uncertainty is included for determining the extrapolated parameters using the measured distance uncertainty. e An estimated distance of 9.4±1.3 pc is given assuming a low metallicity and using the Cushing et al. (2009)    h For the retrieval, the distance-calibrated spectrum from the SED was used; thus it was set to a distance of 10 pc. Distance uncertainty is included for determining the extrapolated parameters using the measured distance uncertainty. i Parallax from Faherty et al. (2012) was used. j Parallax from Dupuy & Liu (2012) was used.

Discussion
Burrows alkali opacities for both J1416A and J1416B. Here we use the C/O and [M/H] ratios determined from using only the gases that both J1416A and J1416B have in common (H 2 O, CH 4 , and CO). Agreement in the expected behavior of the alkali abundances was the primary deciding factor on the preferred cross sections.

Addressing the Differences in Alkalies
The alkali abundances retrieved for J1416AB are listed in Table 11, using both alkali opacity models. The Allard opacities are able to produce consistent alkali abundance between J1416A and J1416B, only when J1416A is parameterized with the deck cloud. Alkali abundances do not necessarily need to be consistent between J1416A and J1416B because they are condensing out at around the T eff of J1416B Zalesky et al. 2019). However, the Burrows opacities result in J1416AB having a higher alkali abundance than J1416A, which is not expected to occur in T dwarfs due to rainout. To check for correlations or degeneracies between alkali abundance and the cloud parameters of both cloud models for J1416A, we created a corner plot using the Allard  opacity retrieval results and found no correlations for either cloud model. Because the Allard alkalies produce the expected alkali abundance behavior between J1416A and J1416B with the deck cloud and not the slab cloud, this is evidence that the deck cloud produces a more realistic fit to the data over the slab cloud for J1416A.

C/O Ratio
To compare the C/O ratio between J1416A and J1416B, we have derived an atmospheric C/O ratio that only considers the gases in common between both sources (H 2 O, CO, and CH 4 ), due to the differing gas assumptions in the L and T dwarf retrievals. For the L dwarf, there will be a small contribution from VO missing in the oxygen total. However, as VO has a very small abundance it does not make a large impact on the overall C/O ratio. Using this C/O ratio, we find that J1416AB are consistent within 1σ, which points toward evidence in favor of formation and evolution as a pair. Both J1416A and J1416B are approximately solar in C/O and have slightly subsolar metallicities. Considering the various methods to determine the C/O for J1416AB, all methods are consistent within 1σand do not differ based on which alkali opacities are used.
As a note, the C/O ratios in Table 11 do not include the rainout correction, as we have not made any corrections to the C/O ratio of J1416A. The rainout correction applied to J1416B accounts for oxygen that should be in the atmosphere above any deep cloud not detected in the retrieval. For J1416A, because the retrieval prefers a cloudy model, we have an entirely different situation to consider. If a correction is necessary for J1416A it would likely be a smaller amount, because a much smaller fraction of the total atmosphere is above the cloud (i.e., for the median slab or deck cloud we would be accounting for oxygen above ∼0.1 bar) than is the case for J1416B. In addition, we should consider oxygen tied up in SiO gas in J1416A. Considering this, the correction for J1416A could range from 0.5% to 12%, which is well within our 68% confidence interval of our C/O ratio.

Metallicity Differences?
To compare the metallicity between J1416A and J1416B, it is important to remember that the gases used to derive the individual atmospheric metallicities differed between the L and T dwarf atmospheres. To account for this, we take the same approach as for the C/O ratio and determine a metallicity using only the gases in common between the L and T dwarf to determine the elemental abundances in our metallicity calculation. This approach does not include elements that are expected to have a large portion taken up by unobservable sinks such as N 2 or condensation of iron in the L dwarf. However, both nitrogen (for J1416B) and iron (for J1416A) would affect the metallicity determination at the 10% level, well within our 68% confidence intervals.
The 1σconfidence intervals are quite large for both alkali opacity variant metallicities, with the Allard opacities producing consistent [M/H] between J1416A and J1416B regardless of J1416A's cloud model. When using the Burrows opacities, the derived metallicities are inconsistent between J1416A's deck cloud model and J1416B. It should be noted that both alkali models produce a lower median metallicity for the deck cloud compared to the slab; however, only the Allard model is consistent. Additionally, the Burrows model produces a higher median metallicity for J1416A, but a lower median metallicity for J1416B. Only the Allard opacities produce a consistent picture of the comoving pair.

Conclusions
In this work we present the first distance-calibrated SED of J1416A and an updated distance-calibrated SED of J1416B. We present the first retrieval of J1416A and the second retrieval of J1416B. J1416A is best parameterized by a power-law deck cloud model; however, it is indistinguishable from a power-law slab cloud model, while J1416B is best fit by a cloud-free model, agreeing with previous results from Line et al. (2017). For both cloud models of J1416A, we find our retrieval radius is smaller than the evolutionary model radius and inconsistent within 1σ. We also find that the retrieval produces a hotter T eff than the SED to compensate for the smaller radius and to maintain the same flux we observe. We find that relatively small changes in the composition can drive major changes in observed features in the spectrum, particularly for low temperature sources.
Examining the retrieval results across the pair, we find that only the Allard alkali opacities produce alkali abundances expected for J1416AB (with the T dwarf abundance lower than that of the L dwarf) and only for the deck cloud model for J1416A. Both J1416A and J1416B have slightly subsolar metallicities that are consistent with each other, no matter the chosen alkali opacity model. J1416AB is consistent with an approximately solar C/O ratio, with the median value slightly supersolar for J1416A and slightly subsolar for J1416B. These results point toward the pair having formed and evolved together. Retrieval results of this binary are the first look from a larger sample that aims to dive deeper into understanding subdwarf atmospheres by asking (1) are subdwarfs cloudless? and (2) how do their PT profiles compare to similar spectral type or T eff sources (Gonzales et al. 2020, in preparation)? Having both cloudy and cloud-free results from this work provides a step in understanding the nuances of metallicity in L and T dwarf atmospheres. Figure 16. J1416A power-law deck cloud posterior probability distributions for the retrieved parameters and extrapolated parameters. One-dimensional histograms of the marginalized posteriors are shown along the diagonals, with 2D histograms showing the correlations between the parameters. The dashed lines in the 1D histograms represent the 16th, 50th, and 84th percentiles, with the 68% confidence interval as the width between the 16th and 84th percentiles. Parameter values listed above are shown as the median± 1σ. Gas abundances are displayed as log 10 (X) values, where X is the gas. T eff , radius, mass, C/O ratio, and [M/H] are not directly retrieved parameters but are calculated using the retrieved R 2 /D 2 and log g values along with the predicted spectrum. Our derived C/O ratio is absolute, where solar C/O is 0.55, while our [M/H] is relative to solar. CO 2 , CH 4 , and TiO abundances are not constrained and thus only provide upper limits. Figure 17. J1416A power-law deck cloud posterior probability distributions for the cloud parameters. The cloud top pressure (log P top ) and the cloud height (dP) are shown in bars, and α is from the optical depth equation τ=τ 0 λ α . Figure 18. (a) Retrieved forward model spectra for the deck cloud model of J1416A. The maximum-likelihood spectrum is shown in dark green, the median spectrum in yellow, and 500 random draws from the final 2000 samples of the EMCEE chain in red. The SpeX prism data are shown in black. For comparison the Sonora grid model solar-metallicity spectra for log g=5.0 and T eff =1600 K, 1700 K, and 1800 K (solid light green, teal, and blue), as well as [M/H]=−0.5 for log g=5.0 and T eff =1800 K and 1900 K (dotted blue and purple), are shown. These T eff values bracket the range of the SED-derived and retrieval-derived T eff . (b) Retrieved uniform-with-altitude mixing abundances for constrained gases compared to solar-metallicity and C/O model abundances. The approximate location of the photosphere is shown in gray.  The maximum-likelihood spectrum is shown in dark green, the median spectrum in yellow, and 500 random draws from the final 2000 samples of the EMCEE chain in red. The SpeX prism data are shown in black. For comparison the Sonora grid model solar-metallicity spectra for log g=5.0 and T eff =1600 K, 1700 K, and 1800 K (solid teal, blue, and purple), as well as [M/H]=−0.5 for log g=5.0 and T eff =1800 K and 1900 K (dotted teal, blue, and purple), are shown. These T eff values bracket the range of the SED-derived and retrieval-derived T eff . (b) Retrieved uniform-with-altitude mixing abundances for constrained gases compared to solar-metallicity and C/O model abundances. The approximate location of the photosphere is shown in gray.

A.2. J1416B Burrows Models
Figures 23-25 show the cloud-free Burrows alkali crosssection model for J1416B, which presents a better fit to the data; however, it produces inconsistent alkali abundances between J1416A and J1416B.  are not directly retrieved parameters but are calculated using the retrieved R 2 /D 2 and log g values along with the predicted spectrum. CO abundance is not constrained and thus only provides an upper limit.