We present an atmospheric retrieval analysis on a set of young cloudy red L dwarfsCWISER J124332.12600126.2 BD60 1417B and WISEP J004701.06680352.1 W0047using the Brewster retrieval framework. We also present the first elemental abundance measurements of the young Kdwarf K0 host star BD60 1417 using highresolution R  50000 spectra taken with the Potsdam Echelle Polarimetric and Spectroscopic Instrument on the Large Binocular Telescope. In the complex cloudy Ldwarf regime the emergence of condensate cloud species complicates retrieval analysis when only nearinfrared data are available. We find that for both L dwarfs in this work despite testing three different thermal profile parameterizations we are unable to constrain reliable abundance measurements and thus the carbontooxygen ratio. While we cannot conclude what the abundances are we can conclude that the data strongly favor a cloud model over a cloudless model. We note that the difficulty in retrieval constraints persists regardless of the signaltonoise ratio of the data examined SN  10 for CWISER BD60 1417B and 40 for WISEP W0047. The results presented in this work provide valuable lessons about retrieving young lowsurfacegravity cloudy L dwarfs. This work provides continued evidence of missing information in models and the crucial need for JWST to guide and inform retrieval analysis in this regime.
['Brown dwarfs', 'Atmospheric composition', 'Substellar companion stars', 'Fundamental parameters of stars', 'L dwarfs', 'Stellar abundances', '185', '2120', '1648', '555', '894', '1577', 'Astrophysics - Solar and Stellar Astrophysics', 'Astrophysics - Earth and Planetary Astrophysics']
Retrieving Young Cloudy L Dwarfs A Nearby Planetarymass Companion BD60 1417B and its Isolated Red Twin W0047
{'Retrieving Young Cloudy L-Dwarfs: A Nearby Planetary-Mass Companion BD+60 1417B and Its Isolated Red Twin W0047': 'Caprice L. Phillips , 1 Jacqueline K. Faherty , 2 Ben Burningham , 3 Johanna M. Vos , 4 Eileen Gonzales , 5, 6, ∗ Emily J. Griffith , 7, † Sherelyn Alejandro Merchan , 2, 8 Emily Calamari , 9, 10 Channon Visscher , 11, 12 Caroline V. Morley , 13 Niall Whiteford , 10 Josefine Gaarn, 14 Ilya Ilyin , 15 and  \n15  \n3  \nKlaus Strassmeier  \n1 Department of Astronomy, The Ohio State University, 100 W 18th Ave, Columbus, OH 43210 USA  \n2  \nDepartment of Astrophysics, American Museum of Natural History, New York, NY 10024, USA  \nCentre for Astrophysics Research, School of Physics, Astronomy and Mathematics, University of Hertfordshire, Hatfield AL10 9AB  \n4 School of Physics, Trinity College Dublin, The University of Dublin, Dublin 2, Ireland  \n- 5 Department of Physics and Astronomy, San Francisco State University, 1600 Holloway Ave., San Francisco, CA 94132, USA  \n6 Department of Astronomy and Carl Sagan Institute, Cornell University, 122 Sciences Drive, Ithaca, NY 14853, USA  \n- 7 Center for Astrophysics and Space Astronomy, Department of Astrophysical and Planetary Sciences, University of Colorado, 389 UCB, Boulder, CO 80309-0389\n- 8 Hunter College, City University of New York, 695 Park Avenue, New York, NY 10065, USA 9 The Graduate Center, City University of New York, New York, NY 10016, USA  \n10 Department of Astrophysics, American Museum of Natural History, New York, NY 10024, USA  \n11 Chemistry & Planetary Sciences, Dordt University, Sioux Center IA 51250  \n12 Center for Extrasolar Planetary Systems, Space Science Institute, Boulder, CO 80301  \n13 Department of Astronomy, The University of Texas at Austin, Austin, TX 78712, USA  \n14 3 Centre for Astrophysics Research, School of Physics, Astronomy and Mathematics, University of Hertfordshire, Hatfield AL10 9AB,  \nUK  \n15 Leibniz-Institut for Astrophysics Potsdam (AIP), An der Sternwarte 16, D-14482 Potsdam, Germany  \n(Received April 26, 2024; Revised June 14, 2024; Accepted June 28 2024)', 'ABSTRACT': 'We present an atmospheric retrieval analysis on a set of young, cloudy, red L-dwarfs - CWISER J124332.12+600126.2 and WISEP J004701.06+680352.1 - using the Brewster retrieval framework. We also present the first elemental abundance measurements of the young K-dwarf (K0) host star, BD+60 1417 using high resolution (R = 50,000) spectra taken with PEPSI/LBT. In the complex cloudy L-dwarf regime the emergence of condensate cloud species complicates retrieval analysis when only near-infrared data is available. We find that for both L dwarfs in this work, despite testing three different thermal profile parameterizations we are unable to constrain reliable abundance measurements and thus the C/O ratio. While we can not conclude what the abundances are, we can conclude that the data strongly favor a cloud model over a cloudless model. We note that the difficulty in retrieval constraints persists regardless of the signal to noise of the data examined (S/N ∼ 10 for CWISER J124332.12+600126.2 and 40 for WISEP J004701.06+680352.1). The results presented in this work provide valuable lessons about retrieving young, low-surface gravity, cloudy L-dwarfs. This work provides continued evidence of missing information in models and the crucial need for JWST to guide and inform retrieval analysis in this regime.  \nKeywords: Brown dwarfs (185); Atmospheric composition (2120); Companion stars (291); Fundamental (894) parameters of stars (555); L dwarfs (1679); Stellar abundances (1577)', '1. INTRODUCTION': "Brown dwarfs are a class of astronomical objects that are not massive enough (13 - 75 M Jup ) to stably fuse hydrogen in their cores (Kumar 1963). As such, they contract and cool as they age and progress through spectral types (Kirkpatrick 2005). As the bridge between lowmass stars and directly-imaged gas-giants, brown dwarfs are key to understanding planetary atmospheres and planetary formation mechanisms (Faherty et al. 2016; Liu et al. 2016). Brown dwarfs that share the same parameter space (e.g., wide orbits/long orbital periods, colors, mass, radius, temperature, and spectral morphology) with directly-imaged planets can be complementary to study (Burgasser 2010; Faherty et al. 2013a,b, 2021). Furthermore, young, low mass brown dwarfs are valuable laboratories for refining models of atmospheric chemistry and cloud formation. Studies have shown evidence for a correlation between clouds and youth (Faherty et al. 2013a, 2016; Vos et al. 2018). At younger ages (low surface gravities), the sedimentation and appearance of clouds are affected, further highlighting the value of studying young, low mass brown dwarfs (Marley et al. 2012; Su'arez & Metchev 2022, 2023). As analogs of young directly-imaged exoplanets they are ideal objects to test and verify models and atmospheric retrieval techniques. The direct detection and atmospheric studies of exoplanets remains a challenge due to the high contrast ratio with the host star (Marois et al. 2010; Chilcote et al. 2017; Rajan et al. 2017; Wang et al. 2020; Ruffio et al. 2021; Whiteford et al. 2023). Brown dwarfs can be studied in far greater detail than directlyimaged exoplanets, as they are easier to observe. Overall, brown dwarfs have been used to inform our observational knowledge of extrasolar atmospheres and advance theoretical understanding of cool worlds (Burrows et al. 2011; Tremblin et al. 2017; Marley et al. 2021).  \nAs far back as the discovery of the first methane bearing brown dwarf, Gliese 229B (Oppenheimer et al. 1995; Nakajima et al. 1995) the exact formation mechanism for brown dwarfs has been debated. While they are largely thought to form like stars via cloud fragmentation, there are alternative pathways which might lead to the lowest mass isolated sources (Bate et al. 2002; Whitworth & Goodwin 2005; Whitworth & Stamatellos 2006; Bonnell et al. 2008). Disk fragmentation and core accretion are two formation mechanisms that overlap with giant, widely separated exoplanets. Distinguishing between  \nformation pathways has been tackled on the theoretical side (e.g., Whitworth 2007; Luhman et al. 2007; Whitworth et al. 2010; Whitworth 2018) and recent spacebased and ground-based data allows us to address it on the observational side with carbon-to-oxygen (C/O), [C/H] and [O/H] ratios (e.g., Madhusudhan et al. 2016; Maldonado & Villaver 2017; Gonzales et al. 2020; Burningham et al. 2021; Gonzales et al. 2022; Calamari et al. 2022; Vos et al. 2023).  \nC/O ratios may help inform what pathway of formation a given object followed, shedding light on whether brown dwarfs form in a disk like planets in our solar system or in a fragment of a molecular cloud like binary stars ( Oberg et al. 2011). For example, at large orbital radii, we expect particularly strong differential evolution between gas and dust in the disk, which may translate to large observational differences in composition between substellar objects and their primary stars. Pertaining to planets, a stellar C/O ratio is expected for objects formed through gravitational instability, while super-stellar C/O indicates accretion of solids. The C/O ratio of the companion can be compared to the host star to provide context for formation and evolution history. With spectroscopic observations of a host star, we can measure the star's photospheric abundances - a chemical fingerprint of the gas the system was born out of and then have the base level for comparison to a companion (Konopacky et al. 2013; Wang et al. 2022; Xuan et al. 2023; Hoch et al. 2023).  \nExoplanetary atmospheric retrieval (or the spectral retrieval method) refers to the inference of atmospheric properties of an exoplanet given an observed spectrum (Madhusudhan 2019). Using this method, one can retrieve the atmospheric abundances and in turn derive bulk abundances which can be compared to their host star values and extract clues to their formation mechanisms (e.g., Wang et al. 2022; Wang et al. 2023; Xuan et al. 2022).  \nWhile a suite of retrieval frameworks have been developed in recent years (see review from MacDonald & Batalha 2023), for this work, we use the Brewster (Burningham et al. 2017) retrieval framework. Brewster is optimized for cloudy L-dwarfs and was initially designed to model and capture the complexity in clouds in extrasolar atmospheres. Brewster has been utilized and tested on brown dwarfs in the L - T spectral type regimes (e.g. Burningham et al. 2017; Burningham et al. 2021; Gonzales et al. 2020; Gonzales et al. 2021; Gonzales et al. 2022; Calamari et al. 2022; Vos et al. 2023; Gaarn et al. 2023).  \nIn this paper, we employ the Brewster spectral retrieval framework to understand the atmospheric structure and composition of two young red L-dwarfs. In Section 2, we present literature data on both systems. In Section 3, we present spectroscopic data for CWISER J124332.12+600126.2 (hereinafter BD+60 1417B) and WISEP J004701.06+680352.1 (hereinafter W0047) along with newly obtained data for the host star BD+60 1417. In Section 4, we discuss SED analysis used in the Brewster framework. We discuss the host star parameters and abundance analysis in Section 5. In Section 6, we discuss the retrieval framework and modification settings used in our work. We present the retrieval analysis results for BD+60 1417B and W0047 in Sections 8 and 9 respectively. We provide key discussions in Section 10, and conclude in Section 11.", '2. YOUNG RED L-DWARF SPECTROSCOPIC TWINS: BD+60 1417B AND W0047': 'The nearby ( ∼ 45 pc) recently discovered L6 - L8 γ companion, BD+60 1417B, is a planetary-mass (15 ± 5 M Jup ) object prime for atmospheric characterization (Faherty et al. 2021). BD+60 1417B lies at a wide separation from its young (50 - 150 Myr) K0 host star ( ∼ 1662 AU). Faherty et al. (2021) found that this object lies in sparsely populated region of mass ratio vs separation parameter space, with only 8 other objects, with an estimated mass < 20 M Jup and an orbital separation > 1000 AU. BD+60 1417B has a mass ratio in line with known exoplanets but appears too widely separated to have formed as a planet in situ . Thus there are open questions regarding its formation history and mechanism.  \nW0047 (Gizis et al. 2015) is an isolated red L7 γ dwarf member of the ∼ 120Myr AB Doradus moving group. Faherty et al. (2021) found that it was a spectroscopic twin to BD+60 1417B. BD+60 1417B and W0047 have similar masses, temperatures, and extremely red colors of (J-K s ) = 2.72 ± 0.34 and (J-K) = 2.55 ± 0.2 respectively (see Table 1). Additionally, both objects show signatures of condensate cloud driven variability (e.g. Vos et al. 2018, Vos et al. in prep). Spitzer photometric monitoring was obtained for W0047, which showed high amplitude variability of 1.07 ± 0.04%. This variability measurement put W0047 amongst the highest Spitzer variability amplitudes detected (Vos et al. 2018). These results are consistent with near-IR variability amplitudes found for W0047 using HST (11%) (Lew et al. 2016). The similarities between W0047 and BD+60 1417B therefore drives our intercomparison of their resultant retrieved properties.', '3.1. Spectroscopic Data for BD+60 1417': "BD+60 1417 is a nearby ( ∼ 45 pc; Gaia Collaboration et al. 2023) solar type star (K0 spectral type). As summarized in Faherty et al. (2021), BD+60 1417 has numerous youth indicators including: a Li I abundance, measurement of a fast rotation rate, a position on color-magnitude and color-color diagrams indicating youth all of which point to an estimated age range of 50 - 150 Myr. However, until this paper, BD+60 1417 lacked a high resolution spectrum used to calculate solar abundances. Therefore, we obtained high-resolution optical spectra of BD+60 1417 (V = 9.364 ± 0.11) from the Potsdam Echelle Polarimetric and Spectroscopic Instrument (PEPSI; Strassmeier et al. 2015) on the Large Binocular Telescope (LBT) on 9th February 2023. Observations were taken with the 300 µ m fiber for a resolution R = λ/ ∆ λ = 50000 and with cross-dispersers (CD) III and VI to obtain high signal to noise spectra in the wavelength ranges 4800 - 5411 and 7419 - 9140 ˚ A. This wavelength regime and resolution contains multiple stellar absorption features for C, O, Mg, Si, Ca, and Al. BD+60 1417 was observed for 1260 s, reaching a signalto-noise ratio (S/N) of 795 in CD III and and 517 in CD VI. We plot selected regions of the spectra in Figure 1.  \nWe employ the spectroscopic data systems for PEPSI pipeline (SDS4PEPSI) to reduce the spectra of BD+60 1417, as described in Strassmeier et al. (2018). In short, this pipeline applies bias and flat-field corrections and estimates the photon-noise contribution to the spectra. It removes contamination from scattered light as well as cosmic rays. The pipeline conducts flux extraction and wavelength calibrates the spectra against Th-Ar lines, as well as corrects for the echelle's blaze function, fringing, and vignetting effects. Finally, the pipeline corrects for the star's radial velocity shifts to bring the spectra to rest-frame wavelengths and applies a continuum normalization by fitting a 2D smoothing spline in the cross-dispersion and dispersion directions.", '3.2. Spectroscopic Data for BD+60 1417B and W0047': 'For BD+60 1417B, we use the published 1 - 2.5 µ m Infrared Telescope Facility (IRTF)/SpeX Prism spectrum from Faherty et al. (2021) with R ∼ 100 and S/N ∼ 9. The spectra for W0047 consist of a 1 - 2.5 µ m IRTF/SpeX spectrum from Gizis et al. (2012) with R ∼ 150 and S/N ∼ 37 and a 1.10 - 1.69 µ m HST/WFC3 spectrum from Manjavacas et al. (2019) with R ∼ 130 and S/N ∼ 370 at 1.25 µ m.  \n4. FUNDAMENTAL PARAMETERS FROM SPECTRAL ENERGY DISTRIBUTIONS FOR BD+60 1417B AND W0047  \n<!-- image -->  \n<!-- image -->  \nFigure 1. PEPSI spectra of BD+60 1417 in selected wavelength regions of 7687 -7695 ˚ A (left), 7770 -7778 ˚ A (center), and 8770 -8777 ˚ A (right). Line features used in the abundance analysis have been labeled.  \n<!-- image -->  \nIn order to acquire fundamental parameters, we constructed and analyzed the SEDs for BD+60 1417B and W0047. We compiled all relevant data including the parallax measurements of BD+60 1417B and W0047 along with spectra and photometry to acquire a distance-calibrated SED. WISE W1 and W2 photometry was used from the catWISE2020 reject table, as BD+60 1417B was rejected from the main table because it was near a diffraction spike from the primary. We follow Faherty et al. (2021) and find the photometry to be usable for construction and analysis of the SED. We used the open-source package SEDkit (Filippazzo 2020) presented in Filippazzo et al. (2015) to construct the spectral energy distribution (SED) and derive the fundamental parameters of these twin analogs. Integrating under the flux calibrated SED, we calculated the L bol for each target. Pairing the L bol values and ages of 100 ± 50 Myr (Faherty et al. 2021) and 120 ± 50 Myr (Gizis et al. 2015) for BD+60 1417B and W0047, respectively. Along with the Marley et al. (2021) evolutionary models, we semi-empirically calculate the T eff and infer the radius, mass and log g . The spectral energy distribution of BD+60 1417B and W0047 are shown in Figure 2 and fundamental parameters derived from the SEDs are shown in Table 1.  \nThe SED study and derived fundamental parameters for BD+60 1417B and W0047 emphasize the similar nature of these sub-stellar objects, as indicated by their T eff being within ∼ 50K, identical radii and mass measurements that extend into the nominal space between brown dwarfs and planets.', '5. SPECTRAL PARAMETER AND ABUNDANCE DETERMINATION FOR BD+60 1417': "5.1. Stellar Parameter Determination  \nWe used the PEPSI reduced spectrum to determine stellar parameters and stellar abundances with MOOG,  \nTable 1. Observed properties and derived fundamental parameters for BD+60 1417B and W0047  \na Epoch J2000, ICRS  \nReferences -F21: Faherty et al. (2021); TW: This Work; C03: Cutri et al. (2003); C12: Cutri & et al. (2012) G15: Gizis et al. (2015), M21: Marocco et al. (2021); L16: Liu et al. (2016); G21: Gaia Collaboration et al. (2021)  \na LTE radiative transfer code (Sneden 1973), implemented through the iSpec wrapper (Blanco-Cuaresma et al. 2014; Blanco-Cuaresma 2019). Stellar parameters T eff , log( g ), [M / H], [ α/ M], and v micro are fit with equiv-  \nFigure 2. SEDs of BD+60 1417B and W0047. The distance-calibrated SEDs of BD+60 1417B and W0047 are shown in orange and blue, respectively. The SpeX/Prism spectrum for BD+60 1417B is shown in orange, with a shaded region indicating the flux uncertainties. Similarly for W0047 the SpeX/Prism and HST/WFC3 spectrum are shown in blue. Photometric points for W0047 are shown in a light blue and in light orange for BD+60 1417B. Circles are used for 2MASS, squares for WISE, and a down triangle for Pan-STARRS. The wavelength coverage for the photometric filters are indicated with horizontal lines. Error bars on the photometric points are smaller than the point size.  \n<!-- image -->  \nlent width analysis. Through iSpec, we fit the Fe i and Fe ii lines in the Gaia-ESO version 6 line list (Gilmore et al. 2012; Heiter et al. 2021) with Gaussian or Voigt profiles and measure their equivalent widths. We remove the strongest and weakest lines by requiring that the reduced equivalent width of the lines lie between 4.3 and -6.0 and by requiring that the lower-energy state of the lines lie between 0.5 and 5.0, as recommended by Mucciarelli et al. (2013), Blanco-Cuaresma et al. (2014), and Blanco-Cuaresma (2019). Using MARCS model atmospheres (Gustafsson et al. 2008) and Asplund et al. (2009) solar abundances, we implement MOOG equivalent width parameter determination within iSpec, which iteratively fits models of T eff , log( g ), [M / H], [ α/ M], and microturbulence velocity ( v micro ), searching for a model that brings the Fe lines into ionization and excitation equilibrium. We initiate the parameter search at T eff = 4993 K, log( g ) = 4 . 55 dex, and [M / H] = 0 . 0 dex, based on the stellar parameters in the TESS input catalog (Stassun et al. 2018) adopted in Faherty et al. (2021). iSpec returns best-fit stellar parameters of T eff = 5382 ± 164 K, log( g ) = 4 . 88 ± 0 . 06 dex, [M / H] = 0 . 37 dex, [ α/ M] = 0 . 0 dex, and v micro = 1 . 026 km/s.  \nIn addition to the above stellar parameters, MOOG and iSpec require values of macroturbulence ( v macro ), rotational velocity ( v sin( i )), spectral resolution ( R ), and  \nlimb-darkening coefficient to conduct stellar abundance determination with spectral synthesis. We set the resolution to that of the instrument ( R = 50000) and the limb-darkening coefficient to 0.6, as recommended by Blanco-Cuaresma et al. (2014). Because the broadening terms v macro and v sin( i ) are degenerate at R = 50000, we set v macro = 0 km/s and find the best fit v sin( i ) through a spectral synthesis fit to Fe i and Fe ii line windows with all other stellar parameters fixed to the values listed above Buder et al. (2021). We find that a v sin( i ) = 5 . 40 km/s best minimizes the residual between the synthetic and observed spectra.  \nTo calculate the error on our stellar parameters, we run the equivalent width and spectral synthesis fits described above at the upper and lower bounds of T eff to find the correlated deviations in log( g ), [M / H], v micro , and v sin( i ). We adopt an error of 0.12 dex on log( g ), 0.04 dex on [M / H], 0 . 26 km/s on v micro and 0 . 35 km/s on the broadening term v sin( i ). We list all best-fit stellar parameters and their errors in Table 2.  \nWe note that the photometric stellar parameters adopted as the initialization are not consistent with our calculated spectroscopic parameters. While the photometric T eff and log( g ) are derived from observables and used to calculate physical parameters such as stellar radii and mass, their adoption inhibits achieving ionization and excitation equilibrium in the Fe i and Fe ii  \nabundances. We adopt the spectroscopic parameters in further analysis as they are optimized for stellar abundance determination, but stress that these values should not be used to derive the stellar mass or radius (e.g., Jofr'e et al. 2019).  \nTable 2. Stellar parameters and abundances for BD+60 1417  \nNote -Unless noted, all [X/H] abundances are 1D LTE, relative to solar abundances in Asplund et al. (2009). We provide LTE and NLTE [O/H] abundances. The reported C/O ratio is calculated from the NLTE [O/H] value.", '5.2. Stellar Abundance Determination': "We determine the stellar abundances of C, O, Mg, Si, Ca, Al, and Fe with MOOG spectral synthesis implemented through iSpec. We use a curated line list created from referencing line lists from Fulbright (2000), GALAH DR3 (Buder et al. 2021), and GaiaESO (Gilmore et al. 2012; Heiter et al. 2021), and the solar spectrum (Moore et al. 1966) comprised of 5 C I lines, 3 O I lines, 16 Mg I lines, 27 Si I lines, 6 Ca I lines, 4 Ca II lines, 6 Al I lines, 265 Fe I , and 10 Fe II lines. While the full spectrum was normalized in the Spectroscopic Data systems for PEPSI pipeline, we fit and normalize the continuum in the 8 ˚ A region around  \neach line window with a spline to ensure local continuum normalization.  \nFor each element, we create windows of 5 ˚ A around all relevant line features fit for the elemental abundance, with iSpec iteratively creating synthetic stellar spectra and finding the template with the abundance value that minimizes the residual between the synthetic and observed spectrum. We further fit each line of each element individually, finding consistent abundances between the lines of different ionization states for Si and Fe. We report the best fit 1D LTE [C/H], [O/H], [Mg/H], [Si/H], [Ca/H], [Al/H], and [Fe/H] values in Table 2. We repeat the stellar abundance fits for each element with the upper and lower limits of the stellar parameters to calculate the abundance errors. The errors range from 0.01 dex (for Mg) to 0.18 dex (for O) and are also reported in Table 2.  \nIn this paper, we calculate 1D LTE abundances for BD+60 1417. 3D and non-LTE (NLTE) effects are known to be a source of abundance uncertainty and can often cause abundance corrections > 0 . 1 dex (e.g. Jofr'e et al. 2019). While we cannot integrate NLTE grids into iSpec, we use NLTE correction lookup tables from the Max Planck Institute for Astronomy (MPIA) 1 to estimate the 1D NLTE corrections for O I (Bergemann et al. 2021), Mg I (Bergemann et al. 2017), Si I (Bergemann et al. 2013), Ca I , and Ca II (Mashonkina et al. 2007). Tabulated NLTE corrections for C and Al are not yet available and are not accounted for in our reported errors. We find NLTE corrections of < 0 . 01 dex for all available Mg and Si lines. While some Ca lines have significant NLTE corrections (-0.1 to 0.04 dex), the median value is < 0 . 01 dex and does not impact our reported [Ca/H] abundance. The O triplet lines used in this analysis have a mean NLTE correction of -0.02 dex, leading to an NLTE [O/H] abundance of -0.02 dex.  \nFrom the [X/H] abundances, we can calculate the C/O, Mg/Si, and Ca/Al ratios. Each ratio represents the ratio of the number of atoms of element X to the number of atoms of element Y where X / Y = N X /N Y = 10 log N X / 10 log N Y , and log N X = log 10 ( N X /N H )+12. We calculate log N X as [X / H] + log N X , ⊙ , using solar abundances log N X , ⊙ of 8.43 for C, 8.69 for O, 7.60 for Mg, 7.51 for Si, 6.34 for Ca, and 6.45 for Al (Asplund et al. 2009). We report the resulting ratios in Table 2.", '6. BREWSTER FRAMEWORK FOR RETRIEVING BD+60 1417B AND W0047': 'To extract comparative abundances for the brown dwarfs BD+60 1417B and W0047, we used the atmospheric retrieval code called Brewster (Burningham et al. 2017). Below we give a short description of the Brewster framework we utilized in this work however we encourage the reader to see a detailed description of Brewster in Burningham et al. (2017, 2021). Most notably, we differ from Burningham et al. (2017, 2021) with the use of a version update (Calamari et al. 2022) that utilizes nested sampling via PyMultiNest (Buchner et al. 2014) instead of emcee .', '6.1. The Forward Model': 'The Brewster forward model 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 Marley et al. (1996); Saumon & Marley (2008); Morley et al. (2012). We set up a 64 pressure layer (65 levels) atmosphere with geometric mean pressures from log P = -4 to 2 . 4 bars spaced in 0.1 dex intervals.', '6.2. Thermal Profiles': "There are a range of thermal profiles used in retrieval analyses ranging from a five-parameters joint exponential power law (Madhusudhan & Seager 2009; Burningham et al. 2017) to a 17-parameters free profile (Line et al. 2015). For this work, we aim to test a variety of thermal profiles for BD+60 1417B and W0047 to determine the reliability of retrieved abundances, pressure-temperature (P-T) profiles, bulk properties from ground-based near-infrared data, and the effect different thermal profile parameterizations have on retrieved abundances due to the notorious difficulty to model low gravity objects (see thesis work by Niall Whiteford and 2 GRAVITY Collaboration et al. 2020). A similar exploration was conducted by Gonzales et al. (2022) & Rowland et al. (2023) for L-dwarfs in a different regime.  \nIn this work we explore three different thermal profiles for our sources.  \n- ■ Case 1 - Five point spline: The initial thermal profile tested is a computationally simple five point parameterization in which we specify five temperature pressure points: the top ( T top ), bottom ( T bottom ) and middle of the atmosphere ( T middle ) and two midpoints between the top and middle ( T q1 ), and bottom and middle ( T q 3 ). These points are calculated in order, beginning with T bottom ,  \nwhich is selected in range between zero and the maximum temperature defined in our prior which is 4000 K in this work. Then T top is chosen between zero and T bottom , T middle chosen between T top and T bottom , and the remaining two midpoints chosen between T top and T middle T q1 , T middle and T bottom (T q 3 ). A uniform prior is assumed for each temperature within its respective range. This does not allow for temperature inversions but can result in 'irregular ' profiles.  \n- ■ Case 2 - Lavie profile: We introduce a second thermal profile parametrization developed and implemented by Lavie et al. (2017) in the HELIOS retrieval code and subsequently implemented into other retrieval frameworks (Whiteford et al. 2023). The Lavie profile is a less flexible profile with a simpler parameterization consisting of two free parameters, κ 0 and T int . The profile used in Lavie et al. (2017), originates from Heng et al. (2014), using a reduced form of Equation 126:  \nT 4 = T 4 int 4 · ( 8 3 +3˜ mκ 0 ) (1)  \nwhere T int is the internal temperature and κ 0 the constant component of the infrared opacity. ˜ m is column density determined via P 0 = ˜ m · g , where g is the surface gravity.  \n- ■ Case 3 - Molli'ere profile: Lastly, we explore a modified version of the thermal profile parameterization implemented in (Molli'ere et al. 2020; hereafter Molli'ere 'Hybrid' Profile). In the Molli'ere 'Hybrid' Profile, the atmosphere is split into three distinct regions: the photosphere (middle altitudes), high altitude, and troposphere (low altitudes). In the photosphere region, the temperature is set via an Eddington approximation:  \nT ( τ ) 4 = 3 4 T 4 0 ( 2 3 + τ ) (2)  \nwhere T 0 is a free parameter and τ = δP α . For τ , δ and α are free parameters as well.  \nThe high altitude region is similar to that of Molli'ere et al. (2020). In the troposphere, our Molli'ere 'Hybrid' Profile differs from the implementation of Molli'ere et al. (2020), in that a dry adiabatic atmosphere (for H 2 /He) is used instead of a moist adiabatic. For a more detailed description of the Molli'ere 'Hybrid' Profile we refer the reader to Molli'ere et al. (2020).", '6.3. Gas Opacities': 'We consider the following absorbing gases in our analysis: H 2 O, CO, CO 2 , CH 4 , CrH, FeH, Na, and K. From Kirkpatrick (2005), these gases are chosen as they have been previously identified as important absorbing species in the spectra of L6 - L8 spectral types.  \nOptical depths due to absorbing gases for each layer are calculated using opacities sampled at a resolving power R = 10 , 000 from the Freedman et al. (2008, 2014) collection with the updated opacites from Burningham et al. (2017). Line opacites in our pressure-temperature range are tabulated in 0.5 dex steps for pressure and for temperature in steps ranging from 20 K to 500 K as we move from 75 K to 4000 K, which is then linearly interpolated to our working pressure grid.  \nWe include 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), and include Rayleigh scattering due to H 2 , He and CH 4 but we neglect the remaining gases. Free-free continuum opacities are included for H -and H -2 as well as boundfree continuum opacity for H -(Bell & Berrington 1987; Bell 1980; John 1988).', '6.4. Determining Gas Abundances': "We use the vertically constant mixing ratios method to determine gas abundances. The vertically constant method is sometimes referred to as a 'free' retrieval method which enables us to directly retrieve individual gas abundances under a computationally simple assumption of a single abundance for each gas throughout the 1D modelled atmosphere. It is considered 'free' in the sense that there are no constraints imposed by thermochemical self-consistency. However, the simplicity of this method can not capture important variations with altitude that are expected for some species which can vary by several orders of magnitude in the photosphere (Rowland et al. 2023).", '6.5. Cloud Model': "We use the same simple cloud model utilized in previous works (e.g., Burningham et al. 2017, Gonzales et al. 2020, and Vos et al. 2023) where the cloud is parameterized as a 'deck' or 'slab'. Both cloud types are defined where opacity due to the cloud is distributed among layers in pressure space and the optical depth defined as either grey or as a power-law ( τ = τ 0 λ α , where τ 0 is the optical depth at 1 µ m). As done in Gonzales et al. (2020), the single scattering albedo is set to zero thereby assuming an absorbing cloud.  \nThe deck cloud is defined to always become optically thick at some pressure top, such that we only  \nsee the cloud top and the vertical extent of the cloud at pressures lower than P top . There are three parameters for the deck cloud: (1) a cloud top pressure P top , the point at which the cloud passes τ = 1 (looking down) for the first time ( P deck ), (2) the decay height ∆log P , over which the optical depth falls to lower pressures as dτ/dP ∝ exp(( P -P deck ) / Φ) where Φ = ( P top (10 ∆log P -1)) / (10 ∆log P ), and (3) the cloud particle single-scattering albedo ( α ). At pressures P > P top , the optical depth increases following the decay function until ∆ τ layer = 100. Deep below the cloud top there is essentially no atmospheric information as the deck cloud can rapidly become opaque with increasing pressure. Therefore, we stress that the pressure-temperature (P-T) profile below the deck cloud top extends the gradient at the cloud top pressure.  \nThe slab cloud differs from the deck as it is possible to see the bottom of the slab. The slab parameters include deck parameters (1) and (3). Because it is possible to see to the bottom of the slab, parameter (2) is now a physical extent in log-pressure (∆ log P ). We also include an additional parameter for determining the total optical depth of the slab at 1 µ m ( τ cloud ), bringing the total number of parameters to 4. The optical depth is distributed through the slab's extent as dτ/dP ∝ P (looking down), reaching its total value at the highest pressure (bottom) of the slab. The slab can have any optical depth in principle, however, we restrict our prior as 0 . 0 ≤ τ cloud ≤ 100 . 0.  \nIf the deck or slab cloud is non-grey, an additional parameter for the power ( α ) in the optical depth is included.", '7. RETRIEVAL MODEL': "The retrieval process consists of optimizing the parameters of the forward model such that the resultant spectrum provides the best match to the observed spectrum. As described by Burningham et al. (2017), we use a Bayesian framework to optimize the model fit to the data by varying the input parameters. Following Calamari et al. (2022), we use a modified sampler to the original version of Brewster , that uses PyMultinest (Buchner et al. 2014) instead of emcee (Foreman-Mackey et al. 2013).  \nBrewster applies Bayes' theorem to calculate the 'posterier probability', p ( x | y ), the probability of a set of parameters' ( x ) truth value given some data ( y ), in the following way:  \np ( x | y ) = L ( x | y ) p ( x ) p ( y ) (3)  \nwhere L ( x | y ) is the likelihood that quantifies how well the data match the model, p ( x ) is the prior probability  \non the parameter set and p ( y ) is the probability of the data marginalized over all parameter values, also known as the Bayesian evidence.  \nThe posterior probability space is explored using the PyMultiNest sampler (Buchner et al. 2014) which utilizes nested sampling to discover the set of parameters with the maximum likelihood given the data. The samples in an n-dimensional hypercube, or state-vector, are translated into parameter values via the 'prior-map'. The prior-map function is how prior probabilities are set for each parameter and are then transformed into appropriate parameter values to be used in the forward model. This algorithm is equipped to handle a parameter space that may contain multiple posterior modes and/or degeneracies in moderately high dimensions.", '7.1. Model Selection': 'A variety of models were tested in our retrievals of BD+60 1417B and W0047. These models differed in their cloud prescriptions and temperature-pressure profile parameterizations, which are explained in Section 6. In order to compare all of our retrievals, we calculate the logarithm of the Bayesian evidence (logEv), where the highest logEv is preferred. For each of the models tested, the ∆logEv is calculated by subtracting the logEv of the winning model of the Bayesian evidence.  \nWe use the following significance intervals from Kass & Raftery (1995) to distinguish between two models, with evidence against the lower logEv as:  \n- · 0 < ∆logEv < 0.5: no preference worth mentioning\n- · 0.5 < ∆logEv < 1: positive\n- · 1 < ∆logEv < 2: strong\n- · ∆logEv > 2: very strong  \nWe began by building from the least complex model (cloud-free) to the most complex (power law slab cloud model) for the three thermal profile parameterizations tested in this work for BD+60 1417B and W0047 (see Table 4).', '8. BD+60 1417B RETRIEVAL RESULTS': 'We show the list of models tested for BD+60 1417B with the number of parameters and ∆logEv for each model and thermal profile parameterization in Table 4.', '8.1. Best-fit Models': "The top-ranked model for BD+60 1417B was the power law deck cloud model when using either the simple five point spline or the Molli'ere 'Hybrid' thermal profile. When using the Lavie thermal profile, the grey deck cloud model was top-ranked and the power law deck cloud was a close third. With ∆logEv greater than 2 in all three cases, it is clear that the cloudless model is strongly rejected, indicating that the data strongly favor a cloud model to fit the spectroscopic features observed in BD+60 1417B.", '8.2. Thermal Profiles': "We show the retrieved thermal profiles and cloud pressures for the winning models compared to Sonora Diamondback self-consistent grid models (Morley et al. 2024) and condensation curves for BD+60 1417B in Figure 3. For the grid models shown, the base model has a T eff = 1300 K, log g = 4.5, and a sedimentary efficiency, f sed = 1-3. These values are selected from the SED-derived parameters for BD+60 1417B.  \nFor the five point spline profile (Figure 3a) there exists a slight thermal inversion near the top of the atmosphere, which we conclude is unphysical. This phenomenon is a byproduct of the five point 'irregular' thermal profile (Calamari et al. 2022). Additionally, we note that the retrieved thermal profile is not in agreement with the Sonora Diamondback grid models at the predicted temperature, especially within the photosphere (indicated by light grey shading) where our retrieved profile exhibits primarily warmer temperatures than the models predict. In the bottom half of the observable photosphere, our retrieved profile is nearly isothermal. Near the top of the visible photosphere region, the f sed =1model is near the same temperature as our retrieved profile. The location of the median cloud deck (1 σ ) is within the observable photosphere for the five point thermal parameterization. Generally we find that the best fit retrieval using the five point spline thermal profile for BD+60 1417B looks physically implausible (Figure 3a)  \nFigure 3b shows the retrieved thermal profile for the winning model for the Lavie thermal profile. The Lavie profile is a less flexible profile, which forces an almost isothermal profile above the photosphere. In the photosphere region, the thermal profile for BD+60 1417B is cooler than Sonora Diamondback grid models predictions. The cloud top is near the bottom of the photosphere.  \nIn Figure 3c we show the retrieved thermal profile and cloud pressure for the winning model for BD+60 1417B for the Molli'ere thermal parameteriza-  \nTable 3. Priors for BD+60 1417B and W0047 retrieval models  \ntion. In the photosphere region, the retrieved profile is nearly isothermal and the median cloud lies within the observable photosphere, similar to the five-point spline model. In the regions extrapolated below the observable photosphere, the f sed = 2 and f sed = 3 models most closely trace the extrapolated profile.", '8.3. Contribution Functions and Retrieved Spectra:': "The contribution function in a layer is defined as  \nC ( λ, P ) = B ( λ, T ( P )) ∫ P 2 P 1 dτ exp ∫ P 2 0 dτ (4)  \nwhere 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 contribution function indicates the location in the atmosphere where observable emission features are produced from and where spectral features become muted. The τ = 1 gas and cloud contribution for winning models are shown as well.  \nFigure 4a shows the contribution function for the five point spline thermal profile model for BD+60 1417B for the winning model of the power law cloud deck. The τ = 1 gas and cloud contributions are shown. The majority of the flux contribution to the observed spectrum comes from approximately 0.004 to 2 bar region, corresponding to the photosphere. The K band is shaped by the cloud opacity, which has a more significant contribution. In the J and H bands, the gas opacity contributes more, which is the wavelength regions corresponding to H 2 O  \nmolecular band features. There are windows in the J ( ≈ 1.2 to 1.4 µ m) and H ( ≈ 1.6 to 1.8 µ m) bands where the cloud opacity dominates.  \nFigure 4b, shows the contribution function for the Lavie thermal profile for BD+60 1417B for the winning model, the power law cloud deck. The bulk of the flux between ≈ 0.7 and 20 bars contributes to the observed spectrum. Across the spectrum, the cloud opacity becomes more important, with the exception of a small region in the H band ( ≈ 1.8 - 2.0 µ m), where the gas opacity dominates.  \nFigure 4c shows the contribution function for the Molli'ere thermal profile. The contribution function is similar to the five point spline thermal profile. One of the key differences is the dominance of the cloud opacity in the entire J band region. Between 1.3 - 1.5 µ m is the only region where the gas opacity has a larger contribution.  \nWe show all the retrieved spectra for BD+60 1417B in Figure 5 for our winning fits. The retrieved forward model spectra consist of both the maximum-likelihood spectrum and median spectra. The maximum likelihood spectrum is the spectrum that maximizes the likelihood (i.e., best fit from every single iteration). In the retrieval, the median spectrum is the spectrum from the median value of each retrieved parameter which is put into the forward model. Utilizing both the maximumlikelihood and the median spectra explores the parameter space well and ideally these two retrieved fits should be near identical.  \nTable 4. List of models and thermal profiles tested in this work for BD+60 1417B and W0047 along with the corresponding ∆logEv. The lowest ∆logEv (e.g., 0.0) corresponds to the 'winning' model for each respective thermal parameterization for BD+60 1417B and W0047.  \nFigure 5a shows the retrieved spectrum for BD+60 1417B for the five point spline thermal profile. The retrieved spectrum shows good agreement through the J and H spectral bands. However, the retrieved spectrum shows a worse fit in the K -band which holds the temperature dependent CO molecular bands. The retrieved spectrum appears too faint for the K -band overall, which may explain our poorly constrained CO abundance in the five point spline thermal parameterization.  \nFigure 5b shows the retrieved spectrum for the Lavie profile for BD+60 1417B. We note that the retrieved spectrum is fainter than the flux peaks in the H spectral band which contains the gravity sensitive triangularband shape, this is in contrast to the five point spline thermal profile. Similarly, the K spectral band is poorly fit as well.  \nFigure 5c shows the retrieved spectrum for the Molli'ere thermal profile for BD+60 1417B. Similar to the Lavie retrieved profile (Figure 5b) the K band is too faint, which contains the CO band head features.", '8.4. Retrieved Gas Abundances and Fundamental Parameters:': "Figure B1 shows the posterior distributions for the radius, mass, metallicity, C/O, and T eff for the retrieval case of BD+60 1417B. The derived log(g) is calculated from the retrieved radius and mass. The scale factor  \n( R 2 /D 2 ) is calculated from the retrieved radius and parallax. T eff is then determined using this derived scale factor and by integrating the flux in the resultant forward model spectra extrapolated to cover 0.5 - 20 µ m.  \nDetermining the carbon-to-oxygen ratios (hereafter C/O) of substellar objects involves retrieving the following abundances ( X ):  \nC O = X CO + X CO 2 + X CH 4 X CO +2 X CO 2 + X H 2 O (5)  \nGiven our limited spectral coverage and retrieved abundances, when appropriate the C/O ratio is calculated under the assumption that all of the oxygen exists in H 2 O and CO and all of the carbon exists in CO. Our ratio is considered as a CO to H 2 O ratio. We note however that broader wavelength range retrievals for L dwarfs (e.g. Burningham et al. (2021) suggest that CO 2 abundances are higher than expected and may need to be taken into account for this total ratio.  \nC O = X CO X CO + X H 2 O (6)  \nThe following equations were used to derived a value for metallicity, [M/H]:  \nf H 2 = 0 . 84(1 -f gas ) (7a)  \nN H = 2 f H 2 N tot (7b)  \n4  \n<!-- image -->  \n(a)  \n<!-- image -->  \n(b)  \nFigure 3. Retrieved thermal profile (black line, light-red shading for 1 σ and 2 σ intervals) and cloud pressure for BD+60 1417B of the winning models for tested thermal profiles. Dashed lines show the condensation curves for possible cloud species. The lightgray region is the approximate photosphere. Self-consistent grid models from Sonora Diamondback (Morley et al. 2024) for SED derived temperature are plotted as solid colored lines. The median pressures of the extent of the cloud is shown in purple shading with gray shading indicating the 1 σ range, and the blue shading shows the vertical distribution of the deck cloud over which the optical depth drops to 0.5. (a) BD+60 1417B retrieved thermal profiles for five point spline, (b) Lavie profile, (c) Molli'ere 'Hybrid' profile  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \nFigure 4. Contribution functions for BD+60 1417B based on top-ranked model for each thermal parameterization case. The black shading shows the percentage contribution of each pressure to the flux at each wavelength. τ = 1 lines are included for gas phase opacities (cyan), and a cloud (magenta) (a) BD+60 1417B contribution function for five point spline thermal profile, (b) Lavie profile, (c) Molli'ere 'Hybrid' profile  \n<!-- image -->  \nFigure 6 shows the retrieved gas mixing ratios compared to predictions from thermochemical equilibrium models. The thermochemical equilibrium models are calculated and interpolated for our derived metallicity and C/O ratio for the three thermal profile parameterizations for BD+60 1417B where applicable. In the cases where the C/O ratio and [M/H] are unable to be determined (due to nonphysical results) we assume solar values (Asplund et al. 2009) - this is the case for the five point spline profile and the Molli'ere 'Hybrid' profile.  \n<!-- image -->  \n<!-- image -->  \nFigure 5. Retrieved forward model spectra for BD+60 1417B. The maximum-likelihood spectrum is shown in green and the median spectrum in yellow. The SpeX prism data is shown in black. (a) BD+60 1417B retrieved spectra for five point spline thermal profile, (b) Lavie profile, (c) Molli'ere 'Hybrid' profile.  \n<!-- image -->  \nN element = ∑ molecules n atom f molecule N tot (7c)  \nN M = ∑ element N element N H (7d)  \nwhere f H 2 is the H 2 fraction, f gas is the total gas fraction for all other gases, 0.84 represents the fraction of gas that is H 2 , N H is the number of neutral hydrogen atoms, N element is the number of atoms for each element of interest, n atom is the number of atoms for a given element contained in a single molecule and N tot is the total number of gas molecules.  \nIn our calculation of metallicity, N solar is determined using the same formula as N M , using the sum of the solar abundances (from Asplund et al. 2009) relative to H.  \n[M / H] = log N M N ⊙ (8)  \nIn Table 5, we show the retrieved gas abundances for the winning models for the three thermal profile parameterizations for BD+60 1417B.  \nIn Figure 6a, for the five point spline thermal profile major atmospheric constituents (CO and H 2 O) in the photosphere denoted by the gray box, we see that the gases are not predicted to vary greatly with altitude for our assumed solar abundances. However, our retrieved abundance for H 2 O is implausibly high and does not match thermochemical equilibrium predictions. For the FeH abundance, our uniform with altitude abundance is in disagreement with the model thermochemical abundance predictions, which is predicted to vary with altitude in the photospheric region. Conversely, we see that while the Na+K abundance is uniform with altitude in the atmospheres, our median retrieved abundance is depleted compared to the equilibrium predictions. Additionally, we see that our CO measurement for the five point spline thermal profile shows a depletion compared to predictions. This depletion could be the result of poor fits to the K band region of the spectrum for BD+60 1417B as shown in Figure 5a.  \nFigure 6. Retrieved uniform-with-altitude mixing abundances (solid lines) compared to thermochemical equilibrium grid abundances (dashed lines) for BD+60 1417B. Equilibrium predictions are calculated for assumed solar abundance or calculated based on derived [M/H] and C/O values. Shaded regions around the retrieved abundance values show our median retrieved values and 16th to 84th percentiles, respectively. (a) mixing abundances for the five point spline profile, (b) Lavie profile, (c) Molli'ere 'Hybrid' profile.  \n<!-- image -->  \nTable 5. Comparison off retrieved gas abundances for the preferred models for BD+60 1417B and W0047 for our various thermal profile parameterizations  \nIn Figure 6b we show our retrieved gas mixing ratios along with predictions from thermochemical equilibrium models interpolated for our derived metallicity and C/O ratio for the Lavie thermal profile. For H 2 O, we see a slight disagreement ( > 1 σ ) in the photospheric region for our retrieved abundance in comparison to the thermochemical equilibrium models. Similar to the five point spline profile, we observe a depletion of predicted CO in the photosphere region that might be due to the poor fits in the K band region. The CO is also predicted to vary with altitude in the photosphere, which cannot be accounted for in our uniform-with-altitude abundance retrieval model. The retrieved FeH abundance shows photospheric agreement near the 10 bar region. While an overall poor constraint, for the CO 2 retrieved abundances we see an agreement with the predicted abundance that is within the 1 σ error bar. The  \nNa+K retrieved abundances are depleted compared to thermochemical equilibrium models.  \nFor the Molli'ere thermal profile, similar to the five point spine thermal profile, the major atmospheric constituents (CO and H 2 O) in the photosphere, we see that the gases are not predicted to vary greatly with altitude for our assumed solar abundances. We again see an implausibly high retrieved H 2 O abundance. For FeH, our retrieved abundance is consistent with the predicted abundance at the bottom of the photosphere at 10 bar. The Na+K retrieved abundance is not predicted to vary with altitude and the values for the thermochemical equilibrium models within the 1 σ ranges for our retrieved values  \nTable 7 shows the comparison of the fundamental and derived parameters for BD+60 1417B and the analog object W0047 for the three different thermal profile parameterizations. We find that for all three of our thermal parameterization checks for BD+60 1417B our fundamental parameters are in disagreement with the SED parameters. We discuss this in more detail in Section 10.3.", '8.5. Cloud Properties': "We find that the spectrum for BD+60 1417B is best described by a cloudy model regardless of thermal profile parameterization. This remains our strongest retrieval conclusion in the absence of being able to definitively converge on either a P-T profile or reasonable chemical abundances. A summary of cloud parameters for BD+60 1417B is shown in Table 6. We note a positive slope in the contribution function for the five point spline thermal profile which can be indicative of large particles (1 µ m) in the atmosphere from Mie theory (e.g. Lacy & Burrows 2020). Where appropriate the α parameter is positive for the power law.  \nThe pressure of the cloud in the atmosphere of BD+60 1417B is shown in various shaded regions to the right of Figure 3. The notable difference is that the pressure at which τ > 1 for the Lavie profile is higher than those for the five point spline and Molli'ere 'Hybrid' thermal profile. For BD+60 1417B the cloud height is similar for the five point spline and Lavie profile, while the Molli'ere 'Hybrid' Profile cloud location extends further in the photosphere.", '9. W0047 RETRIEVAL RESULTS: A COMPARATIVE EXAMPLE': 'The original aim of this work was to use retrieval analysis and the Brewster framework to understand the newly discovered young low surface gravity planetarymass companion BD+60 1417B. However, we found that  \ndespite testing three different thermal parameterizations from the literature, often times not only did our fundamental parameters not match the SED analysis but we arrived at implausibly high or unconstrained molecular abundances for this benchmark object.  \nIn Faherty et al. (2021), the isolated L-dwarf W0047 was found to be a close match in spectral morphology and fundamental characteristics to BD+60 1417B (Table 1). W0047 is an isolated object with higher S/N SpeX prism observations (S/N ∼ 5 vs S/N ∼ 13) than BD+60 1417B, and slightly brighter apparent magnitude. We aim to use the same retrieval framework on this isolated, clone object to see if more robust results are obtained on higher quality data. Table 4 shows the list of models tested for W0047, along with the number of parameters for each model and the resulting ∆logEv.', '9.1. Best-fit Models': 'The power law deck cloud is preferred by all three thermal parametrizations tested. For this comparative object to BD+60 1417B it is again evident that based on our selection criterion, cloudless models are strongly rejected for all three different thermal parameterizations. Once again this emphasizes our strongest retrieval conclusions: W0047 is a cloudy L dwarf. W0047 is known to be variable in the infrared, likely driven by nonuniform condensate clouds (Lew et al. 2016; Vos et al. 2018) thus our retrieval conclusion is not unexpected.', '9.2. Thermal Profiles': "Figure 7 shows the retrieved thermal profiles and cloud pressures for the winning models compared to grid models and condensation curves for W0047. For the five point spline (Figure 7a), the retrieved thermal profile is not in agreement with the grid models, as was seen for BD+60 1417B. In the photosphere, the retrieved thermal profile is ∼ 300 - 400K warmer than the grid models. The shape of the retrieved thermal profile is nearly linear throughout the photosphere and the extrapolated regions above and below the photosphere. The grid model that has the closest fit to retrieved thermal profile in the photosphere is the f sed = 1 model. The median cloud location 1 σ range is within the photosphere.  \nFor the Lavie thermal profile parameterization (Figure 7b), we find that the retrieved thermal profile is not well aligned with the grid models. The cloud is located within the observable photosphere. Below the observable photosphere, the closest fitting model is f sed = 1.  \nFigure 7c shows the retrieved thermal profile for the Molli'ere thermal parameterization. The models ( f sed = 1 - 3) are generally not in agreement with the retrieved  \nTable 6. Summary of retrieved cloud properties for the preferred models of respective thermal parameterization for BD+60 1417B and W0047.  \nprofile in the photosphere region. Above the photosphere in the extrapolated region, the models diverge more from the extrapolated thermal profile. For the Molli'ere profile, the cloud deck is prevalent within the photosphere.", '9.3. Contribution Functions and Retrieved Spectra': "For the retrieved spectra for W0047 (Figure 8a), the five point spline approach underestimates the flux in the H and K band peaks hence the CO bands ( ≈ 2.2 2.4 µ m) are not fit well. This is a similar conclusion to what was found for BD+60 1417B. Figure 8b shows the retrieved spectra for W0047 for the Lavie thermal profile. The retrieved spectrum generally traces the observed spectra and has a better fit to the CO feature in the K band and does not fit the flux peaks of the triangular H band region as with BD+60 1417B. Similar to the five point spline profile the Molli'ere thermal profile retrieved spectra underestimates and poorly fits the flux in the H and K bands. For all three thermal parameterizations, the retrieved spectra underestimates the flux peak in the J -band.  \nFigure 9 shows the contribution functions for W0047 for the winning cloud deck power law models. For the five point spline thermal profile, the bulk of the flux between ≈ 10 -2 and 1 bar contributes to the observed spectrum. Across the spectrum for the five point spline thermal profile (Figure 9a), the cloud opacity dominates, except ≈ 1.4 µ mregion, where H 2 Oopacity contributes significantly. For the Lavie thermal profile, the bulk of the flux between ≈ 0.003 and 10 bars contributes to the observed spectrum. The cloud opacity has a large contribution across the full spectrum with the exception of a small region in the H band ( ≈ 1.8 - 2.0 µ m), where gas opacity becomes more important. Notably, the observable photosphere extends to the lowest pressures of the three thermal parameterizations that were tested.  \nFor the Molli'ere thermal profile, the cloud opacity dominates the spectrum except for the ≈ 1.4 µ m region where the gas opacity is dominant. The contribution function for the Molli'ere thermal profile is similar to that of the five point spline where the bulk of the flux between ≈ 10 -2 and 1 bar contributes to the observed spectrum.  \n(a)  \n<!-- image -->  \n(b)  \n<!-- image -->  \nFigure 7. Retrieved thermal profile (black line, light-purple shading for 1 σ and 2 σ intervals) and cloud pressure for W0047 for a power law cloud deck for the Lavie thermal profile. Dashed lines show the condensation curves for possible cloud species. The lightgray region is the approximate photosphere. Self-consistent grid models from Sonora Diamondback profiles for the SED derived temperature are shown as solid colored lines (Morley et al. 2024). The median pressures of the extent of the cloud is shown in purple shading with gray shading indicating the 1 σ range, and the blue shading shows the vertical distribution of the deck cloud over which the optical depth drops to 0.5. (a) W0047 thermal profiles for five point spline, (b) Lavie profile, (c) Molli'ere 'Hybrid' profile.  \n<!-- image -->  \nFigure 9. Contribution functions for W0047 based on topranked model for each thermal parameterization cases. The black shading shows the percentage contribution of each pressure to the flux at each wavelength. τ = 1 lines are included for gas phase opacities (cyan), and a cloud (magenta) (a) W0047 contribution function for five point spline thermal profile, (b) Lavie profile, (c) Molli'ere 'Hybrid' profile  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \nFigure 8. Retrieved forward model spectra for the power law deck models for W0047. The maximum-likelihood spectrum is shown in orange and the median spectrum in green. The SpeX prism data is shown in black. (a) W0047 retrieved spectra for the five point spline thermal profile, (b) Lavie profile, (c) Molli'ere 'Hybrid' profile  \n<!-- image -->", '9.4. Retrieved Gas Abundances and Fundamental Parameters': "Figure 10 shows the retrieved gas mixing ratios compared to predictions from thermochemical equilibrium models for W0047. As before, the thermochemical equilibrium models are calculated and interpolated for our derived metallicity and C/O ratio for the three thermal profile parameterizations. In the cases where the C/O ratio cannot be determined (due to implausibly high abundances) we assume solar values (Asplund et al. 2009).  \nFor the five point spline thermal model (Figure 10a), the thermochemical equilibrium models predict many of the included gases to remain relatively constant with altitude (e.g., H 2 O, CO, CO 2 , Na+K) in the photosphere. However, FeH and CrH are predicted to vary within the photosphere (denoted by the gray box). As with the retrieved abundances for BD+60 1417B for the five point spline, we retrieve an implausibly high H 2 O abundance that does not match thermochemical equilibrium predictions. Similarly, our retrieved CO abundance is depleted compared to predicted values. In Figure B2, we show the posterior distribution for retrieved and derived parameters for W0047.  \nWe show the retrieved gas mixing ratios compared to the thermochemical equilibrium abundances for the Lavie profile (Figure 10b). For this thermal profile parameterization, the retrieved abundances are more constrained, compared to the five point spline and Molli'ere thermal profiles, however none of the retrieved abundances match the predicted thermochemical equilibrium abundances.  \nSimilar to the five point spline thermal profile, for the Molli'ere 'Hybrid' profile retrieved abundance (Figure 10c) we find an implausibly high H 2 O abundance with respect to the predicted abundances assuming solar composition. While not a tight constraint, our CO abundance 1 σ values are within the limits of the model values. FeH and CrH are predicted to vary with altitude, and show agreement with the retrieved values near the 10 bar region. Our retrieved Na+K values are not within the predicted range within the photosphere.  \nFigure 10. Retrieved uniform-with-altitude mixing abundances as compared to thermochemical equilibrium grid abundances for W0047. Equilibrium predictions are shown as dashed lines. The solid straight lines and shading show our median retrieved values and 16th to 84th percentiles, respectively. (a) five point spline thermal profile, (b) Lavie profile , (c) Molli'ere 'Hybrid' profile  \n<!-- image -->", '9.5. Cloud Properties': "We find that the spectrum for W0047 is best described by a cloudy model regardless of thermal profile parameterization. A summary of cloud parameters for W0047 is shown in Table 6. The α parameter is positive for the power law, with the lowest α present for the Lavie profile.  \nThe pressure of the cloud in the atmosphere of W0047 is shown in various shaded regions to the right of Figure 7. The notable difference is that the pressure at which τ > 1 for the Lavie profile is higher than those for the five point spline and Molli'ere 'Hybrid' thermal profiles, indicating a higher altitude cloud and shallower pressure of the corresponding contribution function (Figure 9b).", '10.1. The Complexity of Low Surface Gravity L-dwarfs': "Young L-dwarfs are redder at near-infrared wavelengths than field objects (Faherty et al. 2013a, 2016), in part due to physical conditions brought on by a lower surface gravity. High altitude silicate clouds as a result of inefficient dust settling is a likely cause for low surface gravity object reddening (Looper et al. 2008; Su'arez & Metchev 2022, 2023). The atmospheric physics and chemistry of low surface gravity L-dwarfs brings about changes to the features and near-infrared spectral shapes that can complicate modeling and analysis (Bonnefoy et al. 2013; GRAVITY Collaboration et al. 2020; Whiteford et al. 2023).  \nComplexities in the retrievals of young L-dwarfs like those seen within this work probably arise from compounding factors involving the emergence of condensate clouds. As noted by Rowland et al. (2023), there exist three main challenges in the near-infrared modeling of L-dwarfs: (1) the limited range of the atmosphere that the near-infrared region can probe, which can limit our ability to constrain thermal profiles (2) the different opacity sources and absorbers which prefer uniform or non-uniform with altitude abundances, (3) the presence of clouds that can cause degeneracies with thermal profiles.  \nNotably, Burningham et al. (2021) found that chemical disequilibrium and non-solar elemental abundances are needed to accurately capture the complexity of cloudy L-dwarfs in future cloud models for retrievals. Related to this regime as well for the directly-imaged planets in the HR 8799 system, Lavie et al. (2017) commented that their future retrieval studies could also provide opportunities to explore disequilibrium chemistry through vertical mixing as well.  \nOne clear conclusion from this work is that regardless of P-T profile chosen, both BD+60 1417B and W0047 retrieval approaches strongly preferred the cloudy model. However what was not clear was the parameterization of said cloud with the data on hand. It remains unclear how the interplay of clouds, equilibrium chemistry, and gravity sensitive features impact the spectra for red young L-dwarfs. Similarly it is unclear whether the full physical picture of all those parameters is fully represented within the Brewster retrieval framework. Hence performing the retrieval and achieving a physical or non-physical thermal profile and abundance output is itself a test on whether we have included the full range of needed physics in order to model these complex objects.", "10.2. 'Free Retrievals' vs Chemical Equilibrium Assumption": "Due to consistent, implausibly high abundances retrieved for BD+60 1417B and W0047 we test our retrievals assuming thermochemical equilibrium, which has previously been demonstrated by Burningham et al. (2021) and Gonzales et al. (2022). Instead of retrieving individual abundances, this method retrieves [M/H] and C/O. For simplicity, we test the chemical equilibrium assumption on the 'winning' model for the Molli'ere 'Hybrid' profile for BD+60 1417B.  \nIn this model, we find the following parameters: [M/H] = 1.83 +0 . 10 -0 . 16 , C/O = 0.38 +0 . 16 -0 . 09 , R Jup = 0.51 ± 0.01 and M Jup = 49.66 +16 . 78 -16 . 79 . The chemical equilibrium assumption yields a smaller radius and higher mass than with our free retrieval. While the C/O and [M/H] might be more realistic under these assumptions our models prefer the free retrieval assumption over the chemical equilibrium assumption; ∆logEv = 0 and 6.74 respectively. These results are consistent with Burningham et al. (2021), who also found that free retrievals were preferred for a higher gravity, cloudy L dwarf.  \nFor the chemical equilibrium test (Figure A1a), we find that the retrieved P-T profile is not isothermal in the photosphere unlike the 'free' retrieval runs. In the photosphere, the retrieved P-T profile closely traces the grid model with fsed = 2. The 1 σ median cloud is within the observable photosphere. Notably, we find that the chemical equilibrium model provides a better fit to the K -band region containing the CO band head.  \nWe find that while the chemical equilibrium model does not provide a statistically better fit (∆logEv) for BD+60 1417B, it does provide a more plausible thermal profile. As we find that the chemical-equilibrium model (which allows abundances to vary, but in a fixed way) improves certain aspects of the retrieval; incorpo-  \nting vertically-varying abundances into free-retrievals with Brewster can prove to be a promising next step to take (Rowland et al. 2023).", '10.3. Comparison with SED-derived Fundamental Parameters': "We find discrepancies for the fundamental parameters derived and extrapolated for both BD+60 1417B and W0047 regardless of thermal model parameterization compared to our SED-derived values. Our retrieved radii for BD+60 1417B and W0047 suffer from the 'small-radius problem' known to be present in retrieval analysis for L and L/T transition objects (Kitzmann et al. 2020; Molli'ere et al. 2020; Gonzales et al. 2020; Gonzales et al. 2021; Gonzales et al. 2022; Burningham et al. 2021; Lueber et al. 2022; Vos et al. 2023). The retrieved radii for W0047 were consistently ∼ 0.5 R Jup regardless of the thermal profile parameterization used. For BD+60 1417B, the Molli'ere 'Hybrid' profile has an upper limit of 0.86 R Jup , the highest of the derived radii, but still inconsistent with the SED semi-empirical evolutionary model of 1.29 ± 0.06 R Jup . The high surface gravity for both BD+60 1417B and W0047 is likely driven by the unrealistically small radii. The 'small radii' problem is predominantly seen in L dwarfs therefore it may be linked to poorly parameterized cloud properties given that condensates are a dominating and difficult to model component of their spectra. The high surface gravity is persistent regardless of the thermal profile tested - also likely due to degeneracies between radius and surface gravity or a known issue with metallicity as discussed in Gonzales et al. (2021)  \nAcross the board of thermal parameterizations, the retrieved mass for BD+60 1417B is larger than the SED semi-empirical analysis (13.47 ± 5.3 M Jup ). The Molli'ere 'Hybrid' profile produced the lowest mass of (36.54 +13 . 94 -12 . 61 M Jup ), however, this is more than 5 σ larger than the SED derived mass. A similar observance is found for W0047. This is likely due to the 'small' radius problem driving a large surface gravity and in turn forcing a larger mass.  \nOur derived T eff is consistently lower than the SED derived values for both BD+60 1417B and W0047. Typically the 'small radius problem' tends to push T eff higher -whereas we see lower T eff in both BD+60 1417B and W0047. Our low-luminosity could be due to our high retrieved H 2 O abundance which in turn could be blocking flux for both BD+60 1417B and W0047.", 'BD+60 1417B': 'Table 7. Comparison of Fundamental and Derived Parameters for BD+60 1417B and W0047 winning models for each thermal profile parameterizaton', 'W0047': 'a Values are from the fundamental parameters derived from SED analysis (see § 4)', '10.4. Host Star Chemistry Implications on the Secondary': "Recently, Calamari et al. (2024) introduced a theoretical chemical framework to predict cloud conditions of brown dwarfs based on FGK host stars chemical abundances. Using the approach of Calamari et al. (2024), we utilized the chemical abundances of the host star BD+60 1417 discussed in Section 5 in order to assess the likely cloud composition of BD+60 1417B. Compared to other stars in the local solar neighborhood (e.g. Brewer et al. 2016) we find that BD+60 1417 has a relatively high Mg/Si ratio and a subsolar C/O ratio. The implications of this chemistry would be on the distribution of refractory elements into particular cloud phases (i.e. enstatite, forsterite, or quartz) and the percentage of oxygen into such silicate clouds. According to Calamari et al. (2024), using the BD+60 1417 elemental abundances of Mg, Si, C, and O we anticipate the formation of forsterite (Mg 2 SiO 4 ) and enstatite (MgSiO 3 ) clouds in the atmosphere of the companion. If this companion is indeed relatively silicon depleted, we would not expect a best fit model containing a quartz cloud. As is shown in Burningham et al. (2021) and Vos et al.  \n(2023), we also expect to see a large portion of the Fe inventory to condense as a patchy Fe deck cloud under the silicate cloud layer. Additionally, based on our host star chemistry, we can estimate the percentage of oxygen sunk into these silicate clouds as done in Calamari et al. (2024). For this metal-rich system, we calculate an oxygen sink percentage of approximately 24%, which lies on the higher end of the local neighborhood distribution shown in Figure 2 in Calamari et al. (2024). In future modeling, we would expect this to have implications on retrieved versus bulk C/O ratio.  \nFinally, based on the results of Su'arez & Metchev (2022), we may also expect to see more enstatite than forsterite given that this source is also a low gravity late-type L dwarf. Consequently, a longer wavelength spectrum such as that enabled by JWST's MIRI spectrograph would yield strong insights into the atmosphere and ties to the host stars composition.", '10.5. Predictions for the Mid-Infrared Region': "With a limited number of planetary-mass worlds observed with JWST (Miles et al. 2022), we have learned a wealth of information from their extended mid-infrared coverage, enabled by the MIRI coverage. In an effort to understand the mid-infrared region we utilize the best fit retrieved spectra from our three thermal profiles for both BD+60 1417B and W0047 (Figure 5 & Figure 8 respectively) to produce extended spectra predictions from 1 to 20 µ m. As stated, while we can not constrain abundances with confidence in this work, we strongly predict the presence of clouds in both objects. We take the parameters for the max-likelihood in each best fit retrieval for a given P-T profile and rerun our forward model using the parameters over a wider wavelength range. The extrapolated spectra for BD+60 1417B and W0047 are shown in Figure 11 and serve as predictive power for what might be seen in a JWST MIRI spectrum.  \nWe examine the extended retrieved spectrum relative to available photometric points for each source to compare/contrast the best P-T profile predictions. In the case of BD+60 1417B we find that for both the five point spline and the Molli'ere 'Hybrid' profiles, the predicted extended spectra underestimate the flux from the WISE photometry measurements. The extrapolated spectrum for the Lavie profile most closely traces the WISE photometry, particularly at the W2 photometric point. For W0047, we find that the extracted spectra again underestimates the flux from WISE photometry measurements for both the five point spline and Molli'ere 'Hybrid' profiles. In contrast, the Lavie profile extracted spectrum appears to slightly overestimate the flux compared to the WISE photometry, with the expected of  \nW1 photometric point. While no single thermal profile prediction appears to match current data, it is clear that the JWST MIRI spectrum has the potential to differentiate the P-T profile and parameterize cloud properties.  \n<!-- image -->  \nFigure 11. Extracted mid-infrared spectra for BD+60 1417B and W0047 with important molecular gas highlighted.  \n<!-- image -->", '11. SUMMARY AND CONCLUSION': "In this work we present an atmospheric retrieval analysis of two young red L-dwarfs with similar spectral morphology and fundamental parameters, BD+60 1417B and W0047, using the Brewster retrieval framework. We also present abundances for the K0 host star, BD+60 1417, from high resolution spectroscopy from PEPSI/LBT (Section 5).  \nWith the PEPSI data, we constrain an iron abundance of 0.27 ± 0.03 dex. We constrain measurements of the C/O = 0.23 ± 0.12 and Mg/Si = 1.41 ± 0.19. Given these parameters, we find that forsterite and enstatite clouds are expected to form and not quartz clouds.  \nIn our retrieval analysis we test a number of different models and thermal parameterizations (five point spline, Lavie profile and a Molli'ere 'Hybrid' profile) and find that both BD+60 1417B and W0047 are best described by a cloudy atmosphere regardless of the chosen P-T profile. Regardless of thermal profile parameterization,  \nthe retrieved profile was generally not in agreement with forward models for both BD+60 1417B and W0047. We find that for both BD+60 1417B and W0047 the retrieved spectrum from Brewster had difficulty fitting the K -band region, which contains the temperature sensitive CO feature.  \nLastly, although we test three different thermal profile parameterizations we find that we are unable to reliably constrain the C/O ratio of these planetary-mass objects due to retrieving implausibly high abundances (e.g H 2 O). Generally we find that our abundances are not in agreement with thermochemical equilibrium models. While we experience challenges in modeling these young low-gravity objects, we find that improvements can be made in terms of modeling approach as we discover that non-uniform with altitude abundances provide a more realistic thermal profile for BD+60 1417B. This work continues to provide technical lessons in retrieval analysis of young cloudy L-dwarfs.  \nJWST has the capabilities to provide wider wavelength coverage and simultaneous higher resolution and signal-to-noise observations. Wide wavelength coverage will allow us to probe the pressure-temperature profile from (0.1 - 1 bar) through mid-infrared all the way to 10 bar pressure with near-infrared coverage. Near-infrared and mid-infrared data is crucial to retrieve a reliable PT profile for early-to-mid L dwarfs (Burningham et al. 2021; Vos et al. 2023; Rowland et al. 2023). Future retrievals of BD+60 1417B and W0047 with JWST data would allow us to accurately constrain the P-T profile and cloud parameters to break degeneracies between clouds and composition.  \nFacilities: Exoplanet Archive, LBT (PEPSI)  \nSoftware: Brewster (Burningham et al. 2017), PyMultiNest (Buchner et al. 2014) , corner.py (Foreman-Mackey et al. 2013)", 'ACKNOWLEDGEMENTS': "We thank our anonymous reviewer for providing detailed recommendations that significantly improved the quality of this paper.  \nThis research has made use of the NASA Exoplanet Archive, which is operated by the California Institute of Technology, under contract with the National Aeronautics and Space Administration under the Exoplanet Exploration Program.  \nNASA's Astrophysics Data System Bibliographic Services together with the VizieR catalogue access tool and SIMBAD database operated at CDS, Strasbourg, France, were invaluable resources for this work. This publication makes use of data products from the Two  \nMicron All Sky Survey, which is a joint project of the University of Massachusetts and the Infrared Processing and Analysis Center/California Institute of Technology, funded by the National Aeronautics and Space Administration and the National Science Foundation.  \nThis work benefited from involvement in ExoExplorers, which is sponsored by the Exoplanets Program Analysis Group (ExoPAG) and NASA ' s Exoplanet Exploration Program Office (ExEP). Caprice Phillips thanks the LSSTC Data Science Fellowship Program, which is funded by LSSTC, NSF Cybertraining Grant #1829740, the Brinson Foundation, and the Moore Foundation; her participation in the program has benefited this work. E.J.G. is supported by an NSF Astronomy and Astrophysics Postdoctoral Fellowship under award AST-2202135. E.G. acknowledges support from the Heising-Simons Foundation for this research. J. M. V. acknowledges support from a Royal Society Science Foundation Ireland University Research Fellowship (URF \\ 1 \\ 221932). BB acknowledges support from UK Research and Innovation Science and Technology Facilities Council [ST/X001091/1].  \nThis work benefited from the 2023 Exoplanet Summer Program in the Other Worlds Laboratory (OWL) at the University of California, Santa Cruz, a program funded by the Heising-Simons Foundation and NASA.  \nThis work has made use of data from the European Space Agency (ESA) mission Gaia (https://www. cosmos.esa.int/gaia), processed by the Gaia This publication makes use of data products from the Widefield Infrared Survey Explorer, which is a joint project of the University of California, Los Angeles, and the Jet Propulsion Laboratory/California Institute of Technology, funded by the National Aeronautics and Space Administration. Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/ gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement  \nThis work is based on observations made with the Large Binocular Telescope. The LBT is an international collaboration among institutions in the United States, Italy and Germany. LBT Corporation partners are: The University of Arizona on behalf of the Arizona Board of Regents; Istituto Nazionale di Astrofisica, Italy; LBT Beteiligungsgesellschaft, Germany, representing the Max- Planck Society, The Leibniz Institute for Astrophysics Potsdam, and Heidelberg University; The Ohio State Uni- versity, representing OSU, University of Notre Dame, University of Minnesota and University of Virginia.  \nThe data used in this publication were collected through the MENDEL high performance computing (HPC) cluster at the American Museum of Natural History. This HPC cluster was developed with National Science Foundation (NSF) Campus Cyberinfrastructure support through Award 1925590. This work  \nwas also supported by the National Science Foundation via awards 1909776, 1614527 and AST-1909837.  \nComputation reported in this work was also carried out on the Unity cluster of the College of Arts and Sciences at the Ohio State University. The computational resources and support provided are gratefully acknowledged.", 'A. CHEMICAL EQUILIBRIUM ASSUMPTION MODEL': "In Figure A1, we present the retrieval plots for the chemical equilibrium test for the winning Molli'ere 'Hybrid' profile for BD+60 1417B. While this model is not preferred over the 'free' retrieval we discuss the implications of testing this retrieval set up.  \n<!-- image -->  \n<!-- image -->  \n(a)  \nFigure A1. (a) retrieved thermal profile for the winning Molli'ere 'Hybrid' profile for BD+60 1417B under the assumption of chemical equilibrium. (b) retrieved forward model spectra for BD+60 1417B for chemical equilibrium (c) contribution function for BD+60 1417B for the winning model  \n<!-- image -->", 'B. CORNER PLOTS': "In Figure B1 and B2, we present the posteriors for the gas abundances and derived fundamental parameters for the three different thermal profile parameterizations used in our retrievals for BD+60 1417B and W0047 respectively.  \nFigure B1. Best Fit Models for the five point spline thermal profile, Lavie profile, and Molli'ere profile for BD+60 1417B: power law cloud deck (blue), grey cloud deck (pink), and power law cloud deck (silver) for best fits respectively. Posterior probability distributions for retrieved gas abundances, mass, radius and derived quantities of the best fit models. Far right diagonal plots show marginalized posteriors for each parameter along with 2D parameter correlation histograms. Dashed lines show the median retrieved value (reported value) along with 1 σ confidence intervals. Gas abundances are in units of dex.  \n<!-- image -->  \n+0.24  \n0.02  \n+0.12  \n0.14  \nFigure B2. Best Fit Models for the five point spline thermal profile, Lavie profile, and Molli'ere profile for W0047: power law cloud deck (gold), power law cloud deck (purple), and power law cloud deck (green) for best fits respectively. Posterior probability distributions for retrieved gas abundances, mass, radius and derived quantities of the best fit models. Far right diagonal plots show marginalized posteriors for each parameter along with 2D parameter correlation histograms. Dashed lines show the median retrieved value (reported value) along with 1 σ confidence intervals. Gas abundances are in units of dex.  \n<!-- image -->"}