Atmospheric Retrieval of L Dwarfs: Benchmarking Results and Characterizing the Young Planetary Mass Companion HD 106906 b in the Near-infrared

<jats:title>Abstract</jats:title>                <jats:p>We present model constraints on the atmospheric structure of HD 106906 b, a planetary-mass companion orbiting at a ∼700 au projected separation around a 15 Myr old stellar binary, using the APOLLO retrieval code on spectral data spanning 1.1–2.5 <jats:italic>μ</jats:italic>m. C/O ratios can provide evidence for companion formation pathways, as such pathways are ambiguous both at wide separations and at star-to-companion mass ratios in the overlap between the distributions of planets and brown dwarfs. We benchmark our code against an existing retrieval of the field L dwarf 2MASSW J2224–0158, returning a C/O ratio consistent with previous fits to the same <jats:italic>JHK</jats:italic>                   <jats:sub>                      <jats:italic>s</jats:italic>                   </jats:sub> data, but disagreeing in the thermal structure, cloud properties, and atmospheric scale height. For HD 106906 b, we retrieve C/O <jats:inline-formula>                      <jats:tex-math> <?CDATA $={0.53}_{-0.25}^{+0.15}$?> </jats:tex-math>                      <mml:math xmlns:mml="http://www.w3.org/1998/Math/MathML" overflow="scroll">                         <mml:mo>=</mml:mo>                         <mml:msubsup>                            <mml:mrow>                               <mml:mn>0.53</mml:mn>                            </mml:mrow>                            <mml:mrow>                               <mml:mo>−</mml:mo>                               <mml:mn>0.25</mml:mn>                            </mml:mrow>                            <mml:mrow>                               <mml:mo>+</mml:mo>                               <mml:mn>0.15</mml:mn>                            </mml:mrow>                         </mml:msubsup>                      </mml:math>                      <jats:inline-graphic xmlns:xlink="http://www.w3.org/1999/xlink" xlink:href="ajacfb87ieqn1.gif" xlink:type="simple" />                   </jats:inline-formula>, consistent with the C/O ratios expected for HD 106906's stellar association and therefore consistent with a stellar-like formation for the companion. We find abundances of H<jats:sub>2</jats:sub>O and CO near chemical equilibrium values for a solar metallicity but a surface gravity lower than expected, as well as a thermal profile with sharp transitions in the temperature gradient. Despite high signal-to-noise ratio and spectral resolution, more accurate constraints necessitate data across a broader wavelength range. This work serves as preparation for subsequent retrievals in the era of JWST, as JWST's spectral range provides a promising opportunity to resolve difficulties in fitting low-gravity L dwarfs and also underscores the need for simultaneous comparative retrievals on L-dwarf companions with multiple retrieval codes.</jats:p>

section_content: Mass is typically taken as the discriminator between planets and brown dwarfs, based on the minimum of ∼13 M J needed for sustained deuterium fusion (Spiegel et al. 2011). While one can use mass alone to define the classes of planets and brown dwarfs, there is an alternate definition based on the formation pathway of an object as more "star"-like or "planet"-like (see, e. g., Janson et al. 2012; Currie et al. 2014 Currie et al., 2020;; Pepe et al. 2014; Schlaufman 2018 ). These definitions may produce similar categories of planetary and brown dwarf companions as with the mass definition. However, 13 M J is not known to be a strict upper limit to forming companions as planets (i. e., that the companion forms within a circumstellar disk surrounding a young star)-nor is 13 M J a strict minimum below which objects may not collapse from a molecular cloud. In exoplanet and brown dwarf demographics, there is a local minimum in the observed distributions of companions' masses as ratios to their hosts' masses, as seen in radial velocity and astrometry (Sahlmann et al. 2010), direct imaging (Reggiani et al. 2016; Vigan et al. 2017; Nielsen et al. 2019), and microlensing (Shvartzvald et al. 2016; Suzuki et al. 2016). That is, using the mass definition of planets and brown dwarfs, it is both difficult to form planets at masses as large as ∼1% of their hosts and difficult to form brown dwarfs with mass ratios that small. It is companions in this region that serve as the more ambiguous cases when using formation history as the criterion for distinguishing planets and brown dwarfs. How can we tell the formation pathway for individual companions?

The chemical composition of a companion reflects its formation pathway, especially in the carbon-to-oxygen (C/O) ratio of its envelope relative to those measured in its host star. 10 lanetary formation pathways are themselves generally divided into core accretion versus gravitational instability (see, e. g., Forgan & Rice 2013; Forgan et al. 2015; Kratter & Lodato 2016; Forgan et al. 2018). Core accretion allows a planet's C/O ratio to diverge from its host's C/O based on where and when companions accrete their material in the disk. The H 2 O, CO, and CO 2 ice lines determine the relative fraction of C and O contained in gases versus solids as a function of distance from the host star (e. g., Mousis et al. 2009; Öberg et al. 2011; Madhusudhan et al. 2011). Since disk chemistry also evolves in time, planet compositions will reflect the chemical evolution in the disk over their development (e. g., Booth et al. 2017; Madhusudhan et al. 2017), both in the disk midplane (Eistrup et al. 2016 (Eistrup et al., 2018) ) and vertically (Cridland et al. 2020). For planets formed via gravitational instability the formation mechanism is different, but there is still ample opportunity for the atmosphere to evolve its chemistry away from a stellar-like C/O ratio (in this case, because the protoplanetary fragment has time to stratify its C and O compounds between the core and envelope; see, e. g., Ilee et al. 2017). 

While a nonstellar C/O ratio will certainly be reflected in the companion's observable emission spectrum, the presence of clouds in the companion photosphere requires a careful modeling approach. Many young (100 Myr) companions in the target mass ratio range fall in the L and T spectral types (Kirkpatrick 2005). In the warmer L dwarfs (1300 K  T eff  2000 K), a variety of cloud species become important opacity sources (see, e. g., Morley et al. 2012; Marley et al. 2013; Helling & Casewell 2014; Helling 2021) ; silicates such as enstatite (MgSiO 3 ) and forsterite (Mg 2 SiO 4 ), iron (Fe), aluminum oxides (e. g., Al 2 O 3 ), and quartz (SiO 2 ) can all contribute significantly in column density to L-dwarf photospheres (Helling & Woitke 2006; Gao et al. 2020; Woitke et al. 2020; Burningham et al. 2021). One important limitation in using observed gas abundances alone to constrain a C/O ratio is that oxygen-rich cloud species can condense out a significant amount of the oxygen budget at the pressures at which they reside. This biases the gas-derived C/O constrained from an emission spectrum, as it will be carbon-rich relative to the cumulative envelope C/O of the companion at the time of formation. We discuss our results in the context of this assumption in Section 7. 

Fitting brown dwarf spectra has traditionally relied on interpolations using grids of forward models that rely on specific input physics (e. g., Burrows & Liebert 1993; Allard et al. 1996; Marley et al. 1996; Tsuji et al. 1996; Saumon et al. 2000; Geballe et al. 2001; Hubeny & Burrows 2007; Saumon & Marley 2008; Yamamura et al. 2010; Patience & King 2012). There have been numerous analyses of field brown dwarfs that use such libraries of model spectra to constrain global properties such as effective temperature, metallicity, age, surface gravity, luminosity, mass, and in some cases cloud layers (e. g., Allers et al. 2007; Cushing et al. 2007; Cruz et al. 2009; Stephens et al. 2009; Rice et al. 2010; Allers & Liu 2013; Bonnefoy et al. 2014; Martin et al. 2017). The link between theory and observation for substellar atmospheres has been evolving for quite some time, as reviewed in works such as Burrows et al. (2001), Marley et al. (2013), and Marley & Robinson (2015), which has motivated a second approach to spectral fitting, namely atmospheric retrieval. 11 Retrievals opt to generate forward models in parallel with a parameter estimation technique such as Markov Chain Monte Carlo (MCMC) or nested sampling, rather than interpolate from a precomputed grid. Such an approach is computationally intensive and typically requires one to make simplifying assumptions in the parameterization that may or may not be physically consistent. However, the potential benefit is a more direct and precise constraint on key physical parameters, which may be warranted if the spectra are sensitive to small changes in the parameters, such as those of the temperature-pressure (T-P) profile, molecular abundances, or significant cloud opacities. 

We now have the results of a growing number of L-dwarf retrievals to guide us in interpreting our data. Burningham et al. (2017) retrieve atmospheric properties from the near-infrared spectra of two mid-L field dwarfs using the Brewster retrieval code. This is then expanded into the mid-infrared in Burningham et al. (2021), constraining multiple cloud species including enstatite (MgSiO 3 ), quartz (SiO 2 ), and iron (Fe). Gonzales et al. (2020) apply Brewster to a L + T subdwarf binary and provide evidence for their co-formation, as well as evidence for clouds in the primary. Peretti et al. (2019) use the retrieval code HELIOS-R on a combination of thermal infrared photometry and R ∼ 30 J-band data and place their retrieved chemical composition in context with both astrometric and radial velocity measurements. Molliere et al. (2020) employ the petitRADTRANS code to fit the near-infrared spectrum of the directly imaged planet HR 8799 e, finding an apparent degeneracy between solutions with significant cloud opacity and those with less cloudy atmospheres but with much shallower temperature gradients. This reflects a theoretical prediction from Tremblin et al. (2015 Tremblin et al. (, 2016) ) that the red J-H and J-K colors of many L dwarfs may just as readily be explained by a chemo-convective instability that produces vertical temperature gradients shallower than would be expected in thermochemical equilibrium. Nowak et al. (2020) produce and compare retrievals from both the ExoREM and petitRADTRANS codes on an R ∼ 500 K-band spectrum of β Pictoris b, finding excellent agreement between the retrieved C/O ratios of the two codes. Wang et al. (2022) ʼs retrieval on K-band data of HR 7672 B represents the highest-resolution spectrum used in an L-dwarf retrieval to date, at R ∼ 35,000, and they were able to precisely constrain the H 2 O and CO abundances, finding a C/O ratio consistent with the primary. Lueber et al. (2022) present a systematic retrieval of brown dwarfs across the L and T spectral types at an average resolution R ∼ 100 across the near-infrared, but they do not find any consistency or trends in the retrieved cloud properties for the L dwarfs. However, when considering mid-infrared (here 5-14 μm) spectra at similar resolution, Suárez & Metchev (2022) find that silicate features emerge starting at a spectral type of approximately L2, continuing through the mid-Ls, with the variability of the brown dwarf correlating positively with the presence and strength of silicate absorption. In this work we will present our own retrieval efforts on a widely separated companion classified as an early L dwarf and will use an additional retrieval on a previously studied field L dwarf to compare the results of our code with those of a different retrieval code on an object in a similar spectral class. 

We discuss the available data on the HD 106906 system and its companion in Section 2 and describe the components of our atmospheric forward model and retrieval code in Section 3. We test the ability of our code to converge on consistent results by using synthetic data in Section 4, benchmark the code by modeling a field L dwarf that has been previously retrieved on with a different code (Section 5), and finally show our results for the L-dwarf companion HD 106906 b in Section 6. Finally, we discuss the interpretation and limitations in our retrieval in Section 7 and summarize our key findings in Section 8.
section_name: Introduction

section_name: The HD 106906 System

section_content: The HD 106906 system consists of a pair of nearly identical mass young F-type stars (combined mass 2. 6 M e ) orbiting each other at 0. 36 ± 0. 002 au (Lagrange et al. 2016; Rodet et al. 2017; De Rosa & Kalas 2019). Its membership in the Lower Centaurus Crux (LCC) association (Gagné et al. 2018) places the system's age at 15 ± 3 Myr (Pecaut et al. 2012; Pecaut & Mamajek 2016). A companion, HD 106906 b, has an estimated mass of 11 ± 2 M J from fits to evolutionary models and sits at a projected angular separation of 7 11 ± 0 03 (Bailey et al. 2014). At a distance of 103. 3 ± 0. 4 pc (Brown et al. 2021), this places HD 106906 b at a projected physical separation of 734 ± 4 au (Zhou et al. 2020a). Table 1 provides a compilation of published system properties. 

Given its mass and orbit, HD 106906 b straddles the line between planets and brown dwarfs. Its exact formation pathway remains uncertain, as its mass ratio relative to-and remarkably wide separation from-its binary hosts poses challenges for all possible scenarios. To date, no studies of the HD 106906 system have provided substantial evidence for one formation pathway over another for the companion. On the planet formation side, there are several efforts to understand the dynamical history of the HD 106906 system, given the misalignment of the companion with the observed debris disk (Bailey et al. 2014; Kalas et al. 2015; Lagrange et al. 2016; Bryan et al. 2021). A core accretion pathway would require HD 106906 b to have formed interior to 100 au, which then requires a mechanism to evolve its orbit to the current projected separation in excess of 700 au. Wu et al. (2016) highlight the possibility that HD 106906 b could have been a planet scattered outward by its binary host, though the binary-planet scattering timescale is thought to be longer than the age of the system (Jílková & Zwart 2015). This hypothesis has been tested with efforts to constrain its orbital motion (Rodet et al. 2019; Nguyen et al. 2020). Nguyen et al. (2020) posit that HD 106906 b's orbit could have been excited in both orbital eccentricity and inclination from an unstable resonance with the binary. 12 However, this explanation is unlikely for HD 106906 in particular given the low density (<0. 11 stars per cubic parsec) of the LCC, which makes it unlikely that the companion's current position is the result of flybys scattering an initially closer-in orbit. The most recent study of the dynamical origin of this system as of the writing of this article is Moore et al. (2023), who provide an argument via numerical simulations that HD 106906 b could have been captured into the system as a planetary-mass free-floating object. They estimate that the probability of this scenario occurring within the past 5 Myr is ∼10 -6, which, while still low, is an order of magnitude more likely than in previous estimates. Swastik et al. (2021) demonstrate that the occurrence rate of companions at ∼10-1000 au shows a negative correlation with host metallicity-as opposed to the positive correlation seen in close-in gas giants-for masses greater than about 4 M J. This suggests that the formation histories of both the most massive planets and brown dwarfs may be dominated by gravitational instability, as the theory of formation by instabilities in the disk predicts a negative correlation with host metallicity (see, e. g., Helled & Schubert 2009). Bryan et al. (2021) find that the spin axis of HD 106906 b, its orbital plane, and the plane of HD 106906ʼs circumstellar disk are all mutually misaligned. They conclude that formation via gravitational instability is a plausible mechanism, as it is most consistent with misalignment across all three vectors. This scenario points to a C/O ratio consistent with the hosts, as this could occur either with gravito-turbulent instability or with fragmentation of a selfgravitating turbulent cloud.
section_name: System Properties

section_content: Photometries for the hosts and companion span the optical (F606W from Kalas et al. 2015 and ¢ z from Wu et al. 2016) through the thermal infrared ( ¢ L ; see Bailey et al. 2014). Two sources of photometry exist within the JHK s wavelength range; see Table 2. The first is from Bailey et al. (2014), who published J and K s magnitudes from the Magellan Adaptive Optics (MagAO) Clio2 instrument. The second is from Zhou et al. (2020a), who observed the HD 106906 system in the F127M, F139M, and F153M bands of the Wide Field Camera 3 (WFC3) of the Hubble Space Telescope (HST). The F127M bandpass overlaps with the J bandpass, as does F153M's bandpass with that of H, thus providing an independent (though not precisely congruent) comparison with our estimated H magnitude. F139M's bandpass falls almost entirely within the gap between the J-and H-band data. 

The highest-resolution spectrum of HD 106906 b comes from Daemgen et al. (2017), who present data obtained with the SINFONI integral field spectrograph on the Very Large Telescope (VLT). The data consist of three spectra in J, H, and K s, discontiguous with each other, with resolutions ≈2000-4000, at a signal-to-noise ratio (S/N) ∼20-50 (see Table 2 ). They derive an effective temperature of 1820 ± 240 K and a spectral type of L1. 5 ± 1. 0 based on comparisons with classifications in Allers & Liu (2013) and Bonnefoy et al. (2014). Daemgen et al. (2017) also classify the gravitysensitive features as most consistent with a very low gravity (consistent with the "γ" class, as defined and used in, e. g., Kirkpatrick 2005; Cruz et al. 2009; Allers & Liu 2013; Faherty et al. 2016). There are a few factors that affect the uncertainties in the reduced spectra. First, because the observations in Daemgen et al. (2017) lack a reliable H-band magnitude, we must make an estimate for the H magnitude. We choose to calculate J -H and H -K colors from a selection of lowgravity L dwarfs, published in Table 3 of Faherty et al. (2013). From these we take a weighted average with HD 106906 b's known J and K s magnitudes (Bailey et al. 2014) to obtain an estimate H = 16. 2 ± 0. 2. Second, Daemgen et al. 

(2017) identify regions, mostly at the edges of each spectroscopic band, that suffer large overall telluric absorption, as well as isolated wavelength ranges within each band (though concentrated in the H band) that may suffer from systematic uncertainties from the removal of telluric hydrogen in the data reduction. The reduced JHK s spectrum is normalized to the flux density at a specific reference wavelength in each band. 

A striking disagreement arises when comparing the flux densities inferred from the J versus that of the F127M photometry: the flux density derived from HST photometry is roughly twice as bright. To calculate this, we take the portion of our spectrum within the F127M filter and calculate how bright this object would be given the J magnitude, since the vast majority of the F127M band lies within the J band. From this calculation (done using the pysynphot package; see STScI Development Team 2013) we expect to see an F127M flux density of ≈3. 15 × 10 -13 erg cm -2 s -1 μm -1, versus the ≈6. 28 × 10 -13 that actually was observed with HST and reported in Zhou et al. (2020a). L dwarfs are known to be variable from photometric monitoring, with the notable case of VHS 1256 b with a ∼20% variability across 1. 1-1. 7 μm with a period of ≈21-24 hr (Bowler et al. 2020 ) and variability at the ∼6% level when extending to 5 μm (Zhou et al. 2020b ). The prevailing hypothesis for this variability is rotation (see, e. g., Zhou et al. 2016), with nonuniform cloud cover imparting brightness variations as changes in visible cloud opacity. However, Zhou et al. (2020a) found that HD 106906 b was only variable at ∼1% in the HST WFC3 F127M band, which is far smaller than would be needed to explain the discrepancy. 

One may prefer to adopt the HST photometry as the point of reference for the J and H bands, as space-based photometry provides a greater instrumental precision and does not suffer from systematics from telluric subtraction. However, there are no direct HST comparisons for the K s band, and so one will still need to mix the sources of photometry to calibrate the entire spectrum. We do not resolve the disagreement between these two sources of photometry-instead, we run our retrieval of HD 106906 b with the original J and K s photometric normalization. We use calibration parameters in flux normalization to attempt to capture uncertainties in this flux normalization. These uncertainties stem from the error bars on the J and K s magnitudes used to normalize the spectra in Daemgen et al. (2017). 

We do not have a direct constraint on the C/O ratio of the hosts; such a constraint is needed to compare with the C/O ratio we retrieve to test the hypothesis that the HD 106906 b formed as a substellar companion. One approach is to combine a C/O-to-[Fe/H] relation for planet hosts with metallicity measurements for stars in the galaxy. In the case of HD 106906, one can use metallicities from members of the Upper Centaurus-Lupus (UCL) and LCC associations (see Table 1 erg s -1 cm -2 μm -1 ) 4. 39 ± 0. 04
section_name: Data

section_content: a Estimated from J -H and H -K s colors for low-gravity L dwarfs; see Table
section_name: Notes.

section_content: To model the atmosphere of HD 106906 b from its emission spectrum, we employ the APOLLO code (Howe et al. 2017 (Howe et al., 2022)), a model framework for generating spectra of planets in both transit and emission. 13 The core of the forward model is modeling the combination of thermal emission, absorption from atomic and molecular species, and extinction (scattering and absorption) from clouds. APOLLO uses a hemispherical approximation to the Toon et al. (1989) twostream scattering routine, which is used primarily for cloud scattering. To calculate the emission spectrum, the outgoing radiation is averaged over eight angle divisions in the outgoing hemisphere. The code is designed to be modular in parameterizations of molecular abundances, temperature-pressure, observing modes, and noise models, with particular focus on observing configurations for JWST. APOLLO is equipped with a likelihood sampling routine that serves as the retrieval component of our model (discussed in Section 3. 4). The parameters used in our retrievals, including the bounds on their priors for the parameter estimation routine, are listed in Table 3.
section_name: Atmospheric Modeling and Spectral Retrieval

section_content: The molecular radiative transfer scheme in APOLLO relies on sampling of cross-section tables for a variety of applicable species; these are precomputed from a grid of line-by-line opacities and are derived from the sources in Table 1 of Freedman et al. (2014), with the exception of the alkalis. The opacities for the alkali lines are drawn from the Lupu et al. (2022) catalog, which derives its Na and K profiles from a series of works analyzing the interaction between atomic lines and molecular hydrogen (Allard et al. 2007 (Allard et al., 2012 (Allard et al., 2016 (Allard et al., 2019)). We employ two levels of down-sampling to create the cross sections we use in our forward models, with resolutions of 10,000 and 50,000. For retrievals on real data, we choose the minimum opacity resolution that ensures that the ratio between the mean opacity and data resolution is 100, following a community recommendation to avoid introducing artificial errors from binning effects. 14 We freely retrieve fractional abundances for H 2 O, CO, CO 2, H 2 S, CrH, FeH, TiO, and VO. We assume a solar H/He ratio, and H 2 and He opacities include collisionally induced absorption. Atomic Na and K are included together as a single free parameter, where the ratio of their abundances is fixed to that of solar metallicity (see, e. g., Line et al. 2015). For the molecular abundances we assume a constant mixing ratio and initialize at values corresponding to the chemical equilibrium abundances at the pressure layer closest to the literature effective temperature (here taken to be 1820 K from Daemgen et al. 2017 ). The equilibrium abundances were calculated through a routine in the PICASO atmospheric radiative transfer code (Batalha et al. 2019). 

We visualize the contributions of the gas to the emission spectrum by calculating a contribution function per atmospheric layer, which is given by ( ) where the atmospheric layer spans pressures P to P + ΔP, ( ) T P is the temperature in the layer, and B λ is the Planck function at that temperature. τ sp and τ tot represent the optical depths due to a given gas species and from the entire contents of the layer, respectively. The contribution function is expressed as fractions of the total across an entire vertical column in the atmosphere. The function, when summed across all gas and cloud species, is proportional to the pressure derivative of the "transmittance"

) times the Planck function at the given pressure and temperature; see, for example, Section 3 of Line et al. (2014).
section_name: Molecular Species and Opacities

section_content: Our T-P profile is adapted from the parameterization proposed in Section 4. 2 of Piette & Madhusudhan (2020, see their Figure 8 ), with additional temperature nodes added to the extremes of the profile. The parameterization is designed to be flexible enough to accommodate a wide range of possible vertical thermal structures, including approximating a radiative-convective equilibrium, while also filtering out excessive unphysical behavior, such as the "ringing" that was described in Line et al. (2015 Line et al. (, 2017)). The parameters are the temperatures 2. 5-7. 5

Gases (log 10 number abundance relative to total; see Section 3. 1)

Temperature at 10 0. 5 bars (T 0. 5, K) 75-4000 Single-scattering albedo ( )

Notes. All priors are uniform within the listed bounds. The calibration factors are also only used for retrieval on HD 106906 b. CO 2 is not included as an absorber in the retrieval of the spectrum of 2M 2224. 

a The temperature profile is also constrained to be monotonic; see the discussion on the dependent priors in Section 3. 2. b An additional constraint is imposed such that the cloud layer does not extend beyond the base of the model. 

at 10 pressure levels, representing nodes between which we interpolate the profile. The temperature nodes are spaced in orders of magnitude from the top of the model atmosphere (10 -4 bars) down to a pressure of 1 bar, beyond which we use half-orders until we reach the deepest pressure of the model at 10 2. 5 bars. Here we label the temperatures of each node by subscripts denoting the base-10 logarithm of their corresponding pressure. We follow the recommendation of Piette & Madhusudhan (2020) to use a monotonic spline interpolation with a Gaussian smoothing kernel of width 0. 3 dex in logpressure, as the mechanism by which one can filter out the aforementioned ringing. In the original setup, the temperature at a pressure of 10 0. 5 bars (T 0. 5 ) is taken as a reference temperature at a fiducial pressure approximating the depth of the photosphere for a typical self-luminous brown dwarf. In this setup, the remaining parameters then define the differences in temperature between each successive node. In contrast, we choose to define all our parameters as the temperatures themselves, but we use an iterative process for proposing temperatures by determining the bounds on the uniform priors for each temperature:

1. The bounds of the prior for the photospheric node (T 0. 5 ) are set by the bounds of the temperatures of the opacity tables (75-4000 K). 2. Then, the shallowest (T -4 ) and deepest (T 2. 5 ) temperature prior bounds are each bounded by the proposal for T 0. 5 and by the minimum and maximum opacity temperatures, respectively. 3. This continues with the nodes closest to the middle of the existing nodes being bounded by those already-chosen nodes, subdividing until the whole profile is bounded and all temperatures proposed. 

This ensures that the profile is monotonic in temperature.
section_name: Temperature-Pressure Profile

section_content: Our cloud model is modeled after the "slab" approaches used in Burningham et al. (2017) and Gonzales et al. (2020). The model cloud occupies a fixed region in pressure space, with a minimum pressure where cloud absorption begins (the cloud "top") and some depth in pressure. The vertical opacity profile is restricted to follow ∂τ/∂P ∝ P. The free parameters include the pressure of the cloud top P top, the depth of the cloud in logpressure space ( ) ( ) D º P P P log log 10 cloud 10 base top, and a wavelength-dependent opacity and single-scattering albedo instead of particle-specific parameters. The wavelength dependence is modeled as a power law with exponent α. The opacity at a given pressure depth and wavelength is therefore given as

where τ 0 is the maximum optical depth (at the base of the cloud at a pressure P base ) at a wavelength of 1 μm. This is an empirical approximation to scattering by cloud particles whose sizes are smaller than the wavelengths of observation. Previous efforts at retrievals in this wavelength range indicate that a model that takes into account specific condensates for its opacity calculations is not preferred over the simpler approach used in this work. 

The final free parameter for our cloud model is the singlescattering albedo ω 0. This is the ratio of photons that are scattered versus those extincted overall (either absorbed or scattered). By choosing to model the albedo with a single free parameter, we assume that it is constant across all wavelengths and pressures. To be precise, ω 0 in this work refers to the single-scattering albedo of the clouds alone; in APOLLOʼs implementation of the Toon et al. (1989) radiative transfer model, their ω 0 refers more broadly to the scattering-toextinction ratio of all absorbers and scatterers in the atmosphere. For our purposes this means including the gas as well, for which we model scattering as Rayleigh scattering since the sizes of each molecule are much smaller than the observed wavelengths.
section_name: Cloud Models

section_content: We sample likelihoods in parameter space with a nested sampling algorithm, using the dynesty Python package (Speagle 2019). We choose to set uniform priors on all parameters, the ranges of which are listed in Table 3. Each model was initialized with 1000 live points in the "rwalk" sampling method. Models were run with the built-in default stopping criterion for assessing convergence, which depends on the amount of evidence accounted for in the cumulative samples. 15 The total number of effective iterations in each run varies based on when the stopping criteria are reached, with test retrievals on simulated data (Section 4) using ∼10 5 and retrievals on real data (Sections 5 and 6) requiring 2-3 times as many. 

Once the runs are complete, we then derive the mass, effective temperature, metallicity, and C/O ratio. The mass is calculated directly from the radius and surface gravity. The effective temperature is calculated from an approximation to the bolometric luminosity, using a low-resolution (R ≈ 200) spectrum that covers 0. 6-30 μm. 16 We report metallicity by comparing the mean molecular weight versus that expected for solar metallicity, rather than reporting a metallicity as an [Fe/ H] value. For the mass fraction Z of non-H/He elements, the metallicity is calculated as

, where we take Z e = 0. 0196. We choose this definition for metallicity because our values of metallicity are not tied to a specific atomic species, and the way in which we model the abundancesuniform in pressure but freely variable-means that our model does not require the abundances to be in chemical equilibrium. 

We use Bayes factors to compare the quality of fits to data between two models. The Bayes factor is simply the ratio of the marginal likelihoods (also known as evidences) of each model's retrieval. A higher Bayes factor confers stronger support for a model relative to another; a recommendation originally proposed in Jeffreys (1998) is to interpret a ratio of 10-10 1. 5 as "strong," 10 1. 5 -10 2 as "very strong," and >10 2 as "decisive" confidence that the model with the higher evidence is preferred. Following this, Benneke & Seager (2013) adapted the heuristics in Table 1 of Trotta (2008) to translate the language of Bayes factor comparisons into a "detection significance," usually quoted in units of σ. This is a convenient way to express analogous statistics in both the Bayesian and frequentist frameworks of model analysis; we report both in the following sections for our model comparisons.
section_name: Parameter Estimation Methods

section_content: This work represents the first application of APOLLO to data from an L dwarf. To test the efficacy of the code, we generate a forward model that approximates the object, making simple assumptions about the atmospheric structure. In principle, our retrieval code should be able to converge on a good fit (i. e., reduced chi-square statistic c ñ 1 2 ) to a data set generated from its own forward model and, based on the distributions of the retrieved parameters, should inform us to how well each parameter could be constrained from a near-infrared wavelength range and S/N similar to that of HD 106906 b. 

We use APOLLO to generate a forward model spectrum for a 1. 5 R J, = g log 4. 19 object; the choice of radius is arbitrary, but the surface gravity is taken from the best estimate of HD 106906 b's gravity from the observed luminosity and effective temperature based on the fits made in Daemgen et al. (2017). For the thermal profile, we produce a parameterization that approximates a SONORA profile at approximately 1800 K if no cloud opacity is present. We use PICASO, specifically the Visscher chemical equilibrium code (Marley et al. 2021), to generate equilibrium abundances for the model pressures given the above parameters. This corresponds to a C/O ratio of 0. 54 and a metallicity of 0. 065. We generate data for two cases: one with clouds, parameterized as described in Section 3. 3, and one "clear" case without clouds. For the clouds, we use a layer that spans ∼10 -0. 5 to 10 1 bars, chosen to bound the estimated photospheric pressures, and has enough opacity to yield an effective temperature of ≈1360 K. All models used to generate these data are identical in all noncloud parameter values. Noise is modeled as independent, Gaussian (white) noise at S/ N = 20, approximately the minimum S/N seen in the HD 106906 b spectrum. Comparisons of the forward model spectral fits, T-P profiles, and distributions of model parameters are shown in Figures 1 2 3 4 5 (on cloud-free data) and Figures 6 7 8 9 10 (on cloudy data).
section_name: Retrieving on Simulated Data from Forward Models of Low-gravity L Dwarfs

section_content: The retrieval on simulated cloud-free data has no issue converging to an excellent fit, with the final c n 2 value very close to 1 (Figure 1 and Table 4 ). The retrieved C/O ratio is consistent with the input value of solar (0. 54), with a 68% confidence interval in the posterior distribution of ±0. 003. The only abundance not tightly constrained is CO 2, with the true abundance sitting above the upper limit of the 68% confidence interval. The retrieved T-P profiles show a weak constraint at either end of the pressure range, with pressures smaller than ∼10 -3 bars or larger than ∼10 bars. The best-fit T-P of the cloudy model is closer to the true profile at the shallowest pressures, but the confidence ranges of the two models overlap significantly at these pressures, meaning that the relative fit Figure 1. Spectra and contribution functions of the retrieved forward model fits to data simulated from cloud-free forward models of an L dwarf. The forward model spectra (using the MLE parameter values) are in color, with the data in gray. The retrieved models are nearly identical and overlap nearly entirely. Immediately beneath the spectra are the contribution functions for the two principal carbon-and oxygen-bearing species; CO 2 is included in the model used to generate the simulated data, but its contributions to the emission are well below those of H 2 O and CO. The deepest contours, outlined in solid colors, enclose the regions where the contribution function reaches >1% of the total contribution within the atmospheric column at a given wavelength bin. Each successive contour denotes 2 orders of magnitude smaller fractional contribution (here, 10 -4 and 10 -6 ). The faint gray lines in the J-band (leftmost) sections contribution plots denote the location of the cloud layer as retrieved by the cloudy model on the cloud-free data. The faintness of the lines denotes the low optical depth of the cloud layer, in contrast with the darker cloud contours as seen in the retrieval on cloudy simulated data (Figure 6 ). 

qualities in this region of the atmosphere are not significantly different. The contribution plots show that most of the contribution is concentrated between pressures of a few tenths of a bar to a few bars; therefore, it is not surprising that most of the uncertainty in the thermal profile arises away from these intermediate pressures. 

The nominal expectation is that the cloudy model, when applied to cloud-free data, will effectively "turn off" the cloud opacity. This is largely true; the median opacity at the reference wavelength of 1 μm is very weak (an optical depth of ∼0. 01, yielding an attenuation of at most a few tenths of a percent). However, there is a tail in the distribution of optical depths; in some cases the model will choose nonnegligible cloud opacity, but the increases in optical depth correlate with the depth of the cloud top. This is consistent with the relative lack of contribution to the emission spectra from pressures deeper than a few bars, which is where the distribution of optical depths reaches ∼1. The single-scattering albedo ω 0 is low, particularly in the low-opacity cases, and is weighted toward low power-law exponents, which would allow nonnegligible cloud opacity at wavelengths <1 μm. This may explain why the distribution of effective temperatures for the cloudy model fit peaks near the true value of 1821 K but has a substantial secondary peak, with the median at 1274 K. If one were to extend the forward models to shorter wavelengths, we would see these models diverge from their cloud-free counterparts. 

Regardless of the precise nature of the way in which the cloud model withholds its opacity from the spectrum, the retrieved distributions of the gas species are very similar, and the constraints on the C/O ratio are of nearly identical accuracy and precision. The cloud-free model returns a Bayes factor higher than its cloudy counterpart by a factor of approximately 18, meriting its preference following the interpretation of Jeffreys (1998). In the interpretation of Benneke & Seager (2013), we could say that we "detect" the cloud-free model with a significance slightly less than 3σ. The cloud-free model gains its advantage by virtue of using fewer free parameters to fit the data.
section_name: Retrievals on Cloud-free Data

section_content: We now show the fits to mock data with clouds included, with the retrieved spectra and contributions in Figure 6, T-P profiles in Figure 7, and posterior distributions of parameters in Figures 8-10. The power-law cloud opacity model yields a reduced chi-square statistic c » n 1 2. The retrieved C/O ratio is slightly less accurate when compared with that of the cloudfree case, with the true value lying just outside the 68% confidence interval (but well within the 95% interval). 

The cloud-free model applied to cloudy data returns a slightly worse fit, with the reduced chi-square statistic increasing to 1. 07. In this case the model compensates for a lack of clouds by decreasing the radius, increasing the gravity, and increasing the abundances of all species except that of CO 2 by 0. 05-0. 1 dex. This allows for a C/O ratio distribution that is still marginally consistent with the truth at the 95% confidence level. When comparing the quality of the fits between the cloudy and cloud-free models, we find a Bayes factor of approximately 150. This puts the modest difference in the reduced chi-square statistics in greater perspective; the cloudfree model is strongly disfavored when compared with the model with clouds, with a frequentist translation to a model that is preferred at ∼3. 6σ. This is due to an accumulation of minor differences ( the typical error bar) between the cloudy and clear fits, primarily in the J band, where the simulated cloud opacity is strongest. While the cloudy model better fits the data, it does not reproduce all cloud properties precisely. The cloud top position and extent of the cloud are tightly correlated (see Figure 10 ), and the reference optical depth for the maximum likelihood estimator sits at the high tail of its posterior distribution, with the true value even farther out. The true value for the pressure of the cloud top sits at the high end of the 68% confidence interval for the retrieved posterior distribution, which shows that the model is able to reproduce where in the atmosphere cloud opacity should become significant. Since the reference optical depth is inaccurate, this suggests that there is some minimum opacity that effectively suppresses much of the emission deeper than the cloud; no additional opacity is needed; therefore, the model finds a solution centered around the minimum sufficient opacity. 

Finally, when comparing the retrieved T-P profiles, we see that, as in the test retrievals on spectra generated from a model without clouds, the weakest constraints in temperature arise in the shallowest and deepest parts of the atmosphere. However, we now see an additional divergence between the cloud-free and cloudy fits-namely, the cloud-free profile diverges from the true profile within the simulated cloud layer, while the cloudy profile remains close to the true profile. Then, by the time we reach the deepest extent of the clouds, both profiles 4. 

have begun to diverge from the true profile. This suggests that the cloud-free model is attempting to compensate for its lack of clouds by keeping the temperature gradient shallower, suppressing its thermal contribution in a way that can mimic the effect of a cloud layer.
section_name: Retrievals on Cloudy Data

section_content: Taking the results of self-retrievals on cloudy and clear simulated data together, we demonstrate that our code is able to identify both the correct abundances and the thermal structure in the photosphere in an atmospheric simulation. The results also highlight what we might not expect to constrain precisely owing to theoretical limitations, such as the deepest parts of the T-P profile and the opacity profile of the clouds. Additionally, a cloud-free model may be able to reproduce a C/O ratio consistent with the true value but risks returning an inaccurate radius and gravity and molecular abundances that are almost all consistently too high, all with a temperature profile that is consistent at the same pressure ranges as the cloudy model but with higher uncertainties at the lowest pressures. Additionally, when clouds are present, we expect a cloud-free model to show 4. 

the greatest difference from that of a cloudy model within the cloud layer itself, changing its gradient to compensate for the lack of extinction from condensates. The bias in molecular abundances suggests that, at least in this wavelength range, the shape of the spectrum is determined more by the relative abundances than by the absolute abundances; put another way, we may expect to see a potential degeneracy between the T-P profile, gravity, and key molecular abundances but nevertheless may expect the retrieved C/O ratio to not be significantly biased away from the true value. However, these conclusions are necessarily limited to which physics we choose to include in the model used to simulate the data; we are limited to commenting on the efficacy of the code in terms of the consistency of retrievals with the assumptions we have made.
section_name: Lessons from Self-retrievals

section_content: Our first true retrieval is of the mid-L field dwarf 2MASSW J2224438-015852 (Kirkpatrick et al. 2000), which we refer to here as 2M 2224. This is one of the brown dwarfs studied with the Brewster retrieval code in Burningham et al. (2017) and with mid-infrared data in Burningham et al. (2021). To 4. As there is no cloud opacity in the simulated data, the model can adapt its cloud opacity in several ways to effectively remove its influence on the resulting emission spectrum. The first is to turn the reference optical depth ( ) t l 0 to a very low value (=1); the second is to introduce significant opacity ( ) t l 1 0  but to place the cloud deep into the atmosphere, below where the majority of the thermal emission originates in the spectrum (i. e., below the photosphere). benchmark our code against previously published results in the JHK s spectral range, we limit our reanalysis to the original, R ∼ 75 spectrum from Burgasser et al. (2010). 17 The full wavelength range of the data is 0. 65-2. 56 μm, but to compare with Burningham et al. (2017), we choose to only use the range from 0. 8 to 2. 4 μm, though a retrieval was performed with the full data set. Burningham et al. (2017) ) as a free parameter in our corner plots. It should be noted, however, that Burningham et al. (2017) report their results using a "deck" cloud model, though both that and the slab model were tested in their work. Finally, our sampling method differs in that Burningham et al. (2017) used an MCMC parameter estimation technique and imposed more restrictive priors on gravity to keep their mass below the nominal main-sequence limit of 80 M J. Our choices of a more free T-P profile parameterization and wider priors on the gravity mean that our code can explore a broader range of solutions for the vertical atmospheric structure, but with a concession that our solution has a higher risk of introducing structure to the profile that does not have a feasible physical interpretation. 

Results from the retrieval are shown in Figure 11 for the spectrum and contributions and in Figure 12 for the retrieved T-P profiles, and posterior distributions for parameters are shown in Figures 13-15 and tabulated in Table 5. Our retrieved spectrum does not precisely reproduce the shapes of the local peaks in the J and H bands and generally prefers a "smoother" (though not necessarily better) fit to the spectrum than that retrieved in Figure 8 of Burningham et al. (2017). Our retrieved C/O ratio of -+ 0. 86 0. 02 0. 01 sits entirely within the confidence interval reported in Burningham et al. (2017), despite retrieving a higher gravity and higher abundances, particularly in H 2 O and CO. Our model finds a solution that prefers a higher metallicity (1. 61 ± 0. 14). Our T-P profile mimics the shape of the 2017 paper from ∼0. 01 bars until the location of our retrieved cloud layer, where the models then diverge. The T-P profiles diverge most strongly where the preferred deck cloud model of Burningham et al. (2017) ). The extent of our cloud layer encompasses their median τ = 1 pressure, but our retrieved Figure 6. Spectra and contribution functions of the retrieved forward model fits to data simulated from cloudy forward models of an L dwarf. The forward model spectra (using the MLE parameter values) are in color, with the data in gray. The retrieved models are nearly identical and overlap nearly entirely. Immediately beneath the spectra are the contribution functions for the three principal carbon-and oxygen-bearing species. The deepest contours, outlined in solid colors, enclose the regions where the contribution function reaches >1% of the total contribution within the atmospheric column at a given wavelength bin. Each successive contour denotes 2 orders of magnitude smaller fractional contribution (here, 10 -4 and 10 -6 ). 

power-law dependence in wavelength is much more steeply negative than theirs, with a median a = - -+ 7. 73 0. 78 0. 71. In other words, our model prefers a solution that allows significant cloud opacity in the J-band portion of the spectrum but rapidly diminishes at longer wavelengths. Our model has difficulty finding a solution to the atmosphere from ≈0. 8 to 1. 3 μm, given that the model fit fails to capture the smaller-scale variations in the data in the region where it determines that clouds are most significant. This may mean that the cloud model is instead being used to compensate for an inability to fit this portion of the spectrum, while fitting the remainder in the H-and K s -band ranges more accurately. An earlier retrieval with the full 0. 65-2. 56 μm data set did return a cloud powerlaw exponent of α = -2. 04 ± 0. 02 but also had its own difficulties in capturing the entire spectrum, with a better fit to the 0. 8-1. 2 μm region but a worse fit in the K s band from 2. 00 to 2. 35 μm and a similarly high gravity. 

This agreement in C/O despite disagreement elsewhere is similar to the findings in works such as Molliere et al. (2020), where their tests of models with different cloud models yielded similar C/O ratios despite retrieving disagreeing thermal and cloud profiles. In our case, with the more flexible thermal structure, there is an additional degeneracy between the gravity and the molecular abundances/metallicity. The higher the metallicity, the less deep in the atmosphere a given optical depth will be reached, but the higher the gravity, the smaller a path length for a given change in pressure, meaning that the equivalent optical depth will occur at a higher pressure. Our choice of T-P profile allows flexibility in adapting the shape of the vertical thermal profile to changes in model gravity and metallicity; therefore, we expect gravity and metallicity to be negatively correlated. Mirroring the behavior we saw in the cloud-free versus cloudy models applied to simulated data in Section 4, it is possible to retrieve an accurate C/O ratio by retrieving abundances that are accurate relative to each other but biased in their absolute values. It is difficult to compare the shallow thermal gradient with the behavior suggested in Tremblin et al. (2016), where a shallow temperature gradient driven by a thermochemical instability can mimic some of the spectral behavior attributed to clouds, since in this case both the shallow gradient and significant cloud opacity are present in the model solution. Nevertheless, we keep these findings in mind when interpreting the results of our retrieval on HD 106906 b (Section 6).
section_name: Retrieval on a Previously Characterized L Dwarf

section_name: Retrieved Atmospheric Properties for HD 106906 b

section_content: When moving from the test retrievals on simulated data to retrievals on the actual HD 106906 b data, there are a few differences in the model setup, though the core physical model remains the same. The first is the addition of calibration terms that scale the flux in each band by a multiplicative constant; this is to account for uncertainties in the photometry, as discussed in Section 2. 2. These calibration scales are partly degenerate with the retrieved radius, so when reporting these calibration scales, we normalize the radius such that in each case the effective calibration scale in the K s band is 1. The second is to add a parameter for fractional cloud cover ( f cloud ). Since we are using a 1D (vertical-only) model, the fractional cloud cover is assumed to be isotropic, and the emission flux is simply weighted between the fully cloudy flux one calculates from the given parameters (F cloudy ) and the flux given the same parameters but without clouds (F clear ): The final modification is that the data are down-sampled to a maximum resolution of ≈500, a factor of 4-8 lower than the original. This is because the maximum resolution of our opacity tables is 50,000, and to avoid introducing excess artificial noise from binning effects, we impose a limit of R opacity /R data 100, which requires us to resample the data to a lower resolution. We calculate the uncertainties in the downsampled data as uncertainties in the mean, i. e., since we now have N = 4-8 resolution elements of the original spectrum in each of the down-sampled elements, our uncertainties in each new element are assumed to be smaller by a factor of -N 1. This is a lower limit of the true uncertainties in the new spectrum, as the errors between the original pixels and therefore resolution elements are almost certainly correlated at some level. To make an estimate of the typical correlation length, we use the approach in Section 2. 2 of Line et al. (2015), where one calculates the autocorrelation of the residuals for an initial model fit to the data. Doing this, we find that the autocorrelation drops and subsequently remains at or below ≈0. 25 at a scale of 6-8 pixels. Therefore, our reported fit qualities, such as chi-squared statistics, may be overestimated by roughly a factor of 2-3. However, a larger source of systematic errors comes in the form of telluric contamination, which is typically strongest at the boundaries of each band. 4. 

These errors are more likely to introduce biases in the retrieved atmospheric parameters, and so we perform retrievals to examine the effects of including versus excluding these portions of the spectrum. 

An initial retrieval was run on the full HD 106906 b data set, along with retrievals on data from each band (J, H, and K s ) individually; see Figure 16. The purpose of this initial set of runs was to understand how well the retrieval could fit the data and, depending on whether and how the individual band fits compare with the full-spectrum fit, to determine whether the model is capturing the physics of the companion's atmosphere with adequate flexibility (e. g., this comparison can provide a first-order check to whether the assumption of constant mixing ratios is sufficient). The result is broadly that there are wavelength regions of each band where neither the single band nor the full spectrum fits the data well. The largest discrepancy occurs at the blue end of each band, as well as to a lesser extent at the red end of the H band. These regions are consistent with the wavelength regions identified in the original publication of the data (Daemgen et al. 2017, specifically Figure 2 ), where telluric contamination is thought to affect the data reduction most severely. Given the overlap of the most poorly fit regions with the suspected high-telluric regions, we choose to excise these data from the final retrievals. The final 4. 

fit, using the truncated data across all bands, is overplotted in Figure 16.
section_name: Retrieval Setup: Single-band Trials and Regions of High Tellurics

section_content: Results from the retrieval with our cloud model included are shown in Figure 17 for the spectra and contribution functions and in Figure 18 for the retrieved T-P profiles, and posterior distributions of parameters are shown in Figures 19 20 21 and tabulated in Table 6. The retrieved model captures much of the broad shape of the spectrum but fails to capture the amplitudes of some of the features in the J band. Given the S/N of the data, the c n 2 statistic of the model indicates a poor fit to the data,

. As noted in Section 6. 1, we have scaled the original uncertainties assuming that the errors are uncorrelated; this statistic assumes the most optimistic noise model and therefore the most pessimistic fit quality. The retrieved radius posterior range (1. 30 ± 0. 06 R J ) is smaller than the radius one derives from the best-fit bolometric luminosity and effective temperature in Daemgen et al. (2017), which is 1. 47 R J. However, this is affected by our choice to normalize the spectrum such that the K s -band calibration factor is 1; therefore, there is a bit of ambiguity in whether our retrieved radius is strictly consistent or inconsistent with these previous constraints. Our inferred effective temperature is high given the small radius. The retrieved surface gravity is low compared with that derived from the Daemgen et al. 4. 

confidence interval does overlap with this range. This low gravity, combined with the radius, yields a 68% confidence interval of -+

1. 92 0. 70 1. 48 M J for the mass, with a maximum likelihood estimate (MLE) value of 4. 41 M J. This is smaller than the mass range of 11 ± 2 M J from evolutionary models as presented in Bailey et al. (2014), as well as the estimated 13 M J mass if one were to adopt the mean age of the LCC, at 17 Myr. We are unlikely to be able to disentangle the low retrieved mass from the existing degeneracies that persist between gravity, metallicity, and the T-P profile. Additionally, the range of bolometric luminosities we infer from our results is low compared with the original evolutionary model constraints: our cloudy model returns a 68% confidence interval of Notes. The two models (cloud-free and cloudy) used to retrieve on the simulated data were the same as those used to generate the cloud-free and cloudy data. Both sets of parameter values are identical in all noncloud parameters, and the values used in generating the simulated data are also identical except for the inclusion of clouds. We show results from each possible pairing of clear and cloudy model versus clear and cloudy data. a 1821 K for the cloud-free data, and 1361 K for the data with power-law cloud opacity. 

the original luminosity estimates and those from our cloudy model is the fact that our effective temperature is computed by generating a low-resolution forward model over a longer wavelength range than the data, at 0. 6-30 μm. With such a large negative exponent to the wavelength dependence of our retrieved cloud model, the emission bluer than the data wavelength will be suppressed, cooling the derived effective temperature. However, our cloudy model is consistent with the luminosity ranges derived using subsets of "young" ("YNG" and "YNG2") targets, as presented in 0. 53 0. 25 0. 15 is consistent with the estimated C/O ratio distribution of the stellar association in which HD 106906 resides (0. 52 ± 0. 11; see Equation (2)); the three primary C + O constituents (H 2 O, CO, and CO 2 ) are constrained to within 0. 5-1 dex and show positive
section_name: The Cloudy Model

section_content: a Bulk metallicity and reference optical depth not reported. b We use a different T-P parameterization from that of Burningham et al. (2017). A plot comparing the two retrieved profiles is available in Figure 12. Immediately beneath the spectra are the contribution functions for the two principal carbon-and oxygen-bearing species; CO 2 is not included in the retrieval, as it was excluded from the retrieval in Burningham et al. (2017), inviting a more direct comparison of our results. The deepest contours, outlined in solid colors, enclose the regions where the contribution function reaches >1% of the total contribution within the atmospheric column at a given wavelength bin. Each successive contour denotes 2 orders of magnitude smaller fractional contribution (here, 10 -4 ). correlations among each other, as well as with the surface gravity. The correlation between molecular abundances and gravity is known to be a consequence of a degeneracy where, in fitting absorption features, the flattening effect of higher gravity can be at least partially offset by higher abundances (see, e. g., Todorov et al. 2016). The full posterior distributions for gas abundances are shown in Figure 20. The retrieved H 2 O abundance for the best-fit model (-3. 33 dex) is within 0. 1 dex of the expectation given the retrieved T-P profile, if one assumes chemical equilibrium for an object at solar metallicity and C/O ratio (-3. 35 dex). The CO abundance, at -3. 02 dex, is higher than the equilibrium value of -3. 28, which drives the best-fit C/O to just beyond the 68% confidence interval, at 0. 66. CO's impact is comparable to that of H 2 O but over isolated regions of the spectrum; the bulk metallicity has uncertainties of order 0. 4 dex but is consistent with solar metallicity, as well as the metallicity range of its stellar association. The CO 2 abundance (-

5. 21 1. 00 0. 56 dex, best-fit value -4. 38 dex) is the least constrained of the three major C + O molecular absorbers and has the smallest effect on the C/O ratio. The absorbers least consistent with an equilibrium abundance are the alkalis; with a range of - -+ 9. 38 1. 86 3. 15 and a best-fit abundance of -9. 97 dex, the model essentially ignores the alkali absorption features in its fit. This is surprising since there are prominent absorption lines of potassium (two K I doublets) in the J band. The failure of the model to capture 5. Median values for the equivalent parameters in the run in Burningham et al. (2017) are listed in Table 5. 

these absorption lines appears to be a consequence of reducing the spectral resolution; previous attempts at retrievals were less successful at converging to a global atmospheric fit than the ones shown in this work but had enough of the line shape at the original resolution to fit the abundances, as well as using the ratio of the K I doublet line depths as an additional constraint on gravity. Such models return alkali log abundances ranging from about -5 to -7 but also return infeasibly high gravities, exceeding = g log 6 unless a restrictive prior is used. The preferred T-P profile is shallow in its temperature gradient from the top of the atmosphere to a pressure of several bars, after which point the temperature rapidly increases to approximately 3100 K at several tens of bars. The profile then returns to a nearly isothermal behavior to the base of the model atmosphere. Figure 18 shows the profile along with its cloudfree counterpart and a cloud-free SONORA brown dwarf profile interpolated to match the maximum likelihood gravity and effective temperature from the cloudy model. In contrast with the radiative-convective equilibrium profile from SONORA, with a shallow thermal gradient gradually increasing to a higher adiabatic gradient at the radiative-convective boundary, our retrieved profile can be described as nearly isothermal layers for the log-pressure ranges of ∼ -3 to -2, again from ∼ -1 to 5. Median values for the equivalent parameters in the run in Burningham et al. (2017) are listed in Table 5. 0. 5, and, at least for the cloudy model, a nearly isothermal layer at the deepest ∼1 dex of the model pressure range. These nearly isothermal "layers" are punctuated with comparatively rapid temperature increases. The majority of the contribution from the major absorbers (H 2 O, CO, and CO 2 ) comes from pressures of ∼1 mbar-1 bar. The full posterior distributions for the cloud parameters are shown in Figure 21. The distribution of cloud top pressures ranges from the very top of the model (10 -4 bar) to a few millibars, and the top pressure is strongly correlated with the depth of the cloud. The maximum pressure of the cloud appears to be the most important parameter here, which, when combined with the fact that most of the gas contribution to the emission is beneath this maximum pressure, implies that the model prefers whichever cloud layer can produce some fixed total column optical depth. With a highly negative power-law exponent (a = - -+ 7. 15 0. 73 0. 66, best-fit value -7. 05), clouds produce significant opacity only for the J band. As with the retrieved clouds for 2M 2224, this might not reflect an accurate constraint on actual cloud opacity, but for the model it is a way to suppress emission in the bluest wavelengths without an obvious physical interpretation. The retrieved distribution of the single-scattering albedo ω 0 is tightly distributed and is close to the upper limit. The covering fraction ( f cloud ) distribution is consistent with but not centered at 1, which corresponds to a near-global coverage of a very reflective cloud. 5. Median values for the equivalent parameters in the run in Burningham et al. (2017) are listed in Table 5.
section_name: Notes.

section_content: The spectral fits, contributions, profiles, and posterior distributions are plotted with the cloudy solution in Figures 17 18 19 20 21, with parameter constraints also in Table 6. Excluding clouds entirely in our model returns a worse fit quality, with c = n 61 2. The log-ratio of the Bayes factors is 56, indicating that the cloudy solution provides overwhelmingly stronger evidence than the cloud-free solution. The main reduction in fit quality is in the J band, consistent with where the cloudy solution places most of its cloudy opacity. While the gravity ( = -+ g log 3. 43 0. 49 0. 35, best-fit value 3. 69) is consistent with the cloudy fit, the radius of the cloud-free solution increases to 1. 74 ± 0. 06 R J, and the abundances of all species except H 2 S increase by amounts ranging from 0. 5 to 2 dex. The degeneracy between gravity and the molecular abundances is even stronger than that of the cloudy atmospheres, with the CO abundance rising more than the H 2 O abundance, yielding a C/O distribution of -+ 0. 87 0. 04 0. 03. Since the H 2 O contribution dominates in the J band, it is possible that the relatively poor fit of the cloud-free model in the J band affects the accuracy of the retrieved C/O ratio. The T-P profile shows a similar shape to that of the cloudy profile, albeit shifted by roughly 100-200 K from the model top to a pressure of a few bars, deeper than which the gradient increases in a similar Figure 16. Initial retrieval fits to the HD b spectrum. The data are down-sampled to a maximum resolution of ≈500, a factor of 4-8 lower than the original. The fit using the full data set is shown in light red. The shaded pink regions indicate where Daemgen et al. (2017) identified contiguous or near-contiguous regions of suspected high telluric contamination that could not be reliably fully removed in reduction, thus introducing potential residual systematics. The fits with each band individually are shown in the various nonred colors in each band (purple for J, green for H, and yellow for K s ). The fit using the data across all bands, but without the high-telluric regions, is shown in darker red. Both the full and single-band fits perform most poorly in fitting the data in the high-telluric regions, especially on the blue ends of each band. The retrievals are performed with the data binned down to a resolution 4-8 times lower than its original, to mitigate the potential effects of binning errors from the opacity tables. Prominent alkali lines in the J-band data are marked with vertical dashed lines, including two K I doublets and a smaller Na I line. Immediately beneath the spectra are the contribution functions for the three principal carbon-and oxygen-bearing species. The deepest contours, outlined in solid colors, enclose the regions where the contribution function reaches >1% of the total contribution within the atmospheric column at a given wavelength bin. Each successive contour denotes 2 orders of magnitude smaller fractional contribution (here, 10 -4 ). 

fashion to the cloudy profile, but with a much higher uncertainty that is also consistent with an isothermal profile.
section_name: The Cloud-free Model

section_content: In the comparison of our results from Section 5 concerning 2M 2224 with those of Burningham et al. (2017), we retrieve a C/O ratio whose confidence interval overlaps that of their posterior distributions-but our model disagrees on nearly everything else, including the radius, gravity, T-P profile, cloud properties, and absolute abundances. In the case of 2M 2224, we have the hindsight provided by Burningham et al. (2021) that, with mid-infrared data, they were able to distinguish between specific cloud particle compositions. As a result, the slope of their T-P profile decreased considerably, as well as the retrieved abundances of H 2 O and CO by roughly 0. 4 dex. Therefore, it is not surprising that different retrieval codes disagree when the available wavelength ranges differ. However, the fact that our codes nevertheless converged on the same C/O ratio is promising. A recent analysis from Rowland et al. (2023) shows that in the L-dwarf regime, and especially earlier L dwarfs, the choice of T-P parameterization (or nonparameterization) matters in the near-infrared: more restrictive or smoothed parameterizations may bias the retrieved parameters. Therefore, the difference between our retrieved atmospheric profiles and abundances and those of Burningham et al. (2017) may lie primarily in our differing choices of T-P parameterization. Additionally, both of the retrievals on 2M 2224 assume a constant abundance with pressure for all species; Rowland et al. (2023) find that one must account for nonuniform abundances in FeH, in particular due to rainout chemistry, in order to avoid biasing the retrieved T-P profile. This effect is strongest for early L dwarfs and may not bias the existing results as severely for 2M 2224 at a spectral type of L4. 5. However, it suggests that we may consider nonuniform chemistry in future retrievals for objects such as HD 106906 b, which lies at L0. 5. We discuss our current retrieval results below. 

Our retrieved C/O ratio of -+ 0. 53 0. 25 0. 15 for HD 106906 b is entirely consistent with our estimate for the C/O ratios of fellow members of the Sco-Cen association (0. 52 ± 0. 11, as in Section 2. 2). Therefore, our results with the cloudy model do not rule out a stellar-like, brown dwarf companion formation pathway for HD 106906 b. However, our model returns a T-P profile whose shape is unlikely to be entirely physical, alternating between regions of nearly isothermal behavior with regions of rapid temperature increases. Unlike the shallow wavelength dependence of the cloud opacity in our simulated models, the cloudy fits on the HD 106906 b data show a large negative exponent, indicating that clouds contribute primarily in the J band but relatively little at redder wavelengths. This means that the differences matter more in the retrieved profiles between the cloud-free and cloudy models; both fit the H-and K s -band spectra with similar quality, but the cloudy model adjusts the absolute abundances and temperatures in accordance with the cloud constraints from the J band, "breaking" the degeneracy between the T-P profile, gravity, and absolute abundances. In the retrieval of 2M 2224ʼs atmosphere, we stopped short of invoking the interpretation of Tremblin et al. (2015) to characterize the shallow thermal gradients, as our models retrieved significant cloud opacity across the wavelength range of our data. Here, in contrast, the cloudy model Figure 18. The vertical temperature-pressure profiles of the retrieved forward model fits to the HD 106906 b data. We show the MLE, median, and 95% confidence interval of the retrieved T-P profiles for the cloud-free (blue) and cloudy (red) models, the latter of which is described in Section 3. 3. Also plotted is the (cloud-free) SONORA model for a brown dwarf at the gravity and effective temperature of our best-fit cloudy model, with the change from radiative to convective behavior occurring at a few tenths of a bar. Our retrieved profiles, in contrast, vary much less in temperature with pressure down to the expected radiative-convective boundary for an object at this effective temperature. prefers little to no cloud opacity in H and K s, which keeps viable the interpretation of the data as representing a thermochemical instability driven by disequilibrium chemistry. However, the confidence in any claim, whether in the accuracy of absolute molecular abundances or in characterizing the vertical structure, is limited in the absence of wider spectral wavelength coverage that can increase the precision of the retrieved profile and also capture key signatures of cloud condensate species. 

As mentioned above, the retrieved cloud opacity is heavily biased toward shorter wavelengths, with the cloud opacity only reaching an optical depth of ∼1 in the J band. Such a strong wavelength dependence, with a power-law exponent of -7, is likely not to be attributable to a specific condensate in the atmosphere and may be a combination of some cloud opacity (e. g., SiO 2, as seen in the constraints in the condensate pressures of Burningham et al. 2021 ) and potential remaining systematics in the shortest wavelengths. This potential degeneracy is likely only resolved with broader wavelength coverage, which is proving increasingly invaluable for accurate atmospheric characterization, and/or a more sophisticated treatment of clouds, such as modeling multiple distinct cloud layers. An additional drawback to modeling clouds using a functional form for the opacity, rather than incorporating scattering from model cloud particles, is that we are not able to account for any amount of carbon and oxygen contained within 2. the clouds. Burningham et al. (2021) were able to determine that the choice of cloud model has an effect on their C/O ratio constraints at only around the 1% level, which means that the C/O ratios are consistent within their retrieved uncertainties. In their case it was primarily because they found their oxygenbearing clouds to reside primarily at pressures shallower than the photosphere, meaning that their contribution to the overall oxygen budget was ∼1%. We can make a first-order estimate of the maximum effect of silicate condensation on our C/O ratio by following the prescription of Burrows & Sharp (1999), used in Line et al. (2015 Line et al. (, 2017) ) and Burningham et al. (2021), that assumes that on average 3. 28 atoms of oxygen are sequestered per silicon atom in silicate condensates. Since our retrieved metallicity distribution is consistent with solar, as are our abundance distributions for the major oxygen-containing species when compared with solar-metallicity equilibrium models, we can estimate that a maximum of ∼16% of our atmospheric oxygen may be held in silicate clouds. Our best-fit C/O ratios would then drop to as low as 0. 55 (versus 0. 66), with the retrieved range updated to -+ 0. 45 0. 21 0. 13 (versus -+ 0. 53 0. 25 0. 15 ), still consistent with the association C/O. However, a true accounting of the oxygen budget in condensates will necessitate a more careful treatment of clouds than this work provides. 2. 

The limitations of the near-infrared in characterizing atmospheres in this temperature range are now well established. Therefore, our work serves as a preparation for future retrievals, taking advantage of a broader wavelength coverage, on this and other similar planetary-mass companions. JWST now allows high-resolution, high-S/N emission spectra of spatially resolved, very low mass companions, and the largest benefit to retrievals is its ability to extend to the mid-infrared. In the case of HD 106906 b, GTO observations with JWST are already scheduled that will capture both an R ∼ 1000 spectrum using NIRSpec (G395M, λ = 2. 87-5. 27 μm) and a lowresolution (R ∼ 100) Mid-Infrared Instrument (MIRI) LRS spectrum spanning 5-12 μm. The results of Burningham et al. (2021) have suggested that extending into the mid-infrared not only constrains specific cloud compositions but also significantly increases the range of the spectrum little affected by cloud opacity, which can allow for more accurate constraints on both the T-P profile and gas opacities. The results for 2M 2224 suggest that the relative gas abundances may be robust to limitations in wavelength range, but that one should not expect consistent gravity, T-P profile, or cloud constraints unless one has longer wavelength coverage. This being said, we are still limited to regions of the atmosphere that can be seen in emission; longer wavelengths will tend to probe cooler regions, which for directly imaged companions without thermal inversions will mean shallower pressures. The deepest parts of the atmosphere beneath the photosphere for these wavelengths, and/or beneath optically thick cloud layers, may still be inaccessible. This means that a complete picture of metrics such as the C/O ratio is likely to be still out of reach.
section_name: Discussion

section_content: We use an atmospheric retrieval code, the APOLLO code, in its first application to a cloudy L-dwarf spectrum. Our goal is to constrain formation pathways of companions to young stars. Signatures of their formation as either binary-star-like or planet-like (i. e., formed in a disk) should be imprinted in their chemistry, using metrics such as the C/O ratio and metallicity. 

From the analysis of our model results, we conclude the following:

1. Based on our self-retrieval results, the wavelength range and S/N of the HD 106906 b are sufficient to accurately constrain the C/O ratio for simulated data. Cloud-free models can retrieve a similarly accurate C/O ratio but are not preferred statistically when cloud opacity is present in the data, and in that case they may return inaccurate gravities, radii, and particularly bias toward high molecular abundances. 2. When comparing our retrieval results on the field L dwarf 2M 2224-0158 with those of Burningham et al. (2017), we find a consistency in our C/O ratios but a disagreement in the T-P profiles, cloud properties, and molecular abundances. This warrants a similar interpretation to those in Burningham et al. (2017) and Molliere et al. (2020), where a degeneracy is seen, especially in the near-infrared, between retrieved cloud properties and the T-P profile. 3. Our best-fitting model for HD 106906 b yields a C/O ratio of -+ 0. 53 0. 25 0. 15, consistent with the range of C/O ratios estimated for members of the Sco-Cen association (0. 52 ± 0. 11). This implies that we cannot rule out the hypothesis that HD 106906 b formed via the pathway expected for a brown dwarf companion to HD 106906, in contrast with a planet-like pathway. 4. However, our solution for the atmospheric emission of HD 106906 b yields negligible cloud opacity in the H and K s bands, which, along with a shallow temperature gradient at pressures less than a few bars, suggests that our results point to a cloud-temperature degeneracy. 5. As with many other retrievals of objects at similar masses and temperatures, additional data in the mid-infrared (10 μm) will be helpful in breaking the degeneracies in atmospheric structure and composition. 

JWST is currently observing directly imaged companions, obtaining spectra at resolutions 1000 at wavelength ranges not obtainable from the ground. By expanding the region where clouds are expected to contribute little, such as the thermal infrared region of ∼3-5 μm, we can better constrain the thermal structure and gas abundances, and by extension both the gravity and metallicity. Additionally, R ∼ 100 spectra are available through the MIRI, which extends the wavelength range into the realm where cloud-specific features-such as those from enstatite-are visible in emission. Follow-up observations in these wavelength ranges are planned for HD 106906 b that will allow us to employ a cloud model that more directly models specific cloud condensates. Additionally, while we have not resolved the discrepancy in brightnesses in the near-infrared between ground-and space-based observations, additional wavelength coverage into the thermal and midinfrared with JWST will also help us investigate this apparent disagreement. At the same time, accurate C/O ratios and metallicities of more companion hosts are needed to directly compare with the retrieved chemistry of the companions. In either case, retrievals on either ground-based or space-based data will benefit greatly from a set of standardized intermodel comparisons of results from various retrieval codes, to test how each model's treatment of the physics affects the inferred atmospheric properties.
section_name: Conclusions

section_name: Acknowledgments

section_content: We would like to acknowledge Dr. Natasha Batalha, whose advice and expertize on retrieval and statistical methods have been valuable in improving the rigor of this work. We also would like to thank our referees for the careful thought and feedback. 

We are grateful for support from NASA through the JWST NIRCam project though contract No. NAS5-02105 ( M. Rieke, University of Arizona, PI ). 

Software: APOLLO (Howe et al. 2017(Howe et al., 2022)), Astropy (Astropy Collaboration et al. 2013), Jupyter (Kluyver et al. 2016), Matplotlib (Hunter 2007), Numpy (van der Walt et al. 2011), Pandas (McKinney 2010), PICASO (Batalha et al. 2019), Scipy (Jones et al. 2001), Synphot (STScI Development Team 2013, 2018 ).
section_name: Acknowledgments