Cloud busting: enstatite and quartz clouds in the atmosphere of 2M2224-0158

We present the most detailed data-driven exploration of cloud opacity in a substellar object to-date. We have tested over 60 combinations of cloud composition and structure, particle size distribution, scattering model, and gas phase composition assumptions against archival $1-15 {\rm \mu m}$ spectroscopy for the unusually red L4.5~dwarf 2MASSW~J2224438-015852 using the Brewster retrieval framework. We find that, within our framework, a model that includes enstatite and quartz cloud layers at shallow pressures, combined with a deep iron cloud deck fits the data best. This models assumes a Hansen distribution for particle sizes for each cloud, and Mie scattering. We retrieved particle effective radii of $\log_{10} a {\rm (\mu m)} = -1.41^{+0.18}_{-0.17}$ for enstatite, $-0.44^{+0.04}_{-0.20}$ for quartz, and $-0.77^{+0.05}_{-0.06}$ for iron. Our inferred cloud column densities suggest ${\rm (Mg/Si)} = 0.69^{+0.06}_{-0.08}$ if there are no other sinks for magnesium or silicon. Models that include forsterite alongside, or in place of, these cloud species are strongly rejected in favour of the above combination. We estimate a radius of $0.75 \pm 0.02$ Rjup, which is considerably smaller than predicted by evolutionary models for a field age object with the luminosity of 2M2224-0158. Models which assume vertically constant gas fractions are consistently preferred over models that assume thermochemical equilibrium. From our retrieved gas fractions we infer ${\rm [M/H]} = +0.38^{+0.07}_{-0.06}$ and ${\rm C/O} = 0.83^{+0.06}_{-0.07}$. Both these values are towards the upper end of the stellar distribution in the Solar neighbourhood, and are mutually consistent in this context. A composition toward the extremes of the local distribution is consistent with this target being an outlier in the ultracool dwarf population.

stemmed from their spectroscopic signatures. In emission and transmission spectra alike, clouds mute the depth of features by providing continuum opacity that blocks flux through the otherwise brighter atmospheric windows between absorption bands. The largest and most scientifically mature sample of objects displaying this kind of spectral signature of clouds are the L dwarf population in the Solar neighbourhood, and it is these clouds that are the subject of this work. With temperatures in the range 1200 < ∼ T eff < ∼ 2300K, their spectral sequence is widely understood as arising from the impact of cloud layers first appearing and then sinking beneath the photosphere with decreasing temperature (e.g. Kirkpatrick 2005).
A range of different approaches for modelling clouds in substellar and exoplanet atmospheres are now described in the literature (see e.g. Marley & Robinson 2015;Helling & Casewell 2014, for detailed reviews). These fall into two main categories: those that deal with the detailed microphysics of cloud nucleation and grain growth, and those that do not. The pioneering and popular eddysed cloud model (Ackerman & Marley 2001) is an example of the latter. This cloud model has been used widely in both self-consistent grid models and, more recently, in retrieval frameworks (e.g. Saumon & Marley 2008;Morley et al. 2012;Gravity Collaboration et al. 2020;Mollière et al. 2020) In this model, condensation of a particular species is assumed to take place at pressures and temperatures below condensation curves derived from phase equilibrium calculations and thermochemical model grids. These condensation curves are often plotted alongside thermal profiles in the literature to show potential cloud locations, including on the figures in this paper. The balance of mixing condensable gases upwards and cloud particles downwards is set by the sedimentation efficiency parameter (f sed ), and grain sizes are set by mass balance considerations. A log-normal particle size distribution is assumed, with typical effective radii in the ∼ 1 − 10µm range arising. A similar approach is used by the BT Settl model grid (e.g. Allard et al. 2011). However, in this case the condensation, sedimentation and mixing timescales are not parameterised but instead drawn from radiation-hydrodynamical simulations.
Other approaches are built around a detailed treatment of the microphysical processes that govern cloud particle nucleation and growth along with their sedimentation and diffusion (e.g. Helling et al. 2006Helling et al. , 2008Gao et al. 2020;Powell et al. 2018). These models predict a range of particle sizes, including significant fractions of sub-micron grains, and cloud compositions that diverge from predictions based on simple phase equilibrium.
In emission spectra, such as those of free-floating L dwarfs, the effect of grey cloud can be mimicked by an isothermal temperature profile. This has led Tremblin et al. (2015Tremblin et al. ( , 2016 to suggest alternatives for the cloudy L dwarf paradigm, wherein a chemical convective instability decreases the temperature gradient in the photosphere as the carbon chemistry moves from CO to CH4 dominated with decreasing temperature through the L sequence to the T dwarf sequence. However, while a log-normal size distribution with ∼ 1 − 10µm effective radius typically seen in the Ackerman & Marley (2001) cloud model produces a roughly grey opacity for most species, smaller particles and tighter size distributions produce nongrey opacity that can be distinguished from the cloud-free model, and allows different cloud species to be distinguished.
In Burningham et al. (2017, hereafter B17) we introduced the first retrieval analysis of cloudy L dwarfs, and considered simple cloud parameterisations that allowed for grey and nongrey cloud opacity. The latter was treated as a simple function of wavelength raised to some power, i.e. τ ∝ λ α , where the power index was a retrieved parameter. Here, we build upon this work and extend our cloud model to incorporate specific cloud species via Mie scattering in an effort to make use of the full set of 1 − 15 µm spectroscopy available for a range of L and T dwarfs in archival data, including one of the targets from B17. There are three main goals of this extension : 1) to resolve the significant disagreement that was found between the retrieved thermal profiles and those of the self-consistent model grids in B17; 2) to gain empirical insight into the cloud properties of an L dwarf with respect to species and particle sizes; 3) to explore the impact of cloud treatment on estimates of composition indicators such as C/O ratio.
In this work we take a deep dive into the clouds of the red L dwarf 2MASSW J2224438-015852 (2M2224-0158 from here-on), exploring a wide-range of plausible (and implausible) cloud species and combinations thereof. In Section 2 we review the literature concerning 2M2224-0158. Section 3 outlines our retrieval framework, followed by a discussion of our model selection in Section 4. We provide analysis of the winning model and associated parameter estimates in Section 5, and discuss the implications of our results in Section 6. Conclusions are summarised in Section 7.

THE TARGET: 2M2224-0158
2M2224-0158 was originally identified in a search for L dwarfs in the Two Micron All Sky Survey ). Using a Keck LRIS spectrum the object was optically classified as an L4.5 dwarf with detectable Hα emission. Given the activity signature of the source, 2M2224-0158 was photometrically monitored and found to be variable at I band (Gelino et al. 2002), an early but strong signature of atmosphere dynamics. It was also observed for linear polarization in search of a signature of dust in the atmosphere but no detection was made Ménard et al. (2002); Zapatero Osorio et al. (2005). Dahn et al. (2002) reported a parallax for the object making it an ideal candidate for numerous detailed studies of brown dwarfs. Importantly, Cushing et al. (2005) examined the 0.6 -4.1 µm spectrum and determined that it had an anomalously red spectrum typically interpreted as indicative of unusually thick condensate clouds and/or a low surface gravity. Cushing et al. (2009) followed up on that result, obtaining 5.5 -38 µm data using the Spitzer Space Telescope's Infrared Spectrograph (IRS) and determined that there was a significant deviation from the model prediction at ∼ 9µm indicating the presence of a silicate feature. They speculated that this might be due to enstatite and/or forsterite clouds, while Helling et al. (2006) suggested that the feature could arise from quartz (SiO2) grains, in contrast to phase-equilibrium predictions (Visscher et al. 2010b). Sorahana et al. (2014) investigated the 2.5 -5.0µm AKARI spectrum of 2M2224-0158 and found that their self-consistent grid models were unable to fit the 1 − 5 µm spectrum without additional heating to raise the temperature in the upper (P < 0.1 bar) atmosphere by a few 100K relative to equilibrium predictions. The consensus of literature conclusions on the atmosphere of 2M2224-0158 is that it is an outlier compared to the field population.
The anomalous red spectrum also begged the question as to whether or not the source was young. Several groups investigated youth characteristics of this object in context with a growing sample of unusually red L dwarfs with spectroscopic signatures of a low surface gravity (e.g. 2M0355; Faherty et al. 2013 PSO 318;Liu et al. 2013, etc.). For instance, Gagné et al. (2014), Faherty et al. (2016, Martin et al. (2017) and Liu et al. (2016) examined the kinematics, near infrared spectrum, and/or color magnitude position of 2M2224-0158 and concluded that it was a field gravity L dwarf with space motion consistent with the old field population near the Sun. At present there is no obvious moving group or association that 2M2224-0158 might belong to. Therefore the anomalous observables are speculated to be atmospheric in nature and not a consequence of youth.
B17 performed the first atmospheric retrieval of 2M2224-0158. That work obtained effective temperature and log g results which matched semi-empirical results using just the near infrared spectrum of the source. Importantly for this work, B17 applied 6 different approaches to the cloud formulation of 2M2224-0158 and found that a power law cloud deck was the best fitting model. The optically thick cloud deck which passes τ cloud 1 (looking down) at a pressure of around 5 bar. The temperature at this pressure is too high for silicate species to condense, and B17 argued that corundum and/or iron clouds are responsible for this cloud opacity. The retrieved profiles were cooler at depth and warmer at altitude than the forward grid model comparisons, therefore B17 speculate that some form of heating mechanism may be at work in the upper atmospheres of 2M2224-0158. All of these conclusions were drawn in the absence of the mid-infrared spectrum which this paper now includes in the analysis.

RETRIEVAL FRAMEWORK
Our retrieval framework (nicknamed "Brewster") is an extension of that described in B17, which we have improved in respect to our cloud model and our treatment of gas phase abundances and opacities. For a more complete description of the framework and its validation we refer the reader to B17. We summarise its key features before discussing the new extensions in more detail below.
Our radiative transfer scheme evaluates the emergent flux from a layered atmosphere in the two stream source function technique of Toon et al. (1989), including scattering, as first introduced by McKay et al. (1989) and subsequently used by e.g. Marley et al. (1996); Saumon & Marley (2008); Morley et al. (2012). In Brewster's default arrangement (used here), we set up a 64 layer atmosphere (65 levels) with geometric mean pressures in the range log P (bar) = −4 to 2.3, spaced at 0.1 dex intervals.
We set the temperature in each layer via the multiple exponential parameterisation put forward by Madhusudhan & Seager (2009). This scheme treats the atmosphere as three zones: where P0, T0 are the pressure and temperature at the top of the atmosphere, which becomes isothermal with temperature T3 at pressure P3. In its most general form, a thermal inversion occurs when P2 > P1. Since P0 is fixed by our atmospheric model, and continuity at the zonal boundaries allows us to fix two parameters, we consider six free parameters α1, α2, P1, P2, P3, and T3. If we rule out a thermal inversion by setting P2 = P1 (see Figure 1, Madhusudhan & Seager 2009), we can further simplify this to five parameters α1, α2, P1, P3, T3.
In this work we consider the following absorbing gases: H2O, CO, CO2, CH4, TiO, VO, CrH, FeH, Na, K. These gases were chosen as they have been previously identified as important absorbing species in mid-L dwarf spectra, and which thus are amenable to retrieval analysis.
We calculate layer optical depths due to these absorbing gases using opacities sampled at a resolving power R = 10,000, which are taken from the compendium of Freedman et al. (2008Freedman et al. ( , 2014, with updated opacities described in B17. The line opacities are tabulated across our temperaturepressure regime in 0.5 dex steps for pressure, and with temperature steps ranging from 20 K to 500 K as we move from 75 K to 4000 K. We then linearly interpolate to our working pressure grid. We include continuum opacities for H2-H2 and H2-He collisionally induced absorption, using the cross-sections from Richard et al. (2012) and Saumon et al. (2012), and Rayleigh scattering due to H2, He and CH4 but neglect the remaining gases. We also include continuum opacities due to bound-free and free-free absorption by H − (John 1988;Bell & Berrington 1987), and free-free absorption by H − 2 (Bell 1980). The D resonance doublets of NaI (∼ 0.59 µm) and KI (∼ 0.77 µm) can become extremely strong in the spectra of brown dwarfs, with line profiles detectable as far as ∼3000 cm −1 from the line centre in T dwarfs (e.g. Burrows et al. 2000;Liebert et al. 2000;Marley et al. 2002;King et al. 2010). Under these circumstances, a Lorentzian line profile becomes woefully inadequate in the line wings and a detailed calculation is required. For these two doublets, we have implemented line wing profiles based on the unified line shape theory (Allard et al. 2007b,a). The tabulated profiles (Allard N., private communication) are calculated for the D1 and D2 lines of Na i and K i broadened by collisions with H2 and He, for temperatures in the 500 − 3000 K range and perturber (H2 or He) densities up to 10 20 cm −3 . Two collisional geometries are considered for broadening by H2. The profile within 20 cm −1 of the line centre is Lorentzian with a width calculated from the same theory.
We estimate parameters and compare models in a Bayesian framework, and sample the posterior probabilities using the emcee (Foreman-Mackey et al. 2013) algorithm. As with B17, we use 16 walkers per dimension, and extend our iterations until we have run for at least fifty times the estimated autocorrelation length. We also check that the maximum likelihood is no longer evolving with additional sets of 10,000 -30,000 iterations. Typically our emcee chains run for 70,000 -130,000 iterations. We initialise all of our model cases identically, using tight Gaussians centered on equilibrium predictions for solar composition gas abundances and log g = 5.0. Our clouds are initialised with broader distributions reflecting our ignorance of cloud locations, but roughly corresponding to the locations of condensation curves along a T eff = 1700 K log g = 5.0 Saumon & Marley (2008) grid model thermal profile. The thermal parameters are also loosely initialised around the same profile. Since we are combining data from multiple instruments, we also allow for calibration errors by including relative scale factors between the SpeX data and the GNIRS and Spitzer data, respectively. Data from each instrument also has its own tolerance parameter to reflect the range of SNR exhibited by these data sets. We use uniform and log-uniform priors which generously enclose the plausible parameter space. These are summarised in Table 2. The full set of free-parameters in a particular model proposal is called the state vector.

Gas phase abundances
In this work we incorporate two methods for handling gas phase abundances. The first method is to specify vertically constant volume mixing fractions for each gas, and retrieve these directly, as was done in B17. This is a common approach in so-called "free" retrievals. This method has the advantage of simplicity, but may fail to capture important variations of gas abundances with altitude. Figure 1 shows expected thermochemical equilibrium abundances for likely absorbing gases in the LT regime, calculated along a thermal profile taken from a T eff = 1700K, log g = 5.0, solar metallicity self-consistent Marley & Saumon grid model. It can be seen that abundances for many of our absorbers are expected to vary by several orders of magnitude within the 0.1 -10 bar region from which we expect to have large contributions of flux. Although free retrieval of gas abundances that vary with altitude may be desirable, the large number of parameters involved (e.g. 5+ per gas to provide meaningful flexibility), the necessity of allowing sharp changes in vertical profile, and the potential degeneracy with the thermal profile suggest such an approach to the inverse problem would be ill-posed. The second method we use here attempts to address this Parameter Prior gas fraction (Xgas) log-uniform, log Xgas −12.0, gas Xgas 1.0 thermal profile: α 1 , α2, P 1, P 3, T 3 uniform, constrained by 0.0 K < T i < 5000.0 K scale factor, R 2 /D 2 uniform, constrained by 0.5R Jup R 2.0R Jup gravity, log g uniform, constrained by 1M Jup gR 2 /G 80M Jup cloud top 1 , Ptop log-uniform, −4 log 10 Ptop +2.3 cloud decay scale 2 , (∆ log 10 P ) decay log-uniform, 0 < (∆ log 10 P ) decay < 7 cloud thickness 3 (∆ log 10 P ) thick log-uniform, constrained by log 10 Ptop log 10 Ptop + (∆ log 10 P ) thick 2.3 cloud total optical depth (extinction) at 1 µm 3 uniform, 0.0 τ cloud 100.0 Hansen distribution effective radius, a log-uniform, −3.0 < log 10 a(µm) < 3.0 Hansen distribution spread, b uniform, 0 < b < 1.0 Geometric mean radius (log-normal distribution), µ −3.0 < log 10 µ(µm) < 3.0 Geometric standard deviation (log-normal distribution), σ 1 < σ < 5 wavelength shift, ∆lambda uniform, −0.01 < ∆λ < 0.01µm tolerance factor, b uniform, log(0.01 × min(σ 2 i )) b log(100 × max(σ 2 i )) Table 2. Notes: 1) for an optically thick cloud deck this is the pressure where τ cloud = 1, for a slab cloud this is the top level of the slab; 2) decay height for cloud deck above the τ cloud = 1.0 level; 3) thickness and τ cloud only retrieved for slab cloud.
shortcoming by assuming thermochemical equilibrium, and retrieving [Fe/H] and C/O instead of individual gas abundances. The gas fractions in each layer are then drawn from tables of thermochemical equilibrium abundances as a function of T, P, [Fe/H], C/O ratio along the thermal profile of a given state vector.

Cloud model
Our cloud model requires two levels of specification: 1) the cloud structure and location in the atmosphere; 2) optical properties of the cloud particles that define the wavelength dependence of its opacity. For the former, we adopt the same methods for parameterising cloud structure as we did in B17. Each cloud is approximated as one of two options: a "slab" cloud or a "deck" cloud. Both are defined by the manner in which opacity due to the cloud is distributed amongst the layers of the atmosphere. The slab cloud has the following parameters: a cloud top pressure (Ptop), a physical extent in log-pressure (∆ log P ), and a total optical depth at 1µm (τ cloud ). The optical depth of the slab cloud is distributed through its extent as dτ /dP ∝ P (looking down), reaching its total value at the high-pressure extent of the slab. This cloud can have any optical depth in principle, but we restrict the prior to 0.0 τ cloud 100.0.
The deck cloud differs from the slab cloud in that it always becomes optically thick at some pressure. If the slab cloud can be thought of as a cloud that we might be able to see to the bottom of, the deck cloud is one that we only see the top of, and thus cannot infer the true total optical depth. Instead, it is parameterised by the pressure at which its total optical depth at 1µm passes unity (looking down) for the first time (P deck ), and a decay height ∆ log P over which the optical depth falls to shallower pressures as dτ /dP ∝ exp(( . The optical depth of the deck cloud increases following the same function to deeper pressures until ∆τ layer = 100. Deck clouds can become opaque very rapidly with increasing pressure, such that essentially no information about the atmosphere from deep beneath the deck is accessible. This must be borne in mind when interpreting retrieved (parameterised) thermal profiles. In deck cloud cases, the profile below the deck (and its spread) simply extends the gradient of the profile at the cloud deck (and spread therein).
We consider the wavelength dependent optical properties of different condensate species under the assumption of Mie scattering, and also test the alternative model of Distributed Hollow Spheres (Min et al. 2005) for a small subset of models. We have taken the optical data (refractive indices) for likely condensates from a variety of sources (see Table 3), and have pre-tabulated Mie coefficients for our condensates as a function of particle radius and wavelength.
Wavelength dependent optical depths, single scattering albedos and phase angles in each layer are then calculated at runtime by integrating the cross sections and Mie efficiencies over the particle size distribution in that layer for each cloud species present. The total particle number density for a given condensate in a layer is calibrated to the optical depth at 1µm in the layer as determined by the parameterised cloud structure in the state vector. In this work we explore retrieval models that assume either a Hansen (1971) or a log-normal distribution of particle sizes. Wakeford & Sing (2015) Table 3. Sources for optical data for condensates used in this work For a Hansen distribution the number n of particles with radius r is defined as: where the parameters a and b refer to the effective radius and spread of the distribution respectively, and are given by: .
In Brewster, we are able to simulate opacity arising by combinations of condensate species by combining one of our simple cloud structures with a set of optical properties to define a given "cloud". We can then incorporate any linear combination of these "clouds" into our retrieval model to arbitrary degrees of complexity.

MODEL SELECTION
In an effort to explore as much potential cloud parameter space as possible in limited computation time, we have simulated a range of structural combinations, building up in complexity from a simple no-cloud model to combinations of up to four cloud species. These arrangements typically combine one or more slabs with a single deck cloud. We have not modelled multiple species of deck cloud. This is because any cloud opacity that lies at deeper pressure than a deck cloud's optically thick region will be totally obscured. Since a slab cloud is capable of reproducing much of the behaviour of a deck cloud, we only employ a single deck cloud in each model, which is typically retrieved as the deepest cloud of the set (although it is free to appear anywhere with atmosphere according to the prior). Each cloud structure represents a single cloud species, and different slabs and a deck are free to co-locate in pressure. Our choices of model fall into two groups based on the rationale for selecting the model. These are: (i) expected species based on phase equilibrium and cloud modelling predictions (ii) species selected on basis of Mie-scattering features that overlap with features seen in the spectrum of 2M2224-0158.
It is worth noting that the second group may include implausible species. That is, species that we don't expect given our understanding of substellar atmospheres e.g. Fe2O3, which is very unlikely to be found in such a reducing environment.
In B17 we found that the cloud opacity seen there was most consistent with a Hansen distribution of particle sizes (see Equations 3,4,5). For this reason, we have focused our computational resources on models that assumes this particle size distribution, and treat scattering with Mie theory. However, for completeness we have also tested a subset of models using assuming a log-normal particle size distribution, and also DHS scattering (with Hansen distribution). In addition, for a subset of cloud arrangements focused on the top-ranked models at different levels of complexity, we tested models that assume thermochemical equilibrium gas abundances. Table 4 lists the wide range of combinations of condensate species and arrangements that we have modelled in our retrievals. We rank our models according to the Bayesian Information Criterion (BIC), in order of increasing ∆BIC from the "winning" model. The BIC provides a method for model selection for cases when it is not possible to calculate the Bayesian Evidence for a set of models. The BIC is defined as: where: k = number of parameters; n = the number of data points; L = the likelihood. Kass & Raftery (1995) provide the following intervals for selecting between two models under the BIC: 0 < ∆BIC < 2: no preference worth mentioning; 2 < ∆BIC < 6: postive; 6 < ∆BIC < 10: strong; 10 < ∆BIC: very strong.

RESULTS & ANALYSIS
The top ranked model is one that combines slabs of enstatite (MgSiO3) and quartz (SiO2), with a deeper iron deck. The ∆BIC suggests that this model is strongly preferred over other models, with the closest "runner up" having a ∆BIC = 15. Figure 2 shows the model spectrum for the median set of parameters from the top-ranked model. This shows for the first time a model fit to the 1-15µm spectrum of an L dwarf that successfully fits the feature in 8 − 10µm region of theSpitzer IRS data, suggested to be due to enstatite clouds (Cushing et al. 2008;Stephens et al. 2009). In addition, our retrieval model is able to reproduce the full shape of the spectrum between 1 and 15µm. With 36 parameters, presenting the properties of the preferred model is challenging in its own right. A "corner" or "staircase" plot such as is typically used to show correlations and degeneracies in the posterior distributions of multivariate model fits is impractical for all 41 parameters. For this reason we provide a full  Table 4. List of models tested in this work for 2M2224-0158 along with ∆BIC. Rows with no additional notes alongside the number of parameters are model runs with our default vertically constant gas fractions, Hansen particle size distribution, and Mie scattering. Changes to this are indicated in brackets as follows: CE = chemical equilibrium gas fractions; log-normal = log-normal particle size distribution; DHS = distribution of hollow spheres scattering.
41-parameter corner plot as an online plot for completeness, and here break the parameters into subgroups which are visualised separately and discussed separately. Figure 3 shows the retrieved thermal profile along with pressure depths for the clouds, and comparison to self-consistent grid models and phase-equilibrium condensation curves. The comparison model profiles are for T eff = 1700 K and 1900 K, log g = 5.0 and 5.5 atmospheres taken from the DRIFT-Pheonix (Helling et al. 2008) and Saumon et al. (2007) f sed = 2 model grids. The self-consistent grid models have been selected to bracket both our extrapolated T eff and log g, and those estimated from the bolometric luminosity of the target and radii taken from evolutionary models (Filippazzo et al. 2015). The deep thermal profile (P > ∼ 1 -10 bar) follows the gradient of the comparison model profiles well. This in contrast to what was found when using only 1 − 2.5µm near-infrared data in B17, where there was a significant offset between the retrieved thermal profile and the grid model profiles at nearly all pressures. This highlights the importance of using as broad a wavelength coverage as possible when pursuing retrieval analysis on cloudy atmospheres.

Thermal profile
Several noticeable differences are in the two relatively sharp gradient changes seen in three of the grid models, which are not seen in our retrieved profile. For example, in the Saumon et al. (2007) 1900K grid model profile, the gradient strays from the convective adiabat rising from deep in the atmosphere as thermal transport becomes radiation dominated at around 2500K, log P ≈ 1.0. The gradient shifts again close to log P ≈ 0.0 due to the formation of a detached convective zone driven by cloud opacity and the effect of the peak of the local Planck function overlapping with high gas opacity. This gives way to radiative transport again towards shallower pressures. Our five parameter profile model is unable to capture such detail, so this difference cannot be interpreted as significant. However, since these features depend on details of how and where different opacities, such as clouds, come into play (see e.g. Marley & Robinson 2015) it will be an exciting prospect to develop methods for investigating such features in the coming years.
At shallow pressures we can see significant differences between our retrieved profile and the comparison profiles. Specifically, our retrieved profile is considerably warmer at pressures < ∼ 0.1 bar, with an offset of nearly 500K from the warmest grid models at the top of the atmosphere. This is a similar result to that found in B17, and will be discussed further in Section 6.2.

Cloud properties
The retrieved parameters for our "winning" model are summarised in Table 5. Our preferred model frames the cloud opacity as arising from three condensate species: MgSiO3, SiO2 and Fe. The first two of these are parameterised in slab clouds, while the Fe is parameterised as a deck cloud (see Section 3.2). The median locations of these "clouds" in the final 2000 iterations of our emcee chain (1.1 million samples) are indicated by shaded bars to the right of Figure 3.
The two silicate clouds overlap in pressure location, and are found high in the atmosphere at pressures below 0.1 bar. The posterior distributions for their cloud bases show spreads of ∼ 1dex in log P , reflecting relatively weak constraints on their depths in the atmosphere. The Fe cloud deck becomes optically thick at deeper pressures close to 10 bar. The location of the τFe = 1.0 pressure level appears to be tightly constrained to log P (bar) = 0.95 ± 0.06dex. The pressure change over which the optical depth of the deck cloud drops to 0.5 is constrained to d log P = 4.20 +1.95 −2.25 dex, such that τ deck = 0.5 at log P = −3.25 +2.25 −1.95 . This thus represents a cloud opacity that is quite widely distributed through the atmosphere, despite the apparently tight constraints on its location. This decay height is implemented in linear pressure space (see Equation 2, so this corresponds to a decay pressure scale of just under 10 bar, and thus τ deck = 4.0 at approximately 29 bar (log P = 1.46).
In addition to the profile and median cloud locations, Figure 3 also shows the condensations curves for several species, calculated for solar composition gas using the same methods as Visscher et al. (2010b). Condensation curves such as these do not necessarily provide a prescription for what will condense, but instead provide a useful guide as to what can condense. The pressure depths of our silicate clouds are both consistent with being located to the left of their respective condensation curves on our thermal profile. The τ = 1.0 reference pressure for our deep Fe cloud deck corresponds to a temperature of of approximately 2500K. This is more puzzling, and will be discussed in more depth in Section 6.4.
Our top-ranked models all use Mie scattering clouds with a Hansen distribution of particle sizes, and models using either DHS scattering, or log-normal particle size distributions are very strongly rejected. Our top-ranked models also all include enstatite in its amorphous form, although our 6thranked model includes a mixture of amorphous and crystalline grains. All three retrieved clouds are dominated by sub-micron grains and have a negligible number of particles sized larger than 1µm (see Table 5). Particle distributions dominated by small particles were a feature of the retrieved parameters for all the cloudy models tested.
To gauge the plausibility of our retrieved cloud properties and provide suitable outputs for comparison to physically motivated cloud models, we have estimated the amount of material required to produce the retrieved clouds in our best fitting model. The simplest approach for this is to use the retrieved total optical depth at 1µm and the particle size distribution parameters to estimate the column density for each cloud. We find a MgSiO3 column density of 5.4 +1.1 −0.7 × 10 −4 g cm −2 , or by number 3.2 +0.7 −0.4 × 10 18 cm −2 . For SiO2 we find a column density of 1.4 +0.6 −0.2 × 10 −4 g cm −2 , or 1.4 +0.6 −0.2 × 10 18 cm −2 .
If we take our median cloud locations at face value, we can estimate that such cloud masses would account for all of the oxygen in these atmospheric layers, and would require significant enhancement of magnesium (≈ 5×) and silicon (≈ 10×) compared to solar ratios. However, our retrieved cloud locations are not so precise, and are (at best) a very broad brushed 1-dimensional analogue to the dynamic clouds in our target's 3-dimensional rotating atmosphere. As such, it is may be useful to consider these estimated cloud masses through comparison to the available material above their condensation curves on our retrieved thermal profile. In this view, the   Saumon & Marley 2008). Also plotted are phase-equilibrium condensation curves for various potential condensates. The clouds pressures are indicated in bars to the left of the P/T profile. Blue shading indicates the median cloud location for the two slab cloud, with grey shading indicating the 1σ range. The deck cloud is furthest to the right, with the uniform blue shading indicating the optically thick extent of the cloud deck, and grey shading the 1σ interval for the τ = 1 location. Graduated shading shows the range over which the deck cloud optical depth drops to 0.5. silicate clouds account for approximately 30% of the available oxygen. If we (naïvely) assume that there are no other sinks for magnesium or silicon and that the relative proportions of MgSiO3 and SiO2 reflect the proportions of Mg and Si, we can estimate the this requires (Mg/Si) = 0.69 +0.06 −0.08 . Figure 4 shows the corner plot displaying posterior distributions for the retrieved gas fractions for absorbing gases in our model and log g, along with derived values for radius, mass, atmospheric metallicity and C/O ratio, and extrapolated T eff . The derived and extrapolated values are found as follows. Radius is defined using the retrieved model scaling factor and Gaia parallax. Mass is then found using the derived radius and the retrieved log g, L bol is found by extrapolating the retrieval model to cover the 0.5 − 20 µm range, summing the flux, and scaling it by 4πD 2 , where D is the distance defined by the Gaia parallax. T eff is then found using the extrapolated L bol and inferred radius. The atmospheric C/O ratio is estimated by assuming all carbon and oxygen in the atmosphere exist within the absorbing gases considered   Figure 5) and Sonora-Bobcat cloud-free evolutionary models (Marley & Saumon 2020, ;Marley et al in prep.). Our extrapolated bolometric luminosity for 2M222-0158 is log(L bol /L ) = −4.146±0.003. This is similar to, but slightly higher than, the empirical value found by Filippazzo et al. (2015) of log(L bol /L ) = −4.16 ± 0.01. Our derived radius, however, is considerably smaller than predicted by evolutionary models such as the Sonora-Bobcat grid plotted on Figure 5, or the models used by Filippazzo et al. (2015), and this leads to a higher T eff . This discrepancy will be further discussed in Section 6.3.

Composition
Our retrieved gas abundances provide constraints on the composition of the photosphere of 2M2224-0158. Through consideration of the chemical, condensation and dynamic processes at work in its atmosphere, these in turn can provide useful insights to its bulk composition, as derived from its natal environment. Studies of brown dwarfs such as 2M2224-0158 can serve as essential due diligence for developing and validating such methodologies. The most commonly deployed, and blunt, instrument for discussing composition in astronomy is aggregated metallicity [M/H], for which [Fe/H] is often used interchangeably as a proxy. Whilst our chemical equilibrium model allows us to directly access the metallicity of our target, our winning model requires us to estimate it based on the snapshots of photospheric abundance that our vertically constant gas fractions effectively provide. As discussed in the previous section, many of the species we are considering are expected to condense beneath and within our photosphere, as can be seen from the large decrease in atmospheric abundance of FeH, TiO and CrH with decreasing pressure in Figure 1. Our retrieved gas fractions reflect these condensation processes, and would lead to a significantly sub-solar metallicity estimate if used in isolation. For example, based on our retrieved FeH gas fraction alone, we would estimate a value for [Fe/H] = −2.2 ± 0.06. These gases are, however, relatively minor contributors to our metallicity estimate, which is dominated by carbon and oxygen. As such they can be included or ignored with negligible impact on our metallicity estimate of [M/H] = 0.38 +0.07

−0.06
Of particular interest is the C/O ratio, which is viewed as a key diagnostic for investigating the formation environments and mechanisms for giant exoplanets (e.g. Öberg et al. 2011;Gravity Collaboration et al. 2020). We estimate an atmospheric C/O ratio for 2M2224-0158 of 0.83 ± 0.01 based on the assumption that all carbon and oxygen are tied up within the considered absorbing gases. This atmospheric C/O ratio does not account for any oxygen that may be tied up in condensates such as MgSiO3 and SiO2. Our retrieved silicate cloud locations are at much shallower pressures than most of the gas opacity, < ∼ 1% of the atmosphere exists above the point where the MgSiO3 and SiO2 condensation curves cross our retrieved thermal profile. Any correction to the oxygen abundance due to clouds above this point will be small.
Our standard selection of gaseous absorbers in these experiments did not include SiO, which has very weak opacity and has not previously been identified in L dwarf spectra. However, predictions from thermochemical grids (e.g. Visscher et al. 2010b) suggest that SiO should be a significant home for oxygen in L dwarf atmospheres after CO and H2O. If we account for oxygen tied up in SiO gas at the level predicted by our thermochemical equilibrium grid for the thermal profile of our winning model, we estimate that its inclusion might be expected to increase our oxygen abundance by about 8%, reducing the C/O ratio to approximately 0.75.
To further investigate this, we have run an addition retrieval model based on our winning cloud model for which SiO opacity is included. In this case, the SiO abundance is only weakly constrained with log fSiO = −3.59 +0.33 −6.13 , consistent with its small absorbing cross-section. This results in a C/O ratio of 0.83 +0.06 −0.07 . We thus adopt this wider uncertainty range for our C/O ratio.
Our uncorrected atmospheric C/O ratio is significantly super-solar, and is amongst the highest C/O ratios found in the stellar population (e.g. Nissen 2013; Pavlenko et al. 2019;Stonkutė et al. 2020). However, it has also been found that metallicity and C/O ratio are positively correlated due to oxygen being relatively less abundant in high-metallicity thin disk stars Nissen et al. (2014). In this context, our re- trieved C/O ratio and metallicity appear roughly consistent with one another.

Comparison to B17
The future capabilities of JWST will offer wide and detailed wavelength coverage of a wide range of substellar and exo- planetary atmospheres. This work can provide data-driven insights into how additional wavelength coverage can impact retrieval studies of cloudy atmospheres in the temperature range considered here.
The most notable difference is in the shape of our retrieved thermal profile. Figure 6 compares our "winning" model thermal profile to those of lower-ranked models covering a wide range of cloud assumptions, the retrieved profile from B17, and an example self-consistent grid model selected to match our retrieved parameters for 2M2224-0158. The difference between our winning model profile using the full 1 -15µm range of available spectroscopy and that found using just 1 − 2.5µm Spex data is clear. The profile retrieved here provides a close match to the self-consistent grid model at pressures deeper than around 1 bar, whereas the B17 profile disagrees with the grid profile at all pressures. This represents a shift of over 500 K at the bottom of our modelled atmosphere. Some of this difference may arise due to the more complex cloud model employed here. However, it is striking how similar the profiles for our highly ranked models are, despite differing cloud species, and structures. This is due to the fact that the large wavelength range allows for regions that have little cloud opacity to set the gaseous opacities and define the shape of the thermal profile. In B17 this region was much smaller (since we only covered 1 − 2.5 µm), and a greater portion of the coverage was also impacted by problematic opacities such as the pressure broadened wings of the 7700Å K i resonance line. This highlights the benefit of broad wavelength coverage for remote sensing thermal structure in atmospheres, particularly in the absence of well understood cloud properties.
Along with a new thermal profile, we see a decrease of approximately 0.4 dex in our H2O and CO gas fractions. Since these species dominate the C and O budget and each is shifted by a similar amount, our estimated C/O ratio is largely unchanged. This shift in gas fractions is likely linked to the new profile. For the most part it is a steeper temperature gradient through the photosphere, thus requiring less absorbing gas to achieve a similar reduction in emergent flux. Since H2O and CO are the biggest contributors to gaseous opacity, it is not surprising to see them reduced. This reduction brings the overall estimated metallicity for this object much closer to the expected range for the Solar neighbourhood. This again highlights the importance of wide wavelength coverage to constrain the thermal profile and allow accurate abundance estimates.
We also see a significantly altered alkali abundance, with a -1.2 dex shift in our combined Na and K fraction, bringing it into conflict with the otherwise super-solar metallicity. The K opacities are particularly problematic due to the poorly constrained behaviour of the pressure broadened wing of the 7700Å K i resonance line. This hinders our models' ability to simultaneously fit narrow K features and the shape of the pseudo-continuum around 1µm as set by red wing of the 7700Å K i line, and also may impact our ability to retrieve accurate abundances. Physical interpretation of this low abundance of alkali species would thus be premature. This is an ongoing issue that that drove the decision to exclude wavelengths blue-ward of 1µm in this work, along with others (e.g. Line et al. 2015Line et al. , 2017Burningham et al. 2017;Gonzales et al. 2020). Kitzmann et al. (2020) did not exclude this region and noted problems fitting the 1µm Y band peak in their study of T dwarf. They used an updated description of the 7700Å K i line wings from Allard et al. (2016), suggesting that the challenge of accurate opacities in this region remains.
The retrieved radius is also significantly smaller than in B17, which was consistent with predictions of evolutionary models. This will be discussed in more detail in the Section 6.3.

Stratospheric heating
At shallow pressures, our retrieved thermal profile is significantly warmer than the predictions of self-consistent grid models (see Figure 3). Moreover, Figure 6 demonstrates that this is common feature of our highly ranked models. Although there is some scatter in the model profiles, the preference for temperatures above around 1200K is clear. Our top-ranked thermal profile decreases from around 1600 K at log P (bar) = −0.5, to roughly 1300 K at log P (bar) = −4.0. This is in contrast to the self-consistent models, which exhibit temperatures falling from 1600 K to below 900 K in the same range (see Figure 3). This comparison is for cloudy solar composition self-consistent models, since we do not yet have access to cloudy models that span the metallicity and C/O ratio inferred for our target. However, cloud-free models show very small impact (∼ 10 K) on stratospheric temperatures due to changing composition. For a 1900 K, log g = 5.0 model, an increase in metallicity to [M/H] = 0.5 reduces the temperature at log P (bar) = −4.0 by 15 K compared to solar composition. Increasing the C/O ratio to 0.9, raises the temperature in the same region by 10 K.
It is worth considering if this divergence from the selfconsistent model predictions is perhaps some artefact of our adopted parameterisation of the thermal profile. In B17 we tested the retrieval framework using this parameterisation against simulated data created using thermal profiles from the self-consistent model grid plotted in Figure 3, and did not identify any bias to higher temperatures in the upper atmosphere. We also tested a retrieval using this parameterisation against the benchmark T8 dwarf G570D, and found a retrieved profile consistent with grid models and the retrieval of Line et al. (2015).
We have also tested the parameterisation employed by Line et al. (2015Line et al. ( , 2017. This uses a low-resolution 13-point thermal profile with spline interpolation to the full-resolution pressure scale, but penalizes the second derivative of the final curve in the retrieval. This has the effect of minimizing the jaggedness of the profile unless the benefit to improving the fit is significant. However, the L dwarf spectrum is relatively low-contrast compared to that of the T8 dwarfs for which the method was devised. As we previously found in B17, the result is that the data do not justify anything other than an entirely linear profile, even with the inclusion of the long-wavelength data used here. A more in-depth exploration of alternative ways to parameterise the thermal profile is beyond the scope of this work. We can thus not rule out the possibility that alternate parameterisations may find different profiles that may also fit the data. As we noted in B17, the warm retrieved stratospheric temperature is consistent with previous model fitting by Sorahana et al. (2014), who found evidence for a roughly isothermal profile, with T = 1445 K, at pressures shallower than log P (bar) < ∼ − 0.5 for this object. These results thus tell a similar story for the upper atmosphere of this object, and suggest an energy transport mechanism that is not currently incorporated in self-consistent radiative-convective equilibrium models. The nature of this heating is uncertain, but it may be common amongst L dwarfs (Sorahana et al. 2014, B17). The detection of chromospheric activity via Hα emission in several objects with apparent shallow-pressure heating (including 2M2224-0158) lead Sorahana et al. (2014) to argue that the mechanism is magnetic-hydrodynamic in nature, and it is understood in the context of stellar chromospheric heating. The apparent presence of cloud layers in the heated region, however, is more reminiscent of the unexpectedly warm stratospheres and thermospheres of Solar System gas planets (e.g. Appleby 1986;Seiff et al. 1997;Marley & McKay 1999). These heated regions are not thought to be due to auroral or irradiation, but are currently understood in terms of gravity waves that propagate from deeper in the atmosphere before breaking in the stratosphere and depositing their energy (Schubert et al. 2003;Freytag et al. 2010;O'Donoghue et al. 2016).

Radius
The radius implied by our top ranked retrieval model is significantly smaller than predicted by evolutionary models (see Figure 5. This is a surprising, and difficult to interpret result which is a significant departure compared to our previous study of 2M2224-0158 in B17, which found a radius that was consistent with evolutionary models.
An examination of the retrieved radii across our ensemble of models finds that our winning model is towards the smaller end of the distribution, which has a median of 0.82 RJup, and 16th-and 84th-percentiles of 0.76 RJup and 0.97 RJup respectively. The upper end of this distribution is consistent with the predictions of various evolutionary models. However, the retrieval models with the larger radii are some of the most poorly ranked models in our set, with ∆BIC > ∼ 100. To further investigate this issue we have rerun our winning model with a Gaussian prior on the radius with mean of 0.82 RJup, σ = 0.04 RJup. The retrieved parameter set was nearly identical to that retrieved with a 0.5-2.0 RJup uniform prior on the radius, with all values consistent within 1σ, and most much closer than that. The retrieved radius with the Gaussian prior was 0.77 ± .02 RJup, compared to 0.75 ± 0.02 RJup with a flat prior.
Rerunning all of our models with the Gaussian prior is computationally prohibitive. However, we can nonetheless consider whether this Gaussian prior would be likely to impact our model selection, and rank a model that retrieved a radius closer to the expected value more highly than our winning model. To do so we must assume that the Gaussian radius prior will have similarly small effects on retrieved parameters of the other models, and that their fits to the data are similarly good/bad as with the uniform prior. The large deviation from our expected radius for our winning model was penalised in the likelihood, and increased BIC value by 16. Models with radii closer to the center of the Gaussian prior will be penalised to a lesser degree, and so will have smaller increases to their BIC values. However, since none of the large radius models have ∆BIC < 100, we can conclude that introducing a Gaussian prior would have a minimal impact on our model selection results.
Elsewhere, a similarly small radius was found by Gonzales et al. (2020) for the low-metallicity s/dL7 dwarf SDSS J1416+1348A, also using the Brewster framework. Using the Helios.R-2 framework, Kitzmann et al. (2020) found a smaller than expected radius for the T1 dwarf Indi Ba under both chemical equilibrium and vertically constant composition assumptions, and for the T6 dwarf Indi Bb under the vertically constant composition assumption.
By contrast, retrieved radii for very late type T dwarfs have generally been found to be in agreement with evolutionary models. For example, the T7.5p companion SDSS J1416+1348B was found by Gonzales et al. (2020) to have retrieval radius consistent with low-metallicity Sonora-Bobcat models using Brewster. Using the CHIMERA package, Line et al. (2017) found retrieved radii for 9 out of 11 T7 and T8 dwarfs were consistent with the evolutionary model estimates used in Filippazzo et al. (2015). Sorahana et al. (2013) also found smaller than predicted radii for L dwarfs and consistent radii for T dwarfs via fitting of self-consistent models to AKARI spectroscopy. This suggests that this radius problem is not an artefact of some bias in data-driven retrieval frameworks. Kitzmann et al. (2020) have suggested that the small retrieved radii in their study of Indi BaBb could arise from a heterogeneous atmosphere, presumably with dark regions that contribute little flux and resulting in a smaller effective emitting area. Such a scenario is consistent with the apparent lack of a radius problem in late-T dwarfs that are thought to have essentially cloudless photospheres. In the case of 2M2224-0158, it would require the equivalent of some 30% of the top of the atmosphere contributing zero flux to the spectrum to account for our difference with a typical field age evolutionary model. In the cases studied in Kitzmann et al. (2020), the equivalent zero flux covering area would be much larger, as their inferred radii are as small as 0.5RJup. Such a substantial covering fraction might be expected to produce some signal of rotational variability, as seen across so much of the brown dwarf population. However, neither Indi BaBb, SDSS J1416+1148A, nor 2M2224-0158 have convincing detections of rotational modulation (Hitchcock et al. 2020;Khandrika et al. 2013;Metchev et al. 2015;Miles-Páez et al. 2017, Vos private communication). If dark patches are indeed present, this suggests that either these objects are presenting with unfavourable geometry for detecting rotational signals, or the dark regions are arranged latitudinally in bands.
If such dark regions correspond to differences in cloud cover, then this cloud cover must be quite different to that suggested by our retrieval results. Most of our wavelength region is not strongly affected by the cloud opacity we have retrieved, so the additional cloud cover would need to be essentially grey by comparison and lying higher than around 0.1 bar, in order to effectively reduce the total flux without significantly altering the shape of the SED. This would require much larger particles than those implied by our retrieval results thus far.
It is also possible that our retrieved radius is correct, and that the shortcoming arises in the evolutionary models. The Transiting Exoplanet Survey Satellite (TESS) has revealed several new transiting brown dwarfs in the 40−70 MJup range with ages in the 3 − 6 Gyr range which have radii close to, or smaller, than that predicted for 10 Gyr old objects (Carmichael et al. 2020b,a). Carmichael et al. (2020a) speculate that this could be accounted for by the brown dwarfs reaching their asymptotic minimum radius of 0.75 RJup by the age of 4 or 5 Gyr, rather than 10 Gyr as predicted by most evolutionary models. So far, this has only been seen for objects above about 40 MJup. If it is indeed the case that only higher mass brown dwarfs are affected then this would also give rise to the difference in results between L and late T dwarfs seen via spectroscopy, as at field ages these correspond to higher and lower mass populations. As the sample of transiting brown dwarfs grows the origin of this issue will become clearer.

Cloud composition and location
Our investigation of the clouds in 2M2224-0158 explored a wide range of possibilities to make our study as data-driven as possible. Our parameterisation of cloud opacity is intentionally independent of models of cloud condensation so that we can cleanly assess what input the data alone can make to said modelling efforts. We used this flexibility to try a wide range of condensates including those predicted by cloud models, and some which are not expected but whose optical properties had the potential to match features in the data. Our preferred model combines a deep iron cloud with enstatite and quartz clouds at lower pressures. Whilst other condensates are likely present, this result suggests that these are the condensates that are dominating the cloud opacity in the photosphere. These results provide a powerful test for different models of cloud formation and condensate chemistry. A notable area of divergence between differing model approaches is the compositions of the predicted clouds.
A striking aspect of our preferred model is the presence of quartz (SiO2), with abundance comparable to enstatite. Helling et al. (2006) made a clear prediction for the presence of quartz in brown dwarf atmospheres as a dominant component of the silicate clouds. This was in contrast to the predictions of phase equilibrium modelling that suggested that quartz should not exist in substellar atmospheres. Visscher et al. (2010b) found that enstatite (MgSiO3) removes silicon from the gas phase so efficiently, that quartz can only condense if enstatite formation is suppressed somehow i.e. they found that quartz could not form alongside enstatite.
These predictions were made on the assumption of solar composition for which Mg/Si ≈ 1. In fact, the same mass balance and phase equilibrium considerations suggest that quartz is expected to co-exist alongside enstatite if silicon is more abundant than magnesium, i.e for atmospheres with an Mg/Si ratio of less than 1. In Section 5.2, we estimated Mg/Si = 0.69 +0.06 −0.08 , assuming no other sinks for Mg and Si. Mg/Si has been found to be negatively correlated with metallicity, and our estimate is comparable to outlier values seen for high-metallicity outliers in Adibekyan et al. (2015) and Suárez-Andrés et al. (2018).
Another notable feature of our preferred cloud model is the absence of forsterite (Mg2SiO4). The presence of forsterite is predicted by both microphysical and phase equilibrium models at solar composition (e.g. Helling et al. 2006;Visscher et al. 2010a;Gao et al. 2020). Our model set includes a 4 cloud model comprising quartz, enstatite and forsterite slabs, along with a deeper iron cloud. In this case, the retrieved forsterite cloud was was either pushed out of view into the deeper atmosphere, or found at lower pressure with very low-optical depth i.e. it did not contribute to the emergent spectrum, and the model was down-ranked due to additional parameters.
To investigate whether forsterite is expected in this case, we have performed preliminary calculations of phase equilibrium cloud masses for our estimated values of [M/H], C/O and Mg/Si following Visscher et al. (2010b). We find that forsterite is expected to condense over a narrower temperature range than at solar composition (at 1 bar): 1510 < T < 1668 K for our composition, versus T < 1702 K for solar. It is also predicted to be less abundant than enstatite, but more abundant than quartz. These preliminary calculations appear to present a conflict between our retrieved cloud properties and the predictions of phase equilibrium, and highlight the importance of exploring non-solar elemental abundances and/or disequilibrium effects in future cloud models.
The presence of a deeper iron cloud is consistent with the predictions of both phase equilibrium and microphysical models discussed above. However, the reference pressure at which the iron deck cloud reaches τ = 1.0 at 1µm is slightly shallower than 10 bar, and corresponds to a temperature of approximately 2500K on our retrieved thermal profile -too hot to match any of the condensation curves plotted on Figure 3, including iron. We found that this T-P point was also fairly consistent amongst the lower ranked models regardless of the compositions of shallower cloud species.
There are several things to consider in understanding this apparently problematic result. The wavelength dependence of the iron cloud opacity, and the impact of gaseous opacities, means that this deep cloud deck's impact on the spectrum is greatest near 1µm (See Figure 8). As such, the retrieved reference pressure may be affected by issues with the problematic pressure broadened wing of the 7700Å K i resonance line. It is worth noting the deeper reference pressure for Fe cloud found here compared the cloud deck in B17 occurs alongside a reduction in the K abundance. Additionally, the Fe cloud opacity builds up through a wide of range of pressures, including regions to the left of the iron condensation curve. So, it is possible that the opacity is arising from several species, with iron producing the dominant spectral signature, and this could account for some of the extension beyond the iron condensation region.
Somewhat speculatively, we could also imagine that deeper opacity could arise from molten iron raining out of the atmosphere and falling to deeper pressures before being vaporised. Opacity due to falling, and vaporising, rain can be observed on Earth in the form of virga, and such a phenomena taking place in this atmosphere would explain opacity arising beyond the condensation curve.

Cloud particle properties
Our retrieval has also placed constraints on the particle size distribution for our modelled clouds, as well as testing alternative assumptions for modelling the scattering from our cloud particles. All of our highly ranked models share key features with our top model: a Hansen distribution dominated by sub-micron grains, whose optical properties are modelled according to Mie theory rather than as a distribution of hollow spheres (DHS, Min et al. 2005). Helling et al. (2006) suggested that DHS scattering should be considered for modelling brown dwarf cloud opacity to account for inhomogeneous particle shapes. Our retrieval's preference for Mie scattering appears to rule this out. The preference for a Hansen distribution of particle sizes is consistent with the findings of Hiranaka et al. (2016), and contrasts with the popular log-normal size distribution used some self-consistent models (e.g. Ackerman & Marley 2001;Saumon & Marley 2008) and retrievals (e.g. Gravity Collaboration et al. 2020).
A significant difference between our particle size distribution and those of nearly all equilibrium and microphysical cloud models is that ours is uniform with altitude in the atmosphere. This presents a significant hurdle when comparing our results with cloud model predictions. However, some useful comparisons can be made when we consider predicted particle sizes in the context of the atmospheric levels which contribute most to the shape of the emergent spectrum, as it is these regions that will drive the retrieval estimates for particle sizes.
Figure 7 compares our retrieved particle size distributions to a log-normal distribution with parameters µ = ln(5.8µm) and σ = 2. These parameters are as predicted for enstatite by the eddysed cloud model (Ackerman & Marley 2001) used in the Saumon & Marley model grids, for the τ = 2/3 level in a T eff = 1700 K brown dwarf. All three retrieved clouds are dominated by much smaller grains than those predicted by the eddysed model. In the case of the silicate grains, this difference is what allows our retrieved model to fit the broad spectral feature in the 9 -10 µm region. The Mie scattering features are only apparent for relatively narrow distributions of small particles, whereas the broad log-normal distribution, incorporating large particles, smears out these Mie features. Figure 7 also shows that for the enstatite cloud, with a large value for the Hansen b parameter (width), the distribution tends towards a power law. This is a similar functional form to the interstellar grain size distribution employed by the cond, dusty and Settl models.
A common feature of cloud models such as EddySed and microphysical cloud models is a tendency to predict larger particle sizes deeper in the atmosphere and smaller particles at shallower pressures.
Work by Gao et al. (2018) using the CARMA microphysical cloud model found multi-modal size distributions were predicted for KCl clouds in cool (T eff ≈ 400K) substellar atmospheres, with the main peak in the distribution for very small (∼ 10 −3 µm) particles found alongside a significant peak in the abundance of larger particles (∼ 1µm). Powell et al. (2018) also found multi-modal distributions for silicate clouds in hot-Jupiter atmospheres using the same modelling framework. In that case large (< 1µm) sized grains were predicted alongside sub-micron particles. Whilst neither of these predictions are directly applicable to the case studied here, they highlight the possibility of further complexity that could be missed by a retrieval model that assumes a single particle size distribution, which is uniform with altitude.
We also pursued a limited investigation into the presence of crystalline versus amorphous silicate grains. This was driven by the presence of apparently sharp features and structure in the spectral feature attributed to silicate clouds around 9µm, which would be consistent with crystalline grains. Our highly ranked models all feature amorphous grains of enstatite, which is consistent with the predictions from microphysical modelling by Helling et al. (2006).

Clouds and C/O ratios
Given the ubiquity of clouds and the aforementioned interest in measuring C/O ratios of substellar objects and giant exoplanets, we now devote some discussion to considering how the presence of the former can impact the latter.
The impact of condensation processes on the atmospheric C/O ratio should be considered carefully. Line et al. (2015Line et al. ( , 2017) applied corrections to their retrieved C/O ratios for late-T dwarfs on the assumption of 3.28 oxygen atoms per silicon atom being removed from the atmosphere through condensation in MgSiO3 and Mg2SiO4 Burrows & Sharp (1999). This translates as a roughly 25% addition to the retrieved atmospheric oxygen abundance to be included in the C/O ratio calculation. It should be noted that this 25% correction arises by assuming a solar Si/O ratio, and that this correction may thus need fine-tuning to account for a) different dominant condensates to those assumed for above; and b) non-solar abundance ratios.
We have detected MgSiO3 and SiO2 clouds in the upper atmosphere of 2M2224-0158. As discussed in Section 5.2, our median cloud masses account for around 30% of the oxygen in the atmosphere above the condensation pressures for the detected species. However, as these oxygen bearing clouds are at much shallower pressures than most of the gas opacity, any correction to the oxygen abundance due to clouds above this point will be small. A straightforward scaling of the oxygen sequestered in the detected clouds to account for the < ∼ 1% of the observable photosphere that they occupy reduces the correction to around 0.3%. We thus neglected a cloud correction for the C/O ratio in this L dwarf. For cooler objects, such as late-T dwarfs, for which the entire photosphere lies above the condensation zone for these species we would expect the correction to correspond directly to the 30% oxygen estimated to be taken up by silicate clouds in Section 5.2. This is comparable to the 25% condensation correction estimated by Line et al. (2015) for late-T dwarfs based on equilibrium cloud condensation expectations outlined in Burrows & Sharp (1999).
Another point to consider is what impact the choice of cloud model has on the estimated C/O ratio. In this work, it appears it has surprisingly little impact on our estimated C/O ratio. Most of our retrieval estimates for the C/O ratio lie quite close to that of our preferred model, across a wide range of cloud models, including cloud free, and hence different thermal structures. Taking the mean C/O ratio of all our model runs, we find C/O = 0.85 ± 0.03, compared to 0.83 ± 0.01 for the preferred case. As can be seen from Figure 8, the bulk of our retrieval model's cloud opacity lies beneath the gas opacity for most of the spectral range considered here. A similar result was found by Mollière et al. (2020) for the directly imaged exoplanet HR8799e. In that case they found that their retrieved C/O ratio was consistent between two different cloud models that resulted in quite different thermal structures.

Chemistry
A surprising outcome of this work is that the preferred model is one which assumes vertically constant gas fractions, as opposed to one which assumes thermochemical equilibrium abundances of absorbing gases. All available T eff estimates for our target, whether based on bolometric luminosity or extrapolated from our retrieval, suggest that it lies in the regime where the chemical timescales in its photosphere are expected to be fast compared to the mixing timescales.
As we saw in Figure 1, the equilibrium mixing fractions of several of our absorbing gases are expected to vary by orders of magnitude through the photosphere, mainly due to condensation processes. As a result, vertically constant mixing ratios might be expected to struggle to fit features arising from different pressure levels in the atmosphere. In Figure 9 we highlight such a case. The FeH features in the peak of the H band are well fit by our preferred model, while the FeH features at shorter wavelengths are quite poorly fit. The contribution function, also shown in Figure 9, demonstrates  that these two regions of the spectrum are influenced by gases at different depths, with the J band spectrum dominated by contributions from deeper pressures than the H band spectrum.
In Figure 10 we plot our retrieved gas mixing ratios along with predictions from our grid of thermochemical equilibrium models interpolated for our derived metallicity and C/O ratio. As can be seen, the equilibrium FeH fraction is expected to increase rapidly with pressure through the photosphere, so the failure of our vertically constant mixing ratios to fit both sets of features is not unexpected in this context. We also note that our retrieved FeH fraction intersects the equilibrium pre- Figure 9. Comparison of our preferred model J and H band spectrum with the data for 2M2224-0158. The 1σ spread in the model distribution is shown in green shading, and the errors for the data are shaded in orange. Both are mostly covered by the line width except in a few places. Features responsible for regions of particularly poor fit are annotated. Also plotted, with its y-scale on the right hand side, is the contribution function. The τ = 1 line for gas-phase absorption is indicated with a cyan line. diction at a pressure just less than 1 bar. This is a shallower pressure than the peak in contributions for the FeH-affected spectral regions It follows then that the FeH abundance may be somewhat lower than predicted by the equilibrium grid at the peak contribution pressures for those features. This may contribute to rejection of the thermochemical equilibrium model.
Our model infers a ∼ 1.5 dex lower abundance for the combined Na+K fraction than predicted by the thermochemical equilibrium grid for our extrapolated metallicity and C/O ratio. Our preferred model is also unable to fit the narrow K i features in the near-infrared while successfully fitting the pseudo-continuum around 1µm that is set by the pressure broadened wing of the 7700Å K i resonance line (see Figure 9, despite an expected roughly vertically constant mixing ratio. This suggests that one reason the thermochemical equilibrium version of our model was rejected may lie in the ongoing challenge posed by correctly modelling the pressure broadened wings of the alkali lines, rather than indicating that an assumption of thermochemical equilibrium is incorrect for potassium in this type of object, or that there is some deficiency in our chemical grid.
Our vertically constant gas fractions agree well with the predicted fractions for H2O and CO, though our retrieved estimate for H2O is somewhat smaller than the equilibrium prediction deeper than around 1 bar. Our retrieved CO and H2O fractions are also both lower by 0.5 dex compared to our estimate in B17, which used NIR data only, resulting in better agreement with thermochemical models.
By contrast, our retrieved fractions for CH4 and CO2 are both considerably higher than the equilibrium predictions. These non-equilibrium fractions may arise due to rapid vertical mixing, whereby CO and CH4 are quenched at an abundance from a deeper level in the atmosphere (although CO quenches close to its equilibrium value as the dominant Cbearing gas). In this scenario, CO2 also displays an enhanced non-equilibrium abundance due to continued reaction equilibria with CO, which increases the CO2 fraction in the layers above the quench level (e.g. Visscher et al. 2010a). Quenching of CH4 at temperatures of around 1500 K at 1 bar has been been predicted (Visscher & Moses 2011;Moses et al. 2011), but is not included in our thermochemical grid, and the implied quench temperature in this case is somewhat higher for our inferred [M/H] and C/O values. Although the relative abundances of CH4 and CO2 are very sensitive to [M/H], the C/O ratio (e.g., see Moses et al. 2013a,b), and the quench temperature, the behavior implied from the retrieved abun-dances suggests that rapid vertical transport may drive nonequilibrium abundances even for higher-temperature objects.

CONCLUSIONS
We have presented the results of a detailed retrieval analysis of the red L dwarf 2M2224-0158. After testing some 61 different combination of cloud opacities, we have found that the clouds appear to be dominated by layers of small grained amorphous enstatite (particle radii, r 0.1 µm) and quartz (r ∼ 0.4 µm) at shallow pressures ( 0.1 bar), combined with a deep iron cloud deck (r ∼ 0.1 µm) becoming optically thick at 1 µm at a pressure close to 10 bar. The cloud opacity is best modelled by a Hansen distribution of particle sizes that scatter light according to Mie theory. Our analysis strongly rejects the log-normal distribution and the Distribution of Hollow Spheres scattering model. We estimate a radius of 0.75 ± 0.02 RJup, which is considerably smaller than predicted by evolutionary models for a field age object with the luminosity of 2M2224-0158.
Our retrieved thermal profile matches the grid model predictions as well as its flexibility allows at pressures deeper than about 1 bar. However, we find a stratosphere that is some 500 K warmer than expected, suggestive of additional vertical heat transport not included in the self-consistent models.
All of our highly-ranked models assume vertically constant mixing fractions for our absorbing gases, rather than those predicted by thermochemical equilibrium. This likely reflects a combination of quenched carbon chemistry, and ongoing issues with the pressure-broadened alkali opacities.
The combination of quartz and enstatite implies either that this target has a Mg/Si ratio of significantly less than 1 (∼ 0.7), or that phase equilibrium arguments based upon traditional assumptions (e.g. solar elemental abundance ratios) do not well predict dominant cloud opacities. This lowvalue of Mg/Si is consistent with the stellar distribution, which shows significant spread. We have also estimated a high-metallicity ([M/H] = 0.38 +0.07 −0.06 ) and high-C/O ratio (0.83 +0.06 −0.07 ), both of which lie at the upper end of the stellar distribution in the Solar Neighbourhood. High-metallicity is correlated with both high C/O ratios and low Mg/Si ratios in stars, so our retrieved composition appears to be telling a self-consistent story.
Predictions of cloud compositions have historically focused on solar composition atmospheres, or have assumed solar abundance ratios when scaling for metallicity. Since the spread in Mg/Si for exoplanets can be expected to show a similar spread to that of the stellar population (e.g. Bonsor et al. 2021), our results demonstrate the need to extend both microphysical and phase equilibrium cloud models to a wider range of compositions if the full potential of upcoming JWST observations are to be realised. In addition, it will be necessary to test these models against cases where the ground-truth Mg/Si ratio can be estimated, such as in brown dwarf binary companions to main sequence F and G type stars.
Such a detailed treatment of a target is time consuming both in CPU hours and human effort. However, by exploring such a wide-range of both plausible and implausible cloud properties, we can place some confidence that our results are providing meaningful insight to the conditions in this target's atmosphere, and demonstrates our ability to test the predictions of different cloud models. It also highlights the potential of using cloud opacity to constrain abundance ratios such as Mg/Si. Future explorations using our Brewster framework will not require such a broad set of models, and extension to a larger number of LT dwarfs is already well underway. Furthermore, our demonstrated ability to effectively leverage the wide 1 − 15µm wavelength range makes Brewster especially well suited to exploiting JWST data for self-luminous exoplanets in the coming years.