# Probing the Extent of Vertical Mixing in Brown Dwarf Atmospheres with
Disequilibrium Chemistry

Sagnick Mukherjee1 , Jonathan J. Fortney1 , Natasha E. Batalha2 , Theodora
Karalidi3 , Mark S. Marley4 , Channon Visscher5 , Brittany E. Miles1 , and
Andrew J. I. Skemer1 1Department of Astronomy and Astrophysics, University of
California, Santa Cruz, CA 95064, USA
2 NASA Ames Research Center, MS 245-3, Moffett Field, CA 94035, USA
3Department of Physics, University of Central Florida, 4111 Libra Dr, Orlando,
FL 32816, USA
4 Lunar and Planetary Laboratory, The University of Arizona, Tucson, AZ 85721,
USA
5 Dordt University, Sioux Center IA; Space Science Institute, Boulder, CO, USA

###### Abstract

Evidence of disequilibrium chemistry due to vertical mixing in the atmospheres
of many T and Y-dwarfs has been inferred due to enhanced mixing ratios of CO
and reduced NH3. Atmospheric models of planets and brown dwarfs typically
parameterize this vertical mixing phenomenon with the vertical eddy diffusion
coefficient, $K_{zz}$. While $K_{zz}$ can perhaps be approximated in the
convective regions in the atmosphere with mixing length theory, in radiative
regions the strength of vertical mixing is uncertain by many orders of
magnitude. With a new grid of self-consistent 1D model atmospheres from
$T_{\rm eff}$ of 400 - 1000 K, computed with a new radiative-convective
equilibrium python code PICASO 3.0, we aim to assess how molecular abundances
and corresponding spectra can be used as a probe of depth-dependent $K_{zz}$.
At a given surface gravity, we find non-monotonic behavior in the CO abundance
as a function of $T_{\rm eff}$, as chemical abundances are sometimes quenched
in either of two potential atmospheric convective zones, or quenched in either
of two possible radiative zones. The temperature structure and chemical
quenching behavior also changes with gravity. We compare our models with
available near-infrared and M-band spectroscopy of several T and Y-dwarfs and
assess their atmospheric vertical mixing profiles. We also compare to color-
magnitude diagrams and make predictions for JWST spectra. This work yields new
constraints, and points the way to significant future gains, in determining
$K_{zz}$, a fundamental atmospheric parameter in substellar atmospheres, with
significant implications for chemistry and cloud modeling.

###### Subject headings:

Brown Dwarfs, T dwarfs, Y dwarfs, Atmospheric Composition

##  1\. Introduction

###  1.1. Molecular abundances and non-equilibrium chemistry

The astronomical community is now ready to observe the atmospheres of a
diverse range of exoplanets and brown dwarfs in unprecedented detail with the
James Webb Space Telescope (JWST) (Gardner et al., 2006; Pontoppidan et al.,
2022). These atmospheres are complex systems with a range of physical
processes like radiative/convective energy transport, chemical reaction
networks, clouds/haze formation, and dynamical processes like vertical mixing.
Upcoming high signal–to–noise (SNR) spectra from JWST for these substellar
atmospheres will provide us with a unique opportunity to probe these physical
and chemical processes and compare to state–of–the–art theoretical models.
Therefore, it is crucial to develop and update these models by including
necessary complex processes and assessing how uncertainties in physical and
chemical parameters affect the structure, abundances, and spectra of these
atmospheres.

One such complex process in substellar atmospheres is atmospheric dynamics.
Dynamics is a three-dimensional process. In the case of solar system planets,
atmospheric dynamics can be directly studied with spatially and time resolved
imaging or spectroscopy. Although, we cannot spatially resolve brown dwarfs
and exoplanets, dynamical processes can still be studied by their impact on
molecular abundances and condensate clouds, which affect their spectra and
lightcurves. Surface inhomogeneities can be probed via rotationally modulate
lightcurves (Metchev et al., 2015; Cushing et al., 2016; Karalidi et al.,
2015; Vos et al., 2022), often in two dimensions when observed at a single
wavelength. Spectra, which probe a range of atmospheric pressures, provide a
potential view of dynamics in the additional dimension of atmosphere depth.
Processes like convection, the breaking of gravity waves, and convective
overshoot can impact both the chemical and cloud structure of these
atmospheres by causing net vertical transport of atmospheric gases and
particles across many atmospheric scale heights (Freytag et al., 2010;
Parmentier et al., 2013; Bordwell et al., 2018; Tan, 2022).

Dynamical processes can also alter the abundances of atmospheric gases in the
visible atmosphere by transporting gas molecules through several scale
heights, on a timescale faster than chemical reactions are able to locally
equilibrate to their new surroundings. This leads to departures from the
chemical equilibrium, often assumed as a starting point in exoplanet and brown
dwarf atmosphere models (e.g., Marley et al., 2021). The dynamics induced
transport of gases and particles along the radial dimension of these
atmospheres is called “vertical mixing”. Vertical mixing can lead to
disequilibrium chemical abundances in substellar atmospheres for gases like
CO, CH4, H2O, NH3, CO2, PH3, and HCN (Fortney et al., 2020; Tsai et al., 2017,
2021; Moses et al., 2011). Because most of these gases are major sources of
opacities in substellar atmospheres (H2O, CH4, CO, etc.), their spectra can be
altered significantly.

Vertical mixing also directly affects cloud formation through it’s effect on
molecules like H2O in Y-dwarfs and colder exoplanets (Morley et al., 2014;
Ackerman & Marley, 2001a; Gao et al., 2020). Vertical mixing transports
condensable vapor from the deeper atmosphere to lower pressures and
temperatures, where they condense. Moreover, mixing counteracts gravitational
settling, and therefore helps to dictate the size of cloud particles that can
remain aloft. As a result, the atmospheric vertical mixing significantly
impacts the cloud opacity by controlling the cloud vertical extent and
particle size distribution. As both the chemical and cloud structure of these
atmospheres impact radiative and convective energy transport, vertical mixing
also affects the temperature structure of substellar atmospheres.

There is already decades-deep literature on observations and modeling of
disequilibrium chemistry in Jupiter and Saturn (Beer, 1975; Larson et al.,
1978; Prinn & Barshay, 1977; Larson et al., 1977; Fegley & Prinn, 1985;
Visscher et al., 2010b), and in T-dwarfs and early Y-dwarfs (Fegley & Lodders,
1996; Noll et al., 1997; Oppenheimer et al., 1998; Saumon et al., 2003;
Golimowski et al., 2004; Hubeny & Burrows, 2007; Geballe et al., 2009;
Sorahana & Yamamura, 2012; Leggett et al., 2012; Visscher et al., 2010a;
Visscher & Moses, 2011; Zahnle & Marley, 2014; Miles et al., 2020). We briefly
summarize the key developments in the field in §1.2. But, more recently, Miles
et al. (2020) obtained M-band spectra of four late T and early Y-dwarfs and
used archival M-band spectra of two T-dwarfs and a Y-dwarf to infer the
strength of vertical mixing in their atmospheres as a function of $T_{\rm
eff}$. They determined the CO abundance in each of these objects by fitting
models to spectra. 1D brown dwarf atmospheric models with equilibrium
chemistry were post-processed to include the effect of disequilibrium
chemistry due to vertical mixing to interpret the observed CO abundances and
infer the vigor of vertical mixing in these atmospheres. This demonstrated
that measurements of chemical abundances of various gases along with 1D
forward models can be used to constrain vertical mixing in substellar
atmospheres. However, disequilibrium chemistry can lead to significant changes
in the atmospheric state of these objects compared to equilibrium chemistry
calculations due to the feedback of the modified chemical abundances of gases
on the atmospheric temperature structure (Hubeny & Burrows, 2007; Karalidi et
al., 2021; Phillips et al., 2020).

The Miles et al. (2020) paper and new modeling work points to the possibility
of using upcoming high SNR JWST spectra of brown dwarfs to probe their
atmospheric structure and atmospheric chemistry in a much more detailed way.
The rationale of our study here is to specifically examine how differences in
the strength of vertical mixing, and how it is typically parametrized in
convective and radiative atmospheric regions, gives rise to distinct
differences in model spectra. Comparison with measured spectra now and in the
future can then be used to assess the mixing speed in convective and radiative
regions, yielding novel constraints on atmospheric physics.

###  1.2. Modeling Approaches

Mixing in substellar atmospheres is a 3-dimensional process. Both 2D and 3D
models have been used to study vertical mixing in hotter brown dwarfs and hot
Jupiters (Tan, 2022; Tan & Showman, 2021; Parmentier et al., 2013; Freytag et
al., 2010; Bordwell et al., 2018). However, a full exploration of the relevant
parameter space of brown dwarf mixing properties, is currently not feasible
with 3D models. Moreover, these 3D models often include simplified or gray
radiative transfer, so in detail the calculated temperature structures yield
model spectra that are not yet an adequate match to measured spectra (Tan &
Showman, 2021; Mukherjee et al., 2021). Chemical equilibrium is also often
assumed in GCMs (e.g., Lee et al., 2021; Showman et al., 2009) but some models
also treat chemical species within a passive transport framework (e.g.,
Parmentier et al., 2013; Komacek et al., 2019; Tan, 2022) to explore the
effects of atmospheric mixing on chemical structure of substellar atmospheres.

1D self-consistent models, which iterate to a radiative-convective equilibrium
(RCE) solution, given specified atmospheric physics and chemistry, include
many of the complex physical processes of real atmospheres and can be used to
produce high SNR model observables. Given the relatively fast computational
times, they are important in exploring wide ranges of parameter space (e.g.
Marley et al., 1996; Burrows et al., 1997; Marley & McKay, 1999; Burrows et
al., 2003; Fortney et al., 2008; Morley et al., 2014; Gandhi & Madhusudhan,
2017; Malik et al., 2017; Marley et al., 2021). Although such models only
include one dimension, the effects of dynamics can still be included. This is
typically done by approximating vertical mixing as a diffusive process, and
parameterizing this diffusive “speed’ with a single vertical eddy diffusion
coefficient – $K_{zz}$ (Allen et al., 1981) which has dimensions of [L]2/[T],
where [L] and [T] represent length and time, respectively. However,
theoretical and observational constraints of $K_{zz}$, and its dependence on
temperature and pressure, is uncertain by several orders of magnitude today
(Fortney et al., 2020; Miles et al., 2020; Phillips et al., 2020).

The vigor of vertical mixing could well have a complex behavior, even within
the atmosphere of a single object. Brown dwarf atmospheres have convective and
radiative zones. In the convective zones mixing is expected to be
comparatively efficient with higher values of $K_{zz}$, compared to the
radiative zones which are stable to convection. Theoretically, $K_{zz}$ has
been approximated with mixing length theory in the convective zones of these
atmospheres (e.g., Ackerman & Marley, 2001b; Smith, 1998; Gierasch & Conrath,
1985). However, this approximation is still uncertain due to uncertainties in
key parameters in the mixing length theory like the mixing length itself.
Theoretical understanding of $K_{zz}$ is even more primitive in the case of
radiative zones. Typically, $K_{zz}$ is varied by a factor of a _million_ or
more, with its effects on the temperature structure and atmospheric abundances
inadequately explored. This poor knowledge of $K_{zz}$ leaves large gaps in
our understanding of vertical mixing on chemistry, clouds, and temperature
structure.

Previously, several authors computed RCE models that treated non-equilibrium
chemical abundances self-consistently with the atmospheric temperature
structure (e.g. Hubeny & Burrows, 2007; Phillips et al., 2020; Karalidi et
al., 2021). However, these models were all computed either within the
framework of a constant $K_{zz}$ value in the atmosphere, with an exploration
of $K_{zz}$ varied across many orders of magnitude (e.g. Karalidi et al.,
2021; Phillips et al., 2020), or with $K_{zz}$ computed using mixing length
theory in the convective zones and constant $K_{zz}$ in the radiative zones
(e.g. Hubeny & Burrows, 2007). We build on these works, but aim for a more
precise assessment of how a changing $K_{zz}$ within a given atmosphere, due
to convective- and radiative-zone mixing, alters atmospheric abundances and
spectra. The output from 3D GCMs show that $K_{zz}$ can vary by several orders
of magnitude with pressure in a substellar atmosphere (Freytag et al., 2010;
Parmentier et al., 2013; Tan, 2022), but this has not yet been implemented
within a self-consistent RCE model.

Observations from Miles et al. (2020) showed that T-dwarfs with $T_{\rm eff}$
between 400–800 K have much smaller $K_{zz}$ than expected from convective
mixing. They also found a factor of 100 _increase_ in the inferred $K_{zz}$ in
these brown dwarfs when $T_{\rm eff}$ _decreased_ below $\sim$ 400 K. Miles et
al. (2020) suggested an explanation for this behavior if in the warm objects
with $T_{\rm eff}$ between 400-800 K, gases are “quenched” within a radiative
part of the atmosphere (lower $K_{zz}$), but in a convective part of the
atmosphere in objects with $T_{\rm eff}$$<400$K. Therefore, already available
observations demand models with higher complexity than a constant $K_{zz}$
approach. This is especially needed to prepare these models adequately to
interpret the much higher quality data expected from JWST.

In this work, we enhance the models presented in Karalidi et al. (2021) by
inclusion of temperature- and pressure- dependant $K_{zz}$ profiles within the
1D RCE framework. We use this model to explore the complex interconnections
between $K_{zz}$, the temperature structure, molecular abundances, and spectra
in late T and early Y-dwarfs. We use this model and already available ground-
based spectroscopic data of six brown dwarfs to address the following
questions,

  1. 1.

How does a self-consistent treatment of vertical mixing affect the $T(P)$
profile?

  2. 2.

Does a self-consistent treatment of $K_{zz}$ cause changes in the radiative
and convective zone locations and depths in brown dwarfs compared to
equilibrium chemistry models?

  3. 3.

In what objects can we expect quenching of gases like CO, CH4 and NH3 to occur
in radiative zones instead of convective zones?

  4. 4.

Is photometry and spectroscopy of brown dwarfs between 1-14 $\mu$m sensitive
to the $K_{zz}$ in the radiative and convective zone?

  5. 5.

How much can we learn about the atmospheric properties of late T-dwarfs by
applying self-consistent disequilibrium chemistry models on already available
infrared spectral data?

  6. 6.

Can JWST be used to accurately infer the vertical mixing strength in radiative
zones and convective zones of brown dwarfs?

We briefly discuss our new Python based atmospheric model in §2. We present
our results in §3 followed by discussion of our results in §4 and conclude in
§5

##  2\. Modeling Sub-stellar Atmospheres

We have adapted the Fortran based EGP sub-stellar atmospheric model to Python.
This new Python version of the EGP code is open-source and available as a part
of the widely used Python exoplanet atmospheric modeling tool PICASO (Batalha
et al., 2019), as a part of the PICASO 3.0 (Mukherjee et al., 2022) release.
The numerical and functional details of the code along with a series of
benchmarks with previously available models are presented in Mukherjee et al.
(2022). The EGP model has been used by our group for substellar atmospheres
for over two decades (Marley et al., 1996, 2012; Fortney et al., 2005, 2007,
2008; Morley et al., 2014; Marley et al., 2021; Karalidi et al., 2021).
Recently, Karalidi et al. (2021) updated the Fortran based EGP code with the
capability to include constant $K_{zz}$ and its impact on chemistry self-
consistently. We have converted this version of the code in Python and have
modified it by including the capability to treat pressure-dependent variable
$K_{zz}$ in substellar atmospheres. Here we briefly describe the model
methodology used in this work and we refer the reader to Mukherjee et al.
(2022) for a more detailed discussion of the numerical methodology of our
models.

###  2.1. Python Based implementation of EGP with Disequilibrium Chemistry

We divide the 1D model atmosphere in 90 pressure layers (91 layer boundaries
or levels) with logarithmically-spaced pressure values. From an initial guess
temperature-pressure ($T(P)$) profile, the model iterates on the $T(P)$
profile, atmospheric chemistry, and radiative/convective energy transport
using a Newton-Rhapson method, yielding temperature corrections at each
pressure level. These iterations continue until it finds the converged
atmospheric state in which radiative-convective equilibrium is satisfied
throughout the atmosphere. Self-consistent modeling of effects like
disequilibrium chemistry due to vertical mixing requires simultaneous
calculation of gas abundances and opacities with the iteration on the $T(P)$
profile instead of interpolation on a precomputed grid.

Each chemical reaction in the atmosphere proceeds with a characteristic
timescale $t_{\rm chem}$ which is a strong function of temperature, pressure,
and the chemical composition of the atmosphere. Zahnle & Marley (2014) used 1D
chemical kinetics models to parameterize $t_{\rm chem}$ expression for the
major O-, C-, H-, and N-bearing reactions in substellar atmospheres. These
reactions often are a part of a large network of chemical reactions with many
chemical species but still can be represented by simpler net reactions between
a few species of interest. The net chemical reactions between O-, C-, H- and
N-bearing gases we consider here are,

$\displaystyle{\rm CH_{4}+H_{2}O\rightleftarrows CO+3H_{2}}$
$\displaystyle{\rm CO+H_{2}O\rightleftarrows CO_{2}+H_{2}}$ $\displaystyle{\rm
2NH_{3}\rightleftarrows N_{2}+3H_{2}}$ $\displaystyle{\rm
CH_{4}+NH_{3}\rightleftarrows HCN+3H_{2}}$ $\displaystyle{\rm
2CO+N_{2}+3H_{2}\rightleftarrows 2HCN+2H_{2}O}$ $\displaystyle{\rm
CO+NH_{3}\rightleftarrows HCN+H_{2}O}$

Vertical mixing in the atmosphere also proceeds with a typical timescale,
which is associated with the eddy diffusion coefficient $K_{zz}$,

$t_{\rm mix}=\dfrac{H^{2}}{K_{\rm zz}}$ (1)

where $H$ is the local scale height of that atmospheric layer. $t_{\rm chem}$
is generally much shorter than $t_{\rm mix}$ in high pressure, high
temperature regions of the atmosphere where gas abundances are expected to
follow chemical equilibrium. But as the atmosphere gets colder at lower
pressures $t_{\rm chem}$ rises exponentially and becomes much longer than
$t_{\rm mix}$. The pressure at which the intersection between the two
timescales occurs is called the “quench pressure” or “quench level.” At all
pressure levels lower than the quench pressure (higher up the atmosphere)
atmospheric abundances for the relevant chemical compounds take on the
abundances at the relevant quench pressure, rather than those of equilibrium
chemistry at the local $P$ and $T$. In the brown dwarf context this behavior
was first predicted by Fegley & Lodders (1996) and observed by Noll et al.
(1997).

Within our model, the $T(P)$ profile is used to calculate the $K_{zz}$ as a
function of each model layer’s pressure. $K_{zz}$ in the convective zone is
calculated with mixing length theory (Gierasch & Conrath, 1985),

$K_{zz}=\dfrac{H}{3}{\left(\dfrac{L}{H}\right)}^{{4}/{3}}{\left(\dfrac{RF}{\mu{\rho_{a}}c_{p}}\right)}^{{1}/{3}}$
(2)

where $H$ is the local scale height of the atmosphere, $L$ is the turbulent
mixing length, $R$ is the universal gas constant, $\mu$ is the mean molecular
weight of the atmosphere, $\rho_{a}$ is the atmospheric density, $c_{p}$ is
the atmospheric specific heat at constant pressure and $F$ is the convective
heat flux. The convective heat flux is given by the difference between
$\sigma{T_{\rm eff}}^{4}$ and the radiative flux at each layer within the
convective zones of the atmosphere. The maximum convective heat flux would
then be $\sigma{T_{\rm eff}}^{4}$. Here the convective heat flux is calculated
self-consistently in each layer. The turbulent mixing length $L$ is set to be
equal to the scale height ($H$) in our models. But in order to explore
scenarios where convective mixing is less effective, we also simulate cases
where the mixing length is smaller than the atmospheric scale height $H$ by a
factor of 10. Note that this turbulent mixing length ($L$) was only changed in
Equation 2 and not in Equation 1. The mixing length in Equation 1 remained the
same as the local atmospheric scale height throughout this study. For further
background, we refer the reader to previous work (e.g. Visscher & Moses, 2011;
Smith, 1998; Bordwell et al., 2018; Bézard et al., 2002; Visscher et al.,
2010b) exploring the relationship between mixing length and the scale height
in Equation 1.

Over the past several years various parameterizations of $K_{zz}$ in the
radiative zones of substellar atmospheres have been discussed in the
literature (e.g., Wang et al., 2015; Zhang & Showman, 2018; Parmentier et al.,
2013). As a starting point, we use the parametrization from Moses et al.
(2021),

$K_{zz}=\dfrac{5\times 10^{8}}{\sqrt{P_{bar}}}\left(\dfrac{H_{\rm
1mbar}}{{620}\,{\rm km}}\right)\left(\dfrac{T_{\rm eff}}{1450\,{\rm
K}}\right)^{4}$ (3)

where $P_{\rm bar}$ is the pressure in bars and $H_{\rm 1mbar}$ is the
atmospheric scale height at 10-3 bars. The 5$\times$108/$\sqrt{P}$ dependence
was estimated by Parmentier et al. (2013) using tracer transport in 3D GCM for
hot Jupiter HD209458b. This dependence was found to be valid for this object
between 10-6-1 bars. The P-0.5 dependence captures the increase of $K_{zz}$
with lowering pressure in the radiative zone as mixing can be efficient in low
pressure (P $\sim$ 10-4 bars) parts of the radiative atmosphere due to
mechanisms like gravity wave breaking (Freytag et al., 2010). The factor
$H_{\rm 1mbar}T_{\rm eff}^{4}$ was also included in the parametrization as
this factor shows good correlation with the measured $K_{zz}$ (via
photochemical models) of solar system planets (including Jupiter, Saturn,
Uranus, Neptune, Earth, and Venus) and exoplanets (Moses et al., 2005; Zhang &
Showman, 2018; Moses et al., 2021). But, as the $T(P)$ structure of brown
dwarfs can be quite different from irradiated planets, we also explore
deviations from this parametrization where this relation is increased or
decreased by factors of 100 to explore the relevant parameter space. As the
$T(P)$ profile iterates towards the converged solution, the $K_{zz}$ in
convective zones, given by Equation 2 and the $K_{zz}$ in radiative zones,
given by Equation 3, evolve as well. The mixing timescale in each atmospheric
layer is then calculated in every iteration from this $K_{zz}$ using Equation
1. The quench pressures of each of these chemical species are found by
comparing $t_{\rm mix}$ with $t_{\rm chem}$ expressions from Zahnle & Marley
(2014). The abundances of CO, CH4, NH3, CO2, H2O and HCN are allowed to follow
equilibrium chemistry at pressures higher than their relevant “quench
pressure” and their abundances are quenched above the quench pressure level.
But as the $T(P)$ profile iterates to the converged solution, the quench
pressure of various gases also change at each step as both the $t_{\rm chem}$
and $t_{\rm mix}$ are dependent on the $T(P)$ profile of the atmosphere. This
requires the RCE model to include “on–the–fly” calculations of both chemistry
and opacities in order to include disequilibrium chemistry self-consistently
in the calculations.

As in Karalidi et al. (2021) we use the “on–the–fly mixing” method of
opacities following the methodology of Amundsen et al. (2017). In this work we
focus on the quenching of CO, CH4, H2O, NH3, N2, HCN and PH3. Of these, only
CO, CH4, H2O, and NH3 are major opacity sources. Therefore, we mix the
correlated-k opacities of CO, CH4, H2O, and NH3 with the correlated-k
opacities of all the other sets of gases which follow equilibrium chemistry as
the volume mixing ratio of these gases evolve due to quenching with the $T(P)$
profile of the atmosphere. As in Karalidi et al. (2021), given inaccuracies
inherent in the flexible Amundsen et al. (2017) resort-rebin approach
necessitates the use of a denser grid of 661 wavelength bins, compared to the
196 used for our equilibrium chemistry models (Marley et al., 2021).

We benchmark the code in two ways: 1) by setting $K_{zz}$= 0 and finding that
the on–the–fly mixing method still converges to the equilibrium chemistry
solutions from Marley et al. (2021), and 2) by testing the resulting models
against results obtained by Phillips et al. (2020) with non-zero but constant
$K_{zz}$ values. The benchmarking results are presented in Mukherjee et al.
(2022). The code results matches well with both models. In future work the
code will be further enhanced via coupling with 1D chemical kinetics model
VULCAN (Tsai et al., 2017, 2021) so that it can self-consistently treat both
vertical mixing and photochemistry. This Python based code will be open-
sourced as a part of the PICASO 3.0 package. Next, we discuss results obtained
by applying this code to model disequilibrium chemistry and characterize
mixing strengths in late T and early Y-dwarf atmospheres.

##  3\. Results

We performed a detailed four dimensional parameter space exploration of how
varied vigor of vertical mixing in brown dwarf atmospheres affects the
atmospheric state and spectra. We produced a grid of models by varying four
physical parameters – $T_{\rm eff}$, gravity, radiative zone $K_{zz}$, and
convective zone mixing length $L$. The gravity has been varied across five
values 316, 562, 1000, 1780 and 3160 ms-2 ($\log(g)$ from 4.5 to 5.5 with an
interval of 0.25 dex). The $T_{\rm eff}$ range is 400 K to 1000 K with an
interval of 25 K.

To explore the role of radiative zone $K_{zz}$, in our grids we multiply
Equation 3 by a multiplicative factor of $\times$ 100, $\times$ 1, and
$\times$ 0.01. We refer to these parameter values as “Moses$\times$100”,
“Moses”, and “Moses/100”, respectively, after Moses et al. (2021). The mixing
length parameter in the convective zones ($L$ in Equation 2) has been varied
between the scale height $H$ and $0.1\times{H}$ to explore scenarios where
mixing in the convective zone is less efficient than the typically assumed
mixing length theory in 1D models. Furthermore, for reasons discussed in §4,
we ran two sets of models – one where PH3 is quenched given its appropriate
$t_{\rm chem}$ and another where PH3 follows equilibrium chemistry while other
molecules are quenched. This parameter space has a total of
5$\times$25$\times$3$\times$2$\times$2$=$1500 self-consistent radiative-
convective equilibrium models. All of these models are cloud free and assume
solar metallicity and C/O ratio, following Marley et al. (2021). We haven’t
included clouds in our models as at the temperature range of our grid (400 –
1000 K), clouds are not expected to be a significant opacity source in the
upper atmosphere. We present our key findings from this grid of models below
beginning with the feedback of disequilibrium chemistry on $T(P)$ profiles of
brown dwarfs.

###  3.1. Effect of Disequilibrium Chemistry on the $T(P)$ profile of brown
dwarfs

Disequilibrium chemistry leads to enhanced CO, depleted CH4, NH3, and H2O
abundances in T and Y-dwarfs. This leads to major changes in the wavelength-
dependent optical depths in the atmosphere. The top four panels of Figure 1
shows the $\tau=1$ pressure levels, for specific gases, for a T-dwarf with
$T_{\rm eff}$ = 700 K and log(g) = 5.0 in the presence and absence of
quenching due to vertical mixing. These $\tau=1$ pressure levels are strongly
wavelength dependent. Here $K_{zz}$ has been assumed to follow Equation 2 in
the convective zones and Equation 3 in radiative zones. The disequilibrium
chemistry pressure levels are shown in solid lines and the corresponding
equilibrium chemistry levels are shown as dotted lines. The pressure-dependant
volume mixing ratios of these gases under chemical equilibrium and
disequilibrium are also shown in the top x-axis in each of the top four panels
with dashed and solid black lines, respectively.

Figure 1.— Top panels show the $\tau=1$ pressure levels for specific gaseous
absorbers in a model T-dwarf with $T_{\rm eff}$= 700 K and log(g) = 5.0 in the
presence and absence of quenching due to vertical mixing. $\tau=1$ pressure
levels for H2O, CH4, CO, and NH3 are shown in the first, second, third and
fourth columns, respectively. The equilibrium chemistry pressure levels are
shown in dashed lines whereas the disequilibrium chemistry pressure levels are
shown in solid lines. The volume mixing ratio as a function of pressure of
each gas in the presence (solid lines) and absence (dashed lines) of vertical
mixing induced quenching are also shown with black lines in the top x-axis of
each of the top four panels. Bottom panels shows the logarithm of the pressure
ratio of these $\tau=1$ levels with equilibrium chemistry and disequilibrium
chemistry for the gaseous absorbers. The zero mark represents no change in
these pressure levels. Positive values represents an increase in absorption in
the disequilibrium chemistry models compared to equilibrium chemistry models.
Main Point – Disequilibrium chemistry leads to large opacity changes in the
atmosphere.

The bottom panels in Figure 1 depicts the log of the ratio in the $\tau=1$
pressure level between these two models. Values greater than zero represent
shallower depth probed, due to an increase in an absorber of interest, for
instance CO in yellow in the third column, in particular at 4.5-5 ${\mu}m$.
There are also significant changes due to depleted NH3 and CH4. There is a
more modest effect on H2O. Overall, since major opacity sources H2O, NH3, CH4
are depleted, and they absorb across a broad wavelength range, one sees
deeper, to a higher pressure, in a disequilibrium chemistry model.

A manifestation of these effects are shown in Figure 2, which shows
temperature structure (upper left), CO mixing ratios (upper right), and
spectra (bottom) three models at $T_{\rm eff}$ = 500, 700, and 1000 K and
log(g) = 5. Solid lines include induced quenching and dashed lines represent
equilibrium chemistry. The $K_{zz}$ for these models evolved self-consistently
following mixing length theory in the convective zones and followed Equation 3
in the radiative zones.

Disequilibrium abundances mainly lead to colder $T(P)$ profiles compared to
equilibrium chemistry models as can be seen in Figure 2. In the overall
optically thinner disequilibrium chemistry atmosphere (see Figure 1) one must
see down to a higher pressure to obtain a given $T$. This effect was also seen
in Karalidi et al. (2021) for constant $K_{zz}$ models.

Volume mixing ratio profiles of CO for each of these models are shown in the
top right panel of Figure 2. As described earlier, disequilibrium chemistry
greatly enhances the CO abundances above quench level. However, even deep
down, below the quench level, there are important CO abundance differences due
to differences in the converged $T(P)$ profile. As equilibrium chemistry
models have hotter $T(P)$ profiles, at a given pressures, than disequilibrium
chemistry models, they also have higher CO abundances. This drives home the
point that a self-consistent treatment of disequilibrium chemistry (as in
§2.1) is needed rather than post-processing self-consistent equilibrium
chemistry models with quenched abundances. Even if one obtains the correct
quench pressures during post-processing equilibrium chemistry models, one
would calculate a physically inconsistent and higher abundance of CO at the
upper atmospheres than self-consistent disequilibrium chemistry models.

The thermal emission spectra of these models are shown in Figure 2 bottom
panel. The flux from these objects were normalized by their maximum flux and
are plotted from the hottest model ($T_{\rm eff}$=1000 K) at the top to the
coldest model ($T_{\rm eff}$=500 K) at the bottom. Large differences in
outgoing thermal flux can be seen in the H, K, L and M bands between the
equilibrium chemistry and disequilibrium chemistry model spectra. Comparing
these differences with the differences in the $\tau=1$ pressure levels in
Figure 1 makes it clear that the differences in the H, K, and L bands arises
due to changes in the CH4 abundance profiles. Another major difference is in
M-band. CO enhancement and modest H2O depletion leads to a major suppression
of thermal flux between 4.5-5 $\mu$m in the disequilibrium models.

Although NH3 abundances show large changes between equilibrium chemistry and
disequilibrium chemistry models, it still is not the major opacity source at
the shorter IR wavelengths. It becomes a major opacity source at 11 $\mu$m and
causes a distinct feature in the spectra. Figure 2 bottom panels show a
significant change in this NH3 feature due to disequilibrium chemistry.
Quenching of NH3 causes this feature to become slightly less prominent
compared to chemical equilibrium models.

Other carbon and nitrogen bearing molecules like HCN, N2 and CO2 are also
quenched due to $K_{zz}$. But these molecules do not have readily apparent
absorption features in the infrared at these spectral resolutions in a solar
metallicity mixture. N2 is essentially transparent in solar metallicity
mixtures and detecting CO2 is very difficult from the ground but JWST, like
AKARI (Sorahana & Yamamura, 2012), can detect disequilibrium abundances of CO2
in brown dwarf with higher metallicity. As our models have solar metallicity,
we do not specifically discuss these signatures here.

Figure 2.— Top left panel $T(P)$ profile of brown dwarf models with log(g) = 5
and $T_{\rm eff}$ of 500, 700, and 1000 K. The solid $T(P)$ profiles assume
disequilibrium chemistry self-consistently whereas the dotted profiles show
equilibrium chemistry calculations. The convective zones are shown with
thicker lines for each $T(P)$ profile. To differentiate the convective zones
of the equilibrium chemistry models from the disequilibrium chemistry models,
they are shown with much thicker lines compared to the convective zones of
disequilibrium chemistry models. Top right panel Plotted are the CO volume
mixing ratio profiles for these brown dwarf models. The solid lines show the
mixing ratio with disequilibrium chemistry whereas the dotted lines show
calculations from equilibrium chemistry. Bottom panel shows the normalized
thermal emission spectra of these brown dwarfs. The solid lines show the
disequilibrium chemistry results while dotted spectra represent the
equilibrium chemistry calculations. The J, H, K’, L’, and M’ bands are also
shown as shaded regions of different colors. Main Point– Self- consistent
treatment of disequilibrium chemistry leads to large changes in $T(P)$
profile, abundances, and thermal emission spectra. Figure 3.— Left panel shows
the $T(P)$ profile of a brown dwarf with $T_{\rm eff}$ of 700 K and log(g) of
4.5. The dashed line is the $T(P)$ profile when chemical equilibrium is
assumed whereas the solid line shows the $T(P)$ profile for disequilibrium
chemistry. The convective zones are shown with darker and thicker lines. Right
panel shows the wavelength dependant absorptivity of the different atmospheric
layers for both the cases. The absorptivity is given by 1-e-τ. The
absorptivity for the equilibrium chemistry model is shown with dashed lines
whereas that for disequilibrium chemistry model are shown with solid lines.
The local normalized blackbody spectrum is also shown for each layer with
dashed and solid lines for equilibrium and disequilibrium chemistry models,
respectively. The blue horizontal line for each layer shows the absorptivity
corresponding to ${\tau}$ =1. Main Point – Opacity changes lead to large
($\sim$100 K) differences in $T(P)$ profiles between chemical equilibrium and
disequilibrium models.

Of particular interest to our work here are the changes in $T(P)$ structure,
and the locations of radiative and convective zones show in Figure 2. We
examine the causes for these differences in more detail in Figure 3 for a
brown dwarf with $T_{\rm eff}$ = 700 K and log($g$) of 4.5. The left panel in
Figure 3 shows the $T(P)$ profile of this brown dwarf under equilibrium
chemistry with the dashed line and the $T(P)$ profile of the brown dwarf with
chemical disequilibrium is shown with the solid line. The radiative zones are
shown with narrower lines whereas the convective zones are shown with thicker
lines.

There are significant changes in the location and number of radiative and
convective zones in the two profiles, with the disequilibrium model having one
convective zone at depth and radiative zone at lower pressure. However, the
equilibrium models finds two convective zones and two radiative zones. To
examine why, the right panel in Figure 3 shows the wavelength dependent
absorptivity in several layers of these atmospheres with dashed (equilibrium
chemistry) and solid (disequilibrium chemistry) lines. The absorptivity is
given by A=(1-e-τ). A high absorptivity (close to 1) corresponds to an opaque
wavelength whereas a low absorptivity corresponds to a transparent wavelength
window. The normalized local blackbody Planck function Bλ is also plotted for
each layer for the equilibrium (dashed black line) and disequilibrium
chemistry (solid black line) models for comparison with the absorptivity. We
follow Marley & Robinson (2015) for the explanation of this phenomenon below.

As the atmosphere becomes colder with lower pressure in both models, the peak
of the local Bλ shifts to longer wavelengths. In the deep atmosphere $\sim
100$ bar, gaseous optical depths are very high which causes the absorptivity
to be very close to 1 across all wavelength ranges for both cases. The
atmosphere is opaque to radiative energy transport in these pressure layers
for both cases, and convection ensues at these pressure levels. This is the
case at pressures higher than $\sim$ 40 bars for both models. At about 20 bars
some transparent windows start to appear. Since the equilibrium chemistry
model is hotter than the disequilibrium chemistry model, the peak of its local
blackbody emission coincides more closely in wavelength with the transparent
opacity windows between 0.8-1 $\mu$m in the atmosphere at these pressures.
This causes the equilibrium chemistry model to be radiative in these pressure
ranges whereas the disequilibrium chemistry model is still convective as its
local blackbody emission still peaks at longer wavelengths where the
atmosphere is opaque at $\sim$ 20 bars pressure. As the atmosphere gets colder
with lower pressure, large transparent windows develop between 0.5-2 $\mu$m
but as these layers are much colder, their blackbody peaks are at wavelengths
greater than 2 $\mu$m.

At $\sim$ 2 bars, the blackbody peak of the chemical equilibrium model hits
another opaque range between 2-3 $\mu$m. This causes it to develop another,
detached, convective zone at these pressures. At pressures lower than this
range, a large transparent opacity window develops between 3.5-5 $\mu$m in
both the models. The disequilibrium model is more transparent than the
equilibrium model between 3.5-4.5 $\mu$m (due to less CH4 and H2O) but is more
opaque than the equilibrium model between 4.5-5 $\mu$m (due to more CO). As
the local blackbody function of both models also peaks between 3.5-5 $\mu$m in
these pressure layers, these differences in transparency between the two
atmospheres causes a large change ($\sim$ 150 K) between their $T(P)$ profiles
in these pressures. As the cooler local blackbody peak of the disequilibrium
model more closely overlaps with the opaque window in the disequilibrium model
between 4.5-5 $\mu$m, its atmosphere is still convective in these pressures,
while the equilibrium model is radiative. At pressures lower than $\sim$ 0.8
bars, both atmospheres are radiative as more transparent windows start
appearing in longer wavelengths with decreasing pressure. With these
differences in the $T(P)$ profile and convective/radiative zones between
models with chemical equilibrium and disequilibrium, the next important
question to be investigated is whether gases like CO, NH3, and H2O are
quenched in the radiative or convective zones of models with chemical
disequilibrium. Results from our model grid on this question are presented
next.

###  3.2. Quenching in the Radiative Zone or the Convective Zone ?

Quenching of gases can occur in convective zones or radiative zones. Vertical
mixing in convective zones is expected to be very efficient, yielding a high
$K_{zz}$ and correspondingly low $t_{\rm mix}$. The top three panels of Figure
4 show the $T(P)$ profile of a brown dwarf with $T_{\rm eff}$ of 700 K and
log(g) of 4.5 (left panel), the mixing and chemical timescales ($t_{\rm mix}$
and $t_{\rm chem}$) for the same brown dwarf (middle panel), and the volume
mixing ratios of various molecules (right panel) in the brown dwarf
atmosphere. This brown dwarf has a single convective and radiative zone. The
mixing timescale $t_{\rm mix}$ is shown in dashed black whereas the $t_{\rm
chem}$ for CO and other chemical reactions are shown with different colors of
dashed lines. The $t_{\rm mix}$ profile has been calculated using mixing
length theory in the convective zone and with Equation 3 in the radiative
zone. The mixing timescales in radiative zones are expected to be longer at
depth, due to the lack of efficient dynamical mixing processes. However, as
briefly described in §2, the mixing timescale is expected to be short in the
much lower pressure regions (P $\sim$ 10-4 bars) due to the presence of
gravity wave breaking mechanisms (Freytag et al., 2010; Moses et al., 2005,
2011).

$t_{\rm mix}$ is much longer than $t_{\rm chem}$ at pressures $>100$ bars and
therefore CO and CH4 are expected to follow chemical equilibrium in these high
pressure high temperature parts of the atmosphere. But the $t_{\rm chem}$ for
CO$\leftrightarrow$CH4 reaction is a very strong function of both temperature
and pressure (Zahnle & Marley, 2014). Therefore, it increases by 12 orders of
magnitudes within one order of magnitude in pressure and a change of 500 K in
temperature. As a result, the two timescales cross each other at about $\sim$
30 bars at the quench pressure ($P_{\rm Q}$). The abundances of both of these
molecules, and H2O, are quenched to constant values above $P_{Q}$ as a result,
as shown in the Figure 4 top right panel. As different net chemical reactions
proceed at different rates and have _different_ $t_{\rm chem}$, their $P_{Q}$
will also be different.

Figure 4.— Model $T(P)$ profile, chemical and dynamical timescales, and
chemistry for a brown dwarf with $T_{\rm eff}$ 700 K and log(g) = 4.5 is shown
in the top panels whereas the same properties for a brown dwarf with $T_{\rm
eff}$ 700 K and log(g) = 5.25 is shown in the bottom panels. $T(P)$ profile of
the two brown dwarfs are shown in the left panels – one with a single deep
convective zone while another with a deep ”sandwiched” radiative zone between
two deep convective zones. The middle panels show the comparison between the
mixing timescale (black) and the chemical timescales of different gases like
CO and NH3 as a function of pressure for each of the $T(P)$ profiles shown in
the left columns. The right panels show the volume mixing ratios of various
molecules as a function of pressure for each atmosphere. Main Point – Gases
can be quenched in a deeper “sandwiched” radiative zone in some brown dwarfs,
where mixing timescales can be orders of magnitude larger than in convective
zones.

Self-consistent chemical equilibrium models of T-dwarfs with $T_{\rm eff}$
within the range of 400–800 K have shown the presence of an upper detached
convective zone in addition to the deeper convective zone (Burrows et al.,
1997; Marley et al., 2021), with a radiative zone in between. The $T(P)$
profile in this scenario is shown in Figure 4 bottom left panel for a brown
dwarf with $T_{\rm eff}$ of 700 K (as above) but with a log(g) of 5.25. This
detached zone appears when the $T(P)$ profile of the atmosphere is such that
the radiative energy transport in the deeper atmosphere becomes efficient due
to the opening of a molecular opacity window at near infrared wavelengths
(Marley & Robinson, 2015), as shown in Figure 4. The resulting “sandwiched
radiative zone” is very likely to have a much slower rate of vertical mixing
(longer $t_{\rm mix}$) compared to the convective zones. As these radiative
zones appear at high pressures (P $\geq$ 10 bars), it is possible that
molecules like CO, CH4 and NH3 are quenched in the radiative zone instead of
the convective zone. Note that the mixing timescale profile in the upper
middle panel in Figure 4 shows a “kink” like structure in $t_{\rm mix}$ values
between 10-30 bars. This kink is due to the decrease in the convective heat
flux ($F$ in Equation 2) in the convective zone around that region. This can
happen in convective regions where a non-negligible, but smaller than the
convective heat transport, amount of energy transfer in a convective zone
still occurs via radiation. This leads to a smaller convective heat flux in
these regions.

The bottom middle panel of Figure 4 shows this scenario. Within our modeling
framework, the $t_{\rm mix}$ is low (high $K_{zz}$) in the deep convective
zone and $t_{\rm mix}$ is high (low $K_{zz}$) in the deep sandwiched radiative
zone. The $P_{Q}$ for CO in this case is within the sandwiched radiative zone.
Therefore, if one is trying to probe the $K_{zz}$ of this hypothetical object
with constant $K_{zz}$ models (e.g. Miles et al., 2020), one would infer a
much lower $K_{zz}$ than what is expected from Equation 2. However, the reason
for this is is not because of sluggish mixing in the convective zone but
rather quenching is occurring in a radiative zone where mixing is expected to
be sluggish.

Such a scenario was suggested in fits to observations presented in Miles et
al. (2020). They inferred the $K_{zz}$ for several late T and early Y-dwarfs
and found that the $K_{zz}$ for objects with $T_{\rm eff}$ between 400-800 K
are much smaller than expectations from mixing length theory. Although
equilibrium chemistry models find the sandwiched radiative zone between
$T_{\rm eff}$ of 400–800 K (e.g. Marley et al., 2021; Allard et al., 2012), it
is not clear if the self-consistent treatment of disequilibrium chemistry
preserves these sandwiched radiative zones and also if molecules like CO or
NH3 actually quench in these sandwiched radiative zones. Next, we investigate
this scenario in more detail.

Figure 5.— Appearance of the 2nd radiative zone in the $T_{\rm eff}$ vs.
log(g) parameter space in our grid is shown. The colors show the value of the
quantity log10(Pbot/Ptop) where Pbot and Ptop are the pressures at the bottom
and top of the sandwiched radiative zone, respectively. The top panel shows
this in case of equilibrium chemistry. The subsequent panels show the presence
of the 2nd radiative zone if different vigor of $K_{zz}$ is chosen keeping the
convective zone mixing length fixed at $H$. Dark blue indicates no second
convective zone at depth. Figure 6.— Appearance of the 2nd radiative zone in
the $T_{\rm eff}$ vs. log(g) parameter space in our grid is shown. The colors
show the value of the quantity log10(Pbot/Ptop) where Pbot and Ptop are the
pressures at the bottom and top of the sandwiched radiative zone,
respectively. The top panel shows this in case of equilibrium chemistry. The
subsequent panels show the presence of the 2nd radiative zone if different
vigor of $K_{zz}$ is chosen keeping the convective zone mixing length fixed at
$0.1\times{H}$. Dark blue indicates no second convective zone at depth.

We can expand the modeling work shown in Figure 4 over a range in $T_{\rm
eff}$ and log (g) to assess under what conditions a sandwiched radiative zone
occurs in self-consistent non-equilibrium models. Figure 5 shows the
appearance of the sandwiched radiative zones across this phase space in the
presence and absence of vertical mixing. The colors shown in Figure 5
represents log10(Pbot/Ptop) where Pbot and Ptop are the pressures at the
bottom and top of the sandwiched radiative zone, respectively. If the model
does not find a sandwiched radiative zone, this number will be zero. The top
panel in Figure 5 shows the sandwiched radiative zones in chemical equilibrium
models. The other three panels show the appearance of the second radiative
zone with three different magnitudes of mixing in the radiative zone, all
assuming the mixing length $L$ is $H$ in the convective region. Under the
assumption of equilibrium chemistry it is clear that the sandwiched radiative
zones appear between $T_{\rm eff}$ of 450-800 K for a log(g) of 4.5. However,
for objects with log(g) of 4.75 and higher, the sandwiched radiative zones
only appear between 450-600 K in chemical equilibrium models. For models with
disequilibrium chemistry, the models that show the presence of this sandwiched
radiative zone appear at a different set of $T_{\rm eff}$ and log(g) values.
As can be seen in Figure 5 (second, third, and fourth panels), high gravity
models with log(g) of 5.0 and greater have this sandwiched radiative zone
between $T_{\rm eff}$ of $\sim$ 500-800 K. For models with log(g) of 4.75 the
sandwiched radiative zone appears between 500-750 K. However, for models with
log(g) of 4.5, the sandwiched radiative zone is non-existent when quenching
due to vertical mixing is taken into account. The thickness of these zones is
typically smaller, too, as can be seen from the shading that is less yellow.

Figure 6 shows the appearance of this sandwiched radiative zone if the
convective mixing length $L$ is smaller, only $0.1\times{H}$. The second,
third, and fourth panels in Figure 6 show that larger number of models now
develop this sandwiched radiative zone compared to models with $L=H$ shown in
Figure 5. The sandwiched radiative zone does indeed appear at low gravity
(log(g) = 4.5) when $L=0.1\times H$. Additionally, the radiative zones tend to
appear at higher $T_{\rm eff}$ than for $L=H$ models. This shows that the
strength of vertical mixing in both the radiative and convective zone can
influence the presence and absence of multiple convective zones.

The effect of radiative and convective $K_{zz}$ on the chemical abundances is
more complex than it may appear. The presence of a sandwiched radiative zone
does not guarantee that quenching occurs _in_ the radiative zone. Furthermore,
radiative zone quenching can also occur in atmospheres with a single radiative
zone and single convective zones. We investigate the resulting photospheric CO
abundances in $T_{\rm eff}$ vs. log(g) space in Figure 7. The top panel shows
the CO abundance in models in chemical equilibrium. The second, third, and
fourth panels show models with disequilibrium chemistry. Each of these panels
represents models with different values of $K_{zz}$ in the radiative zone
while the mixing length in the convective zone is fixed to $H$. Models where
CO, CH4, and H2O are quenched in radiative zones are marked with black
hatches. The rest of the models still have CO quenching, but in the convective
zones instead.

The CO abundances in these models with radiative zone CO quenching also show
large differences from CO abundances in models with convective zone quenching.
These differences appear as discontinuities in the CO abundance heat map shown
in the second, third, and fourth panels of Figure 7. The reason behind these
discontinuities is mainly the smaller $P_{Q}$ of CO in models with radiative
zone quenching compared to models with convective zone quenching.

However, Figure 8 shows that CO abundances over $T_{\rm eff}$ and log(g) will
change if the mixing length in the convective zone is smaller. Here, again,
$L=0.1\times H$, leading to smaller $K_{zz}$ values, and longer $t_{\rm mix}$
in the convective zone. This causes quenching to often occur in the top or
sandwiched radiative zone in a significantly larger fraction of the $T_{\rm
eff}$ and log(g) parameter space compared to the nominal $L=H$ scenario in
Figure 7. For models with log(g)=4.75 and greater, all models with $T_{\rm
eff}$ $\geq$ 550 K show radiative zone quenching of CO and CH4 in this case.
Both Figure 7 and 8 makes it clear that the chemical properties of brown
dwarfs that will have either radiative zone or convective zone quenching
depends on the vigor of mixing in both of these types of zones.

Furthermore, since $t_{\rm chem}$, and hence the $P_{Q}$, of different gaseous
species are different, this means that models where CO quenches in radiative
zones are not necessarily models where CO2 or NH3 also quenches in the
radiative zone. This effect can be clearly seen in Figure 9, which shows the
same set of models in Figure 7, but now with the CO2 abundances and its
radiative/convective quench behavior. Because the $t_{\rm chem}$ for the
thermochemical equilibrium of CO2 is orders of magnitude shorter than for CO
and CH4, the models where CO2 quenches in the radiative zone are also
different than those where CO/CH4 quenches in the radiative zone. If the
mixing length in the convective zone is $H$, CO2 radiative zone quenching
occurs mostly in models that are hotter than 700 K and have gravity higher
than log(g) of 5. The same effect for quenching of NH3 has been shown in
Figure 10 in models with $L=0.1\times{H}$. Comparing Figure 10 with 8 shows
that the radiative zone quenching of NH3 occurs in only a small subset of
lower temperature (500–800 K) models with radiative zone quenching of CO ,
when $L=0.1\times{H}$. The $t_{\rm chem}$ for the reaction between
N${}_{2}\leftrightarrow$NH3 is the slowest among the other reactions described
in §2.1. As a result, if $L=H$, NH3 always gets quenched in the deep
convective zone in the parameter space explored by our grid. Even if the
convective mixing length is smaller ($L=0.1\times{H}$), NH3 gets quenched in
the convective zone of the majority of the models except in some where deeper
“sandwiched” radiative zones develop. But in order to measure these quenched
abundances observationally and infer the $K_{zz}$ profile of the atmosphere,
we need to investigate the effect of $K_{zz}$ on observables, such as the
emission spectra of these brown dwarfs, which we discuss next.

Figure 7.—  Top panel shows the equilibrium chemistry abundance of CO at
pressures where $T$= $T_{\rm eff}$ in our $T_{\rm eff}$ vs. gravity grid of
chemical equilibrium models. The second, third and last panels depict the
quenched CO abundance within the $T_{\rm eff}$ vs. gravity parameter space in
our disequilibrium chemistry models. The second, third and last panels
correspond models with radiative zone $K_{zz}$ values of “Moses”,
“Moses$\times$100”, and “Moses/100”, respectively while the convective zone
mixing length L was fixed at the scale height H. Models where CO was quenched
at the radiative zone are marked with black hatches. CO was quenched at a
convective zone in the rest of the models. Models with radiative zone CO
quenching show much smaller CO abundance compared to models with convective
zone quenching. Figure 8.— Same as Figure 7 except with the mixing length in
convective zones fixed at 0.1$\times$H. In these models CO is mostly quenched
in a radiative zone. Figure 9.— Same as Figure 7 but for quenching of CO2.
Given the different $t_{\rm chem}$ for CO2, compared to CO, the quench levels
are entirely different. Figure 10.— Same as Figure 8 but for quenching of NH3.
Given the different $t_{\rm chem}$ for NH3, compared to CO, the quench levels
are entirely different.

###  3.3. Effect of Radiative and Convective Zone $K_{zz}$ on Emission Spectra

Figure 11.— M-band Thermal emission spectra across a range of $T_{\rm eff}$
between 400 K and 1000 K with steps of 100 K are shown in each panel. log
($g$) = 5.0 is on the left and 5.5 is on the right. The colors are
representative of the $T_{\rm eff}$. The models shown in the upper left and
right panels assume the convective mixing length $L$ to be $H$ whereas the
models shown in the lower left and right panels assume $L=0.1\times{H}$. The
three line styles show different radiative zone $K_{zz}$ parametrizations,
including Moses (solid), Moses/100 (dashed), and Moses$\times$100 (dash-dot).
If quenching occurs in a convective zone the three spectra plot atop each
other. The shaded region shows the wavelength range where CO is the most
important absorber. Figure 12.— H-band thermal emission spectra of models from
our grid are shown in the top left and right panels whereas the bottom left
and right panels show the K-band emission spectra. Models with log(g)=5 and
different $T_{\rm eff}$ between 400-1000 K are shown in each panel. From
bottom to top, spectra for $T_{\rm eff}$ values of 425 K, 700 K, 850 K, 900 K,
and 1000 K are shown in all the panels. Top left panel shows the thermal
emission spectra of models which have convective mixing length set at $H$ but
have different $K_{zz}$ parametrizations like Moses, Moses/100, and
Moses$\times$100 which are represented by different line types in the panels
like solid, dashed, and dashed-dotted. Top right panel shows the same but for
models with convective mixing length set at $0.1\times{H}$. Bottom left and
right panels show the sensitivity of the K band spectrum to radiative zone
$K_{zz}$ when $L=H$ and $L=0.1\times{H}$, respectively. The shaded region
shows the wavelength range where CH4 is the most important absorber. Figure
13.— Thermal emission spectra of brown dwarfs between 3 and 3.5 $\mu$m are
shown in the top left and right panels. Emission spectra of the same models in
the L-band are shown in the bottom left and right panels. Models with log(g)=5
and different $T_{\rm eff}$ between 400-1000 K are shown in each panel. From
bottom to top, spectra for $T_{\rm eff}$ values of 425 K, 700 K, 850 K, 900 K,
and 1000 K are shown in all the panels. Top left panel shows the emission
spectra of models which have convective mixing length set at $H$ but have
different $K_{zz}$ parametrizations like Moses, Moses/100, and
Moses$\times$100 which are represented by different line types in the panels
like solid, dashed, and dashed-dotted. Top right panel shows the same but for
models with convective mixing length set at $0.1\times{H}$. Bottom left and
right panels show the sensitivity of the L-band spectrum to radiative zone
$K_{zz}$ when $L=H$ and $L=0.1\times{H}$, respectively. The shaded region
shows the wavelength range where CH4 is the most important absorber.

We have established in §3.2 that gases like CO, CH4 and CO2 can be quenched in
the sandwiched radiative zones or single radiative zones in a subset of late T
and early L-dwarfs. In Figures 11, 12, and 13, we show how thermal emission
spectra are sensitive to the vigor of vertical mixing in the radiative zone of
these objects. These spectra, with a spectral resolution of 600, are
calculated using the thermal emission calculation module of PICASO (Batalha et
al., 2019), using the same opacities that were used for model convergences to
RCE. Sources of these opacities can be found in Marley et al. (2021).

The thermal emission spectra between 4.2-5.0 $\mu$m, across a range of $T_{\rm
eff}$ between 400 K and 1000 K with steps of 100 K, are shown for log(g)
values of 5 and 5.5 in the two upper panels of Figure 11. The choice of the
wavelength range focuses on the role of CO, which has significant impacts on
the spectrum in M-band between 4.4-5 $\mu$m. H2O is the primary absorber
between 4.0-4.4 $\mu$m. These models assume $L=H$ in the convective zones.
Figure 11 lower panels shows the spectrum from the same set of models except
$L=0.1\times{H}$. Spectra arising from different vigors of vertical mixing in
the radiative zones are shown with different line dashes in all the panels.

Sensitivity of the M-band spectrum to radiative zone $K_{zz}$ values is more
common for higher gravity objects with log(g) of 5.0 and 5.5. If the
convective mixing length is $H$, then this sensitivity appears between models
with $T_{\rm eff}$ 550- 700 K for log(g) of 5. However, if the convective
mixing length is $0.1\times{H}$, then this sensitivity appears much more
commonly between models with $T_{\rm eff}$ of 500-1000 K for log(g)= 5.0. At
higher gravity, log(g) of 5.5, the sensitivity of M-band spectra to radiative
zone $K_{zz}$ appears at both convective mixing lengths of $H$ and
$0.1\times{H}$. In Figure 11 upper right panel, this sensitivity appears
between 500-900 K whereas in Figure 11 lower right panel this sensitivity
appears between 500-1000 K.

While the M-band shows the greatest sensitivity to $K_{zz}$ in the radiative
and convective zones, smaller differences in the spectra arise in other
wavelengths as well. Figure 12 top panels shows the sensitivity of the H-band
spectra to the radiative zone $K_{zz}$ across various values of $T_{\rm eff}$
for brown dwarf models with log(g)=5.0. The top left panel shows models with
$L=H$ whereas the top right panel shows models with $L=0.1\times{H}$. The
sensitivity mainly appears between 1.69 $\mu$m and 1.75 $\mu$m in models where
CH4 is quenched in the radiative zone as it is the main opacity source at
these wavelengths. CH4 is also the main opacity source between 2–2.2 $\mu$m,
3.2-3.5 $\mu$m, and 3.5–4 $\mu$m. As a result, the spectra at these wavelength
ranges also show some sensitivity to the radiative zone $K_{zz}$ of the deep
atmosphere. This sensitivity between 2–2.2 $\mu$m, 3.2-3.5 $\mu$m, and 3.5–4
$\mu$m is shown in Figure 12 bottom panels and Figure 13 top and bottom
panels, respectively. Examining different spectral regions for a given object,
that are sensitive to both CH4 and CO, will help in constraining the radiative
zone $K_{zz}$, as both gas species abundances can be used as probes. Apart
from the emission spectra itself, other observables like broadband IR colors
and magnitudes of brown dwarfs are also known to be sensitive to $K_{zz}$
(e.g. Saumon et al., 2003; Karalidi et al., 2021), which we explore next.

###  3.4. Color Magnitude Diagrams

Figure 14.— Color magnitude diagrams (CMD) of observed brown dwarfs in
different filters are shown with overplotted equilibrium and disequilibrium
chemistry model tracks with different strengths of convective and radiative
zone vertical mixing. The top row shows the J vs. J-H where the models with
convective zone mixing length of $H$ and different radiative zone $K_{zz}$ are
shown with various linetypes of red. Models with convective zone mixing length
of $0.1\times{H}$ and different radiative zone $K_{zz}$ are shown with various
linetypes of blue and the equilibrium chemistry models are shown with dashed
black lines. Each column shows the same data but the model tracks with
different values of log(g). The second row shows the the L${}^{{}^{\prime}}$
vs L${}^{{}^{\prime}}$-M${}^{{}^{\prime}}$ whereas the third row shows the W1
vs. W1-W2.

Coupling our atmospheric model grid with brown dwarf evolutionary models from
Marley et al. (2021) enables us to generate color magnitude diagrams (CMDs) of
these brown dwarfs. These evolutionary models track the $T_{\rm eff}$, log(g)
and bolometric luminosity as a function of age for each brown dwarf mass. The
procedure is the same as in Marley et al. (2021), although our models here are
slightly not self-consistent, as the upper boundary condition for the Marley
et al. (2021) cooling tracks assume equilibrium chemical abundances.

Figure 14 shows the observed CMD of brown dwarfs in various standard infrared
filter combinations. We show J vs. J-H, MKO L${}^{{}^{\prime}}$ vs.
L${}^{{}^{\prime}}$-M${}^{{}^{\prime}}$, and WISE W1 vs. W1-W2 in the first,
second and third rows, respectively. The mid-infrared wavelengths are most
sensitive to disequilibrium chemistry effects. Model tracks for three
different log(g) (4.5, 5.0 and 5.5) values are shown in the three columns for
each combination of filters. The database maintained by Trent Dupuy
111http://www.as.utexas.edu/ tdupuy/plx (Dupuy & Liu, 2012; Dupuy & Kraus,
2013; Liu et al., 2016) has been used for the observational data. The red
model tracks in all the panels represents models with convective mixing length
set at $H$ whereas the blue tracks represents models with convective mixing
length set at $0.1\times{H}$. Each line type (solid, dotted, and dash-dotted)
in both of these colors represents models with the different radiative zone
$K_{zz}$ values “Moses$\times$100”, “Moses”, and “Moses/100,” respectively.
The black dashed lines represents CMD tracks from equilibrium chemistry
models.

In the J vs. J-H filter combination CMD, it appears that models treating
disequilibrium chemistry self-consistently provide a slightly better fit to
the observed photometry than the equilibrium chemistry, for models at log(g)
of 5.0-5.5, representative of the field sample. The disequilibrium chemistry
models with different convective mixing lengths shown by the blue and red
tracks almost fall on top of each other as the J and H band fluxes are not
particularly sensitive to CO abundance but has some sensitivity to CH4
abundance.

The chemical equilibrium model CMD tracks follow a much redder sequence in the
MKO L${}^{{}^{\prime}}$ vs. L${}^{{}^{\prime}}$-M${}^{{}^{\prime}}$ CMDs
compared to the observed photometry. The disequilibrium chemistry models are
comparatively much bluer. As discussed in the previous sections, enhanced CO
abundance in disequilibrium chemistry model causes the M′-band and the W2-band
fluxes to be fainter than equilibrium chemistry models. On the other hand,
depleted CH4 abundance in disequilibrium chemistry models causes higher fluxes
in the L′ and W1 bands. Both of these effects causes the disequilibrium model
tracks (shown with red and blue lines) to be significantly bluer than the
equilibrium chemistry model tracks (shown with black line) in the
L${}^{{}^{\prime}}$ vs. L${}^{{}^{\prime}}$-M${}^{{}^{\prime}}$ and W1 vs.
W1-W2 CMD.

It is interesting to note the non-monotonic behavior. At higher gravity
(log(g) = 5.0 and 5.5 ) some disequilibrium chemistry models show radiative
zone quenching instead of convective zone quenching. As a result, these
atmospheres have lower photospheric CO compared to disequilibrium chemistry
with convective quenching, but a much higher CO abundance compared to the
equilibrium chemistry models. Therefore, the models with radiative zone
quenching are redder than disequilibrium chemistry models with convective zone
quenching but bluer than equilibrium chemistry models, “jutting out” to the
red for some models.

The second and third rows in Figure 14 also show that disequilibrium models
with $L=H$ with convective zone quenching of CO are much bluer than the
observed sequence. However, disequilibrium models with $L=0.1\times{H}$
provide a generally better fit to the mid-IR photometry. This suggests that
longer $t_{\rm mix}$ values in the deeper atmosphere provide a better fit to
the photometry from the T and Y dwarf population, than the $t_{\rm mix}$
calculated from assumptions like free convection. However, this result is
empirical in nature and needs further validation from 3D tracer transport
models like in Bordwell et al. (2018). Coupling our models with 1D chemical
kinetics models will also help us in interpreting these trends better, which
is something we will address in future studies. It seems clear that even lower
$K_{zz}$ values for the faintest models may be necessary, as the data are
generally more consistent with equilibrium chemistry models at log (g) of
5.0-5.5 We also note that we have not included clouds within our models and
they might have some impact on the model magnitudes presented here especially
for the faintest end of these tracks as H2O clouds are expected to form in the
colder objects ($T_{\rm eff}$ $\leq 450$ K) (Morley et al., 2014).

###  3.5. Comparing to Spectroscopic Observations

We use our grid of atmosphere models to compare with ground-based spectral
observations of six late T and early Y-dwarfs. We aim to determine the
atmospheric properties of brown dwarfs – Gleise 570 D (Burgasser et al.,
2000), 2MASS J0415-0935 (Burgasser et al., 2002), WISE 0313 (Kirkpatrick et
al., 2011), UGPS 0722 (Cushing et al., 2011), WISE 2056 (Schneider et al.,
2015), and WISE 1541 (Schneider et al., 2015) with our models. Although Miles
et al. (2020) had the Y-dwarf WISE 0855 in their sample, we exclude this
object from our analysis as it is expected to have optically thick water
clouds in its atmosphere due to its cold $T_{\rm eff}$ ($\sim$ 250 K) and our
model grid only extends to a $T_{\rm eff}$ = 400 K. As noted in §3, this
choice for the lower bound on $T_{\rm eff}$ in our model grid was made because
our models lack a self-consistent treatment of clouds with disequilibrium
chemistry, currently.

For comparing our models with observational data, which includes both near-
infrared and M-band spectra, we use the goodness of fit parameter – Gk defined
as (Cushing et al., 2008),

$G_{k}=\dfrac{1}{n-m}\sum_{i=1}^{n}w_{i}\left(\dfrac{f_{i}-C_{k}{F^{k}_{i}}}{{\sigma}_{i}}\right)^{2}$
(4)

where $n$ is the number of observational data points, $m$ is the number of
free parameters being fit, $w_{i}$ are the weights of fitting $i$’th data
point, $f_{i}$ is the $i$’th observed flux, $F^{k}_{i}$ is the $i$’th model
flux point of the $k$’th model, and $\sigma_{i}$ is the error in the $i$’th
observed flux. The factor $C_{k}$ is a renormalization factor for the model
fluxes. If the observed fluxes are flux calibrated then $C_{k}$ represents the
factor $R^{2}/D^{2}$ but if the observed fluxes are not flux calibrated then
$C_{k}$ represents the factor $\alpha{R}^{2}/D^{2}$ where 1/$\alpha$
represents the unknown flux normalization of the observed flux. The $C_{k}$
factor is given by (Cushing et al., 2008),

$C_{k}=\dfrac{\sum_{i=1}^{n}w_{i}f_{i}F^{k}_{i}/{\sigma}_{i}^{2}}{\sum_{i=1}^{n}w_{i}(F^{k}_{i})^{2}/{\sigma}_{i}^{2}}$
(5)

Since there is no reported bias in any of the datasets we use, $w_{i}$ is
always fixed at 1 for our analysis.

####  3.5.1 Gliese 570 D and 2MASS J0415-0935

Gliese 570 D is a T7.5 dwarf (Burgasser et al., 2006) which has a measured
spectrum in the M-band with two different instruments – AKARI (Sorahana &
Yamamura, 2012) and GEMINI/NIRI (Geballe et al., 2009). We used the M-band
spectra of Gliese 570 D reported in (Geballe et al., 2009) because it has a
higher SNR ratio. We also use the NIRC 1-2.5 $\mu$m spectra of Gliese 570 D
(Gl 570D) reported in Burgasser et al. (2002) and Spitzer Infrared
Spectrograph (IRS) 5-14 $\mu$m spectra of Gliese 570 D reported in Suárez &
Metchev (2022). We adopt the methodology described in Sorahana & Yamamura
(2012) for fitting the spectrum of Gl 570D with our grid of models. Firstly,
we use our grid of models to calculate $C_{k}$ using only the M-band
GEMINI/NIRI data for Gl 570D. Using this set of calculated $C_{k}$, we
calculate the $G_{k}$ parameter using our models and the M-band GEMINI/NIRI
data for Gl 570D. Then we select all the models with values of the $G_{k}$
parameter between $min(G_{k}$) and $min(G_{k}$) +1. This step gave us 131
models which best-fits the M-band GEMINI/NIRI data for Gl 570D. The $G_{k}$
parameter for these 131 models are distributed between 0.9899 and 1.9550.

Figure 15.— The top left panel shows the observed M-band spectra of Gl 570D
with red points, the 1-2.5 $\mu$m observations of Gl 570D with black points
and the Spitzer IRS measurements with brown points.. The flux values sampled
by the top 131 models with the least Gk parameter after fitting the M-band
spectra only are shown with the shaded blue region. The top right panel shows
the best-fit model obtained by fitting the M-band spectra, the 1-2.5 $\mu$m
data, and the Spitzer 5-14 $\mu$m data for Gl 570D. The bottom left panel
shows a comparison of the best-fit model with the M-band spectra only. The
bottom right panels shows the $T(P)$ profile of the best fit model, the
$t_{\rm mix}$ and $t_{\rm chem}$ of different quenched species as a function
of pressure, and the rightmost panel shows the volume mixing ratios of CO,
CH4, NH3, H2O, CO2 and PH3. Figure 16.— The top left panel shows the observed
M-band spectra of 2M 0415 with red points, the 1-2.5 $\mu$m observations of 2M
0415 with black points and the Spitzer IRS measurements with brown points. The
flux values sampled by the top 101 models with the least Gk parameter after
fitting the M-band spectra only are shown with the shaded blue region. The top
right panel shows the best-fit model obtained by fitting the M-band spectra,
the 1-2.5 $\mu$m data and the Spitzer 5-14 $\mu$m data for 2M 0415. The bottom
left panel shows a comparison of the best-fit model with the M-band spectra
only. The bottom right panels shows the $T(P)$ profile of the best fit model,
the $t_{\rm mix}$ and $t_{\rm chem}$ of different quenched species as a
function of pressure, and the rightmost panel shows the volume mixing ratios
of CO, CH4, NH3, H2O, CO2 and PH3.

The top left panel of Figure 15 shows the observed M-band spectrum of Gl 570 D
(Geballe et al., 2009) with red points, the observed NIRC spectrum of Gl 570 D
(Burgasser et al., 2002) with black points, and the Spitzer IRS spectrum
(Suárez & Metchev, 2022) with brown points. The range of fluxes covered by the
131 best-fitting models to the M-band spectra only, are shown with the blue
shaded region in the top left panel of Figure 15. All of these models fit the
M-band spectrum well but differ significantly in their fluxes in the shorter
and longer wavelengths. We calculate a new set of $G_{k}$ parameter for these
131 models but now with the NIRC and Spitzer IRS data only. The set of $C_{k}$
parameters for each model, determined by fitting the M-band spectra already,
are re-used for this, following the methodology of Sorahana & Yamamura (2012).
The model among these 131 models that fits this shorter and longer wavelength
data as well (has the least $G_{k}$) is taken to be the overall best-fit model
from the grid. This model spectra is shown in Figure 15 top right panel and it
has a Gk value of 669.06. This large difference between the Gk values from
M-band data and the NIRC+IRS data is mainly because of the large difference in
the uncertainties in the fluxes from these two sets of observations.

The best-fit model has a $T_{\rm eff}$ of 725 K and a log(g) of 5.0. The
$K_{zz}$ in the radiative zone for this model is Moses and the convective
mixing length is $0.1\times{H}$. The M-band spectra of this model compared to
the M-band data for Gl 570D is also shown in the bottom left panel. The $T(P)$
profile, mixing timescale and chemical timescales, as well as the mixing ratio
profiles of various gases in this model are shown in the three bottom right
panels. The $T(P)$ profile and the $t_{\rm mix}$ profile both show that this
model has two convective zones. The $t_{\rm chem}$ and $t_{\rm mix}$ shown in
the middle panel of the bottom right part of Figure 15 also show that CO, CH4,
H2O, HCN, and NH3 are quenched in the sandwiched radiative zone for this best-
fit model whereas CO2 is quenched in the upper detached convective zone. The
$C_{k}$ value of this best-fit model is 1.72$\times$10-19. The measured
parallax of Gl 570D is 171.22 mas (Dupuy & Liu, 2012). This corresponds to a
distance (D) of 5.84 pc. We get a radius of 1.04 RJ for Gl 570D combining
these values together. The radius of a brown dwarf with $T_{\rm eff}$=725 K
and log(g)=5.0 can also be interpolated from evolutionary grids like the
SONORA BOBCAT models (Marley et al., 2021). We find this interpolated radius
to be 0.89 RJ which is somewhat smaller than our calculated radius. But as we
also note in §3.4, these evolutionary tracks were computed assuming chemical
equilibrium.

We use the same technique for fitting the AKARI M-band spectrum (Sorahana &
Yamamura, 2012) and the NIRC + Spitzer IRS spectrum (Burgasser et al., 2002;
Suárez & Metchev, 2022) of 2MASS J0415-0935 (2M 0415). The top left panel of
Figure 16 shows the M-band spectra of 2M 0415 with red points, the NIRC
measurements are shown with black points, and the Spitzer IRS measurements are
shown with brown points. The $min(G_{k}$) obtained by fitting the M-band
spectra only is 0.8755 and the flux range covered by all the models with
$G_{k}$ between 0.8755 and 1.8583 is shown with the blue shaded region in the
top left panel. These models were then used to obtain the best-fitting model
with the low and high wavelength data. The best-fit model has a $G_{k}$ of
514.28. This model has a $T_{\rm eff}$ of 675 K and a log(g) of 4.75. The
$K_{zz}$ in the radiative zone for this model is Moses and the convective
mixing length is $0.1\times{H}$. The measured parallax of 2M 0415 is 175.2 mas
(Dupuy & Liu, 2012) which corresponds to a distance of 5.707 pc. Combining
this with the best-fit model $C_{k}$ value (1.84$\times$10-19) provides a
radius of 1.05 RJ for 2M 0415. The interpolated radius from the SONORA BOBCAT
grid for this object is 0.98 RJ. This shows much better agreement with the
calculated radius of the object compared to Gl 570D. The M-band spectrum of
this best-fit model compared to the M-band data for 2M 0415 is also shown in
the bottom left panel. Atmospheric properties of the best-fit model are shown
in the three bottom right panels. This best-fit model has a second sandwiched
radiative zone and CO, CH4, H2O, and HCN are quenched in this sandwiched
radiative zone whereas CO2 is quenched in the upper detached convective zone.

####  3.5.2 WISE 0313, UGPS 0722, WISE 2056 and WISE 1541

Figure 17.— The left panel shows a heat map of the $\chi^{2}$ value calculated
using model photometry and observed photometry (Miles et al., 2020) in various
infrared bands as a function of model $T_{\rm eff}$ and log(g) for WISE 0313.
The left panel indicates models where radiative $K_{zz}$ is fixed at Moses and
the convective mixing length is ${H}$. The middle and right panels also show
the $\chi^{2}$ heat maps using models with radiative $K_{zz}$ of Moses/100 and
Moses$\times$100, respectively. The best-fit $T_{\rm eff}$ is near 550 K.

Miles et al. (2020) has obtained M-band observations of the late T and early
Y-dwarfs WISE 0313, UGPS 0722, WISE 2056, and WISE 1541. However, these M-band
spectra are not flux calibrated. Using photometry collected from various
studies and the SONORA evolutionary tracks, Miles et al. (2020) found ranges
of $T_{\rm eff}$ and log(g) values for these objects. We use our models along
with these photometry datasets (see Table 2 in Miles et al. (2020) for
photometry obtained from Skrutskie et al. (2006); Kirkpatrick et al. (2011);
Lucas et al. (2010); Leggett et al. (2013); Schneider et al. (2015); Geballe
et al. (2001); Knapp et al. (2004); Leggett et al. (2015); Golimowski et al.
(2004); Leggett et al. (2017); Patten et al. (2006); Kirkpatrick et al.
(2012); Cutri et al. (2021); Luhman & Esplin (2016); Schneider et al. (2016);
Wright et al. (2014); Esplin et al. (2016)) to find a similar range of $T_{\rm
eff}$ values for each of these sources. Figure 17 shows the map of
${\chi}^{2}$ values in a $T_{\rm eff}$ vs. log(g) parameter space obtained by
comparing our models with observed photometry of WISE 0313. It is clear that
the ${\chi}^{2}$ values are sensitive to the $T_{\rm eff}$ of models but there
is not much sensitivity to the log(g) of the models. This method helps us in
identifying a range ($\delta$T = 100 K) of $T_{\rm eff}$ values for each
object. Then using models within this $T_{\rm eff}$ range we try to determine
the best-fit atmospheric model by comparing our model spectra shapes with the
shape of the M-band spectra of these four objects. We use the same metric
$G_{k}$ and $C_{k}$ defined in Equations 4 and 5 for finding the best-fit
models. However, as the observed spectra lacks flux calibration, the $C_{k}$
now includes the $\alpha$ parameter as well where 1/$\alpha$ represents the
unknown flux calibration of the observed spectrum. For each object, we use the
models within the determined $T_{\rm eff}$ ranges and allow all the other
parameters such as log(g), radiative zone $K_{zz}$ and mixing length $L$ to
vary while calculating the $G_{k}$ and $C_{k}$ values of the models using the
M-band spectrum. Then the model with the least $G_{k}$ value is chosen as the
best-fit model for the object.

Figure 18.— The best-fit atmospheric models and their comparison with data are
shown in this figure. Each row corresponds to a different object. The left
most column shows the best-fit M-band spectra along with the observed data
from Miles et al. (2020). The second left column shows the best-fit $T(P)$
profile for these objects from the grid of models presented in this paper. The
third column shows the best-fit mixing timescale in blue along with the
chemical timescale for CO in red as a function of pressure. The volume mixing
ratios of various atmospheric gases are shown in the fourth column for each
object.

Figure 18 shows the results of this fitting exercise for all the four objects.
The four rows in Figure 18 shows the results for WISE 0313, UGPS 0722, WISE
2056, and WISE 1541, respectively. The left columns show the best-fit spectra
compared with the observational data for each object. The second column from
left shows the best-fit $T(P)$ profile for each object. The third column shows
the best-fit mixing timescale and chemical timescale for CO as a function of
pressure for these objects and the fourth columns depicts the best-fit volume
mixing ratios for each object.

We find that the best-fit model describing the WISE 0313 and UGPS 0722 data
have radiative zone quenching of CO/CH4 whereas the best-fit models describing
the WISE 2056 and WISE 1541 data has CO quenching in the convective zone.
Table 1 shows the values of all the parameters relevant to fitting the
spectral data for all of these objects. We note that the best-fit log(g)
values obtained in this work are not from fitting any spectral shape/feature
directly. As can be seen in Figure 5, log(g) controls the presence/absence of
a sandwiched radiative zone which in turn dictates whether gases can be
quenched in regions with low or high $K_{zz}$. These effects of log(g) on
atmospheric chemistry lead to preference of certain log(g) values over others
while fitting the data.

Source |  $T_{\rm eff}$ Range[K] |  Best-Fit $T_{\rm eff}$[K] |  Best-Fit log(g) |  Radiative $K_{zz}$ |  Convective L |  M-band Gk |  J–H–K+_Spitzer_ IRS Gk |  Ck
---|---|---|---|---|---|---|---|---
Gl 570D |  700–800 |  725 |  5.0 |  Moses |  $0.1\times{H}$ |  1.414 |  669.06 |  1.72$\times$10-19
2M 0415 |  600–700 |  675 |  4.75 |  Moses |  $0.1\times{H}$ |  0.925 |  514.28 |  1.84$\times$10-19
WISE 0313 |  575–675 |  675 |  4.75 |  Moses/100 |  $0.1\times{H}$ |  3.849 |  – |  6.477
UGPS 0722 |  500–600 |  575 |  4.75 |  Moses$\times$100 |  $0.1\times{H}$ |  23.884 |  – |  11.636
WISE 2056 |  450–550 |  525 |  4.75 |  – |  $0.1\times{H}$ |  6.651 |  – |  14.575
WISE 1541 |  400–450 |  450 |  4.75 |  – |  $0.1\times{H}$ |  5.243 |  – |  20.275
Table 1Results from fitting observed spectra of brown dwarfs with our model
spectra. The radiative zone $K_{zz}$ for models which do not quench at
radiative zones are marked with a ’–’ as they couldn’t be constrained.

We choose our best-fit model to be the model with the least $G_{k}$ value but
we also find that there are quite a few models which have $G_{k}$ values only
differ from the minimum only by a small amount. This shows the presence of
degeneracy in our models with the current precision of available data for
these objects. This degeneracy will be broken with much higher precision flux
calibrated M-band and shorter wavelength spectra expected from JWST which we
discuss in §3.6.

####  3.5.3 Trends in Quench Timescales

Figure 19.— The best fit quench $K_{zz}$ of CO for each object is shown as a
function of their best-fit $T_{\rm eff}$ with the red crosses. The median
best-fit quench $K_{zz}$ vs. $T_{\rm eff}$ inferred by Miles et al. (2020)
with lower gravity models are shown with the faded red crosses. The inferred
quench $K_{zz}$ for NH3, CO2, HCN, and PH3 from the best-fit models for each
object are also shown as a function of their best-fit $T_{\rm eff}$. The black
points show the theoretically maximum $K_{zz}$ allowed from mixing length
theory for these best-fit models.

We present the trend of the best-fit $K_{zz}$ at $P_{Q}$ for CO in these
sources in Figure 19. As we mainly fit the M-band spectra of these objects in
our analysis, we are only sensitive to CO and H2O abundances and therefore we
can only fit the quench $K_{zz}$ of CO. Our best-fit quench $K_{zz}$ for CO as
a function of the best-fit $T_{\rm eff}$ for each object are shown with red
crosses in Figure 19. The best-fit quench $K_{zz}$ determined by Miles et al.
(2020) are shown with faded red crosses as a function of the best-fit $T_{\rm
eff}$ obtained by Miles et al. (2020) for each object. Although the overall
trend of much lower CO quench $K_{zz}$ than the theoretical maximum in objects
with $T_{\rm eff}$ between 600–800 K, found by Miles et al. (2020), still
stands in this analysis, there are some differences in the best-fit $K_{zz}$
and $T_{\rm eff}$ values between the two studies. These differences are due to
multiple factors which we discuss briefly here.

Our work uses RCE models that include disequilibrium chemistry self-
consistently to fit the data while equilibrium chemistry RCE models were post-
processed by Miles et al. (2020) to include quenching of gases to fit the
data. Figure 2 top left panel shows that including disequilibrium chemistry
self-consistently can lead to large changes in the $T(P)$ profile of brown
dwarfs compared to equilibrium chemistry models. This can be the reason behind
the $\sim$ 50 K differences in the $T_{\rm eff}$ values for these objects
between the two studies. For 5 of the 6 analyzed objects in this work, our
assessment of CO quench $K_{zz}$ are somewhat smaller than previously found.
As some of our best-fit models have higher $T_{\rm eff}$ than the best-fit
models of Miles et al. (2020), CO is more abundant in the deeper atmosphere of
our best-fit models than those of Miles et al. (2020). Figure 2 shows that CO
abundance falls very rapidly with decreasing pressure in the deep atmosphere
(which is under chemical equilibrium) in these colder brown dwarfs. Therefore
colder best-fit models in Miles et al. (2020) need to quench CO deeper in the
atmosphere than our hotter best-fit models. So, the best-fit models presented
in Miles et al. (2020) need higher CO quench $K_{zz}$ compared to our best-fit
models to produce the same CO abundance in the upper well mixed atmosphere of
these brown dwarfs. Additionally, the best-fit $K_{zz}$ presented in Miles et
al. (2020) has been shown to be dependant on the gravity of the model used for
fitting the data. In this work, we allowed gravity to be a free parameter
while fitting the data and compared our best-fit $K_{zz}$ in Figure 19 with
the median $K_{zz}$ obtained by Miles et al. (2020) by trying out various
values of log(g).

We can use our best-fit models with CO quenching to infer the quench $K_{zz}$
of other quenched gases like CO2, NH3 or HCN in these best-fit models. The
inferred quench $K_{zz}$ values for the other gases in these objects are also
shown in Figure 19. As each gas has a different $t_{\rm chem}$, they also are
quenched at different depths of the brown dwarf atmosphere as shown in Figure
7 and 9 for CO and CO2. This means that measuring the abundance of each gas
can help us in measuring the height dependant $K_{zz}$ profile in these brown
dwarf atmospheres. For example, in Figure 19, the inferred quench $K_{zz}$ for
CO2 are very different than the inferred quench $K_{zz}$ of CO (or NH3) as CO2
is quenched at different pressures (having different mixing strengths)
compared to CO. But inferring the height dependant $K_{zz}$ profile of a
single brown dwarf by using different gases as probes of $K_{zz}$ at different
pressures is only possible with self-consistent disequilibrium chemistry
models with height dependant $K_{zz}$ profiles such as the model grid
presented in this work.

The theoretical maximum limit on convective zone $K_{zz}$ for each best-fit
model is also shown with black points in Figure 19. This maximum limit on
convective $K_{zz}$ was calculated using Equation 2 with the convective heat
flux fixed at $\sigma$$T_{\rm eff}$4 and the mixing length $L$ fixed at the
scale height of the atmosphere. It is clear from Figure 19 that all the gases
in all these objects quench at lower $K_{zz}$ than this theoretical limit. But
even though all the objects analyzed in this work have much smaller CO quench
$K_{zz}$ than the maximum convective limit shown with the black line, this
does not mean that all these “quenchings” occurred in the radiative zones in
our models. For Gl 570D, 2M 0415, W0313, and UGPS 0722, CO is quenched in the
radiative zones in our best-fit models but for W 2056 and W 1541 CO along with
all the other molecules are quenched in the convective zone with a much
smaller mixing length than the scale height of the atmosphere as was also
suggested by the CMD analysis detailed in §3.4. Finally, we note that our
best-fit models, especially for objects which lack flux-calibrated spectra,
are often degenerate with the present data, and better data from observatories
like JWST will be instrumental in giving us better constraints on the $K_{zz}$
of these objects, which we discuss next.

###  3.6. Constraining Radiative Zone $K_{zz}$ with JWST

Figure 20.— Simulated data from JWST instruments are shown for one of our
model brown dwarfs with $T_{\rm eff}$ of 700 K. The three colors in the
simulated data correspond to models with different radiative zone $K_{zz}$.
The top left panel shows the simulated NIRSpec spectra of a log (g)=5.5 object
with an exposure time of 10 minutes while the bottom left panel shows a zoomed
version of the top left panel to emphasize how JWST NIRSpec M-band
observations can clearly distinguish between different $K_{zz}$ in the
radiative zones for these T-dwarfs. The top right panel shows the NIRSpec
simulated data for a brown dwarf with $T_{\rm eff}$ of 700 K and log(g) of
4.5. The two colors represent different $K_{zz}$ in the convective zone but
the same $K_{zz}$ in the radiative zone. The bottom right panel shows MIRI LRS
simulated data for the same models shown in the top right panels. To compare
the differences between the simulated data from different models to the
typical JWST noise, the ratio of the flux differences ($\delta$F) between
various models and the simulated JWST noise ($N$) is shown below each of these
panels. This shows that both convective and radiative zone $K_{zz}$ can be
constrained with JWST.

An important takeaway from Miles et al. (2020) and our work is that detailed
analysis of atmospheric abundances can constrain $K_{zz}$ in brown dwarf
atmospheres. There is strong evidence for mixing orders of magnitude weaker
than expectations from free convection with $L=H$. The launch of JWST begs the
question of whether JWST data will be able to give us better constraints on
the vigor of vertical mixing in these atmospheres. We can provide a test by
simulating the JWST spectra using the JWST Exposure Time Calculator
(Pontoppidan et al., 2016) with some of the emission spectra from our model
grid. As an illustrative example, we use the spectrum of a 700 K object. In
the left panels of Figure 20 we examine three models at log(g) of 5.5, where
the radiative zone $K_{zz}$ is “Moses$\times$100”, “Moses”, and “Moses/100”
and $L=H$. We assume that the source has a radius of 1RJ and is at a distance
of 5 pc. The temperature and distance chosen for this exercise are somewhat
representative of objects like Gl 570D and 2M 0415.

The NIRSpec panels (which is 3 of the 4 figures) assume fixed slit mode with
the G395M grating and the F290LP filter for these simulations. The results are
impressive, with only a total of 10.39 minutes exposure time, simulated with
10 groups per integration, 2 integration per exposure, and 5 exposures. It is
clear from the top left panel that the radiative zone $K_{zz}$ can be
constrained with high precision with about a total $\sim$ 10 minutes exposure
time per object in this JWST observing mode. This time includes all the 5
exposures used for this simulation. In order to emphasize the difference
arising in the synthetic counts from the three models in the 4.5-5 $\mu$m
range, the bottom left panel shows a zoomed in view of the 4.3-5 $\mu$m
wavelength range, focusing on CO absorption, where the main sensitivity of the
data to radiative zone $K_{zz}$ lies. The top right panel shows a lower-
gravity model (log (g)=4.5), where a lower pressure photosphere favors CO at
the expense of CH4. Here the blue model is simulated from a model where the
convective mixing length $L=H$ whereas the salmon colored data is simulated
from a model where the convective mixing length is $L=0.1\times{H}$. The
radiative zone $K_{zz}$ of both the cases in this panel are
“Moses$\times$100”. Both of these models have the quenching of all relevant
gases occurring in the convective zone. The large differences in the synthetic
counts between the two models in the 3.5–4.25 $\mu$m range shows that JWST
data can also be used to probe differences in the convective zone $K_{zz}$
with exposure times of $\sim$ 10 minutes.

Radiative and convective zone $K_{zz}$ can be probed with other wavelengths as
well. The bottom right panel in Figure 20 shows the simulated spectrum of the
two cases shown in the top right panel but when observed by the MIRI
instrument. The MIRI LRS mode was used for this data simulation along with a
60 minutes exposure time. Differences between the two model spectra arise
between 5–9 $\mu$m. A comparison with Figure 1 shows that these differences
mainly arise from different CH4, H2O and NH3 abundances due to differences in
the convective zone $K_{zz}$.

##  4\. Discussion

###  4.1. PH3 Quenching

Figure 21.— Effect of PH3 quenching on the spectrum has been shown by
comparing them with disequilibrium chemistry model spectra where PH3 is not
quenched. Models with different $T_{\rm eff}$ values ranging from 400 K to
1000 K, with a step of 100 K, are shwon from bottom to top. All of these
models have log(g)=5. Figure 22.— Top panel shows the equilibrium chemistry
abundance of PH3 at pressures where $T(P)$= $T_{\rm eff}$ in our $T_{\rm eff}$
vs. gravity grid of chemical equilibrium models. The second, third and last
panels depict PH3 abundance within the $T_{\rm eff}$ vs. gravity parameter
space in our disequilibrium chemistry models. The second, third and last
panels correspond models with radiative zone $K_{zz}$ values of “Moses”,
“Moses$\times$100”, and “Moses/100”, respectively while the convective zone
mixing length $L$ was fixed at the scale height $H$. Models where PH3 was
quenched in a radiative zone are marked with black hatches. PH3 was quenched
at a convective zone in the rest of the models.

Under chemical equilibrium conditions the major phosphorus-bearing molecule
PH3 is expected to be removed by oxidation at low temperatures in the visible
atmospheres of giant planets and colder brown dwarfs. A number of oxidation
products have been considered (see Visscher, 2020), including P4O6 (e.g.,
Fegley & Lodders, 1994; Visscher et al., 2006), P4O10 (Borunov et al., 1995),
H3PO4 (Wang et al., 2016), or condensation of NH4H2PO4 (Fegley & Lodders,
1994; Visscher et al., 2006; Morley et al., 2018). Following the work of
Visscher & Fegley (2005); Visscher et al. (2006); Miles et al. (2020), as a
representative example of PH3 oxidation we may consider the conversion of PH3
into P4O6 via the net thermochemical reaction,

${\rm P_{4}O_{6}+12H_{2}\rightleftarrows 4PH_{3}+6H_{2}O}$ (6)

However, the $t_{\rm chem}$ of PH3$\rightarrow$P4O6 is such that it can be
quenched in the deeper atmosphere, leading to a large enhancement of PH3 in
the photosphere. This effect makes PH3 the major phosphorus bearing gas in
atmospheres of Jupiter and Saturn (Lodders, 2004; Visscher & Fegley, 2005).
PH3 is expected to be found in colder brown dwarfs as well (Visscher et al.,
2006; Morley et al., 2018).

Due to expected small affect of P-bearing molecules on brown dwarf temperature
structures, we did not include PH3 as a molecule to be mixed on–the–fly for
converging to a final atmospheric state. However, in order to test the effect
of PH3 on the emergent thermal emission spectra, we post-process our converged
disequilibrium chemistry model with quenched PH3 abundance. The $t_{\rm chem}$
for PH3$\rightarrow$P4O6 is given by Visscher & Fegley (2005),

$t_{\rm chem}=\dfrac{1.9\times 10^{12}}{[\rm OH]}{\rm
exp}\left(\dfrac{6013.6}{T}\right){\rm s}$ (7)

where [OH] represents the number density of OH molecules in units of cm-3.
With our converged disequilibrium chemistry modules, we calculate the mixing
ratio of OH molecules using Visscher & Fegley (2005),

$X_{\rm OH}=\dfrac{KX_{\rm H_{2}O}\sqrt{X_{\rm H_{2}}}}{\sqrt{P}}$ (8)

where XOH, X${}_{\rm H_{2}O}$, and X${}_{\rm H_{2}}$ represent the volume
mixing ratios of each of these components and $P$ represents the pressure in
bars. The constant $K$ is temperature dependant and based on fits to
thermodynamic data, and is given by,

$K=3.672-\dfrac{14791}{T}$ (9)

where $T$ is the temperature in K. Using Equations 7, 8, and 9 we determine
the $t_{\rm chem}$ of PH3 for all our converged models. Then using our
converged $t_{\rm mix}$ profiles we determine $P_{Q}$ for PH3 and quench it
creating another set of our models which has quenched PH3 abundances. Using
these modified chemical compositions, we create another set of thermal
emission spectra models with quenched PH3 abundances. Figure 21 compares the
differences in spectra caused by the quenching of PH3 with disequilibrium
models where PH3 is not quenched. This quenching of PH3 increases the PH3
volume mixing ratios dramatically, from $\leq$ 10-30 under chemical
equilibrium, to $\sim$ 10-6 to 10-7 due to quenching. The main impact of this
increase appears between 4-4.5 $\mu$m in the thermal emission spectrum as
shown in Figure 21. This also suggests that PH3 may become an important enough
atmospheric opacity source to have a small effect on their $T(P)$ profiles as
well, which we will aim to explore in future work.

Figure 22 shows the parameter space of models where PH3 quenches at the
radiative zone and convective zones for different levels of radiative zone
$K_{zz}$ when the convective mixing length is $H$. This also shows that
depending on the properties of the object in question, PH3 can also be used as
a probe for determining radiative or convective zone $K_{zz}$. The region
around 4.0-4.4 $\mu$m is difficult to probe from the ground (Miles et al.,
2020), but it is clear from Figure 21 that these PH3 should be readily
detectable with _JWST_. Interestingly, Morley et al. (2018) didn’t detect PH3
between 4.0-4.4 $\mu$m in WISE 0855 even though it is expected to be present
in it’s atmosphere theoretically. Determining the strength of these features,
or their absence, will have important implications in understanding $K_{zz}$
and phosphorus chemistry.

###  4.2. Clouds

The models presented in this work are cloud-free. However, H2O condensation in
the upper atmospheres of objects can begin at $T_{\rm eff}$ $\leq$ 450 K
(Morley et al., 2014) and these clouds can start becoming an important opacity
source colder than $\sim$ 400 K (Morley et al., 2014). Once water clouds
become optically thick, they can change the overall energy balance of the
atmosphere impacting the $T(P)$ profile. The presence of disequilibrium
chemistry causes a cooling down of the upper atmosphere, compared to
equilibrium chemistry, but the opacity of clouds has the opposite effect of
heating up the $T(P)$ profile.

Importantly, the cloud particle size and vertical extent are also quite
sensitive to the mixing parameter $K_{zz}$ as mixing helps to keep cloud
particles aloft against their gravitational settling (Ackerman & Marley,
2001b). Cloud particle sizes are an extremely important parameter in dictating
their scattering properties (Mukherjee et al., 2021; Ackerman & Marley, 2001b)
and hence $K_{zz}$ has a large effect on the cloud opacities in these
atmospheres. Therefore, it is clear that information gleaned on $K_{zz}$ from
disequilibrium chemistry studies will be important in improving cloud models.
Future studies should aim for a self-consistent treatment of disequilibrium
chemistry and cloud formation with a self-consistent treatment of $K_{zz}$
within a given model atmosphere. Recently, Mang et al. (2022) showed that
water clouds have a large impact on the M-band spectra of colder Y-dwarfs,
with $T_{\rm eff}$ $\leq$ 400 K, just below the lower $T$ edge of our study.
Somewhat different than this work, Mang et al. (2022) calculated the radiative
zone $K_{zz}$, needed for cloud modeling, by setting the mixing length to
one–tenth of the scale height in the radiative zones. This approach also leads
to smaller $K_{zz}$ in the radiative zones compared to the convective zones.
When we extend our disequilibrium chemistry model to colder objects, inclusion
of clouds in this framework is an important milestone.

##  5\. Summary and Conclusion

In this work, we have used a newly developed Python–based 1D radiative-
convective equilibrium model with the capability to self-consistently capture
disequilibrium chemistry in brown dwarf and exoplanet atmospheres. Using this
model we have created a self-consistent grid of brown dwarf models between
$T_{\rm eff}$ of 400-1000 K with an increment of 25 K and log(g) of 4.5-5.5
with an increment of 0.25. Instead of using a constant $K_{zz}$ prescription
throughout the entire atmosphere as done in earlier work (e.g. Karalidi et
al., 2021; Phillips et al., 2020) or self-consistent $K_{zz}$ in the
convective zones and constant $K_{zz}$ in the radiative zones (e.g. Hubeny &
Burrows, 2007), we use a self-consistent $T(P)$ structure dependant $K_{zz}$
profile for the entire atmosphere for our grid. In the convective zone, we use
mixing length theory with two different mixing lengths whereas in the
radiative zones we use various multiples of a $T$ and $P$-dependent $K_{zz}$
prescription. The idea behind these models was to assess whether
disequilibrium chemical abundance can be used as a probe for temperature
structure and mixing processes below the visible atmosphere in brown dwarfs.
With this grid in hand, we conclude on the following points.

  1. 1.

As found by others, a self-consistent treatment of disequilibrium chemistry
causes cooling of $T(P)$ profiles, especially in the upper atmospheres of late
T and early Y-dwarfs.

  2. 2.

Post-processing equilibrium chemistry models by quenching various molecules
instead of using self-consistent models with disequilibrium chemistry should
be avoided because this can lead to significant errors in the calculated
abundances of molecules like CO, CH4 and NH3. This can also lead to incorrect
$T(P)$ profiles. Both of these results are especially relevant in the JWST era
with significantly higher S/N mid-infrared data.

  3. 3.

Self-consistent models including disequilibrium chemistry also show a large
change in the location, depths, and number of radiative and convective zones
in brown dwarfs compared to self-consistent models with equilibrium chemistry.

  4. 4.

Depending on the strengths of radiative zone $K_{zz}$ and convective mixing
length, gases like CO, CH4, NH3, CO2 and PH3 can be quenched in the radiative
zones of brown dwarfs in a large fraction of the $T_{\rm eff}$ and log(g)
parameter space.

  5. 5.

Since different molecules have different $t_{\rm chem}$ values, which leads
them to quench at different pressure levels. This means that some molecules
may quench in radiative zones, while other may quench in convective zones, in
the same atmosphere, such that different molecules may be quenched by very
different $K_{zz}$ values.

  6. 6.

The M-band thermal spectra of objects where radiative zone quenching occurs,
is very sensitive to the $K_{zz}$ in the radiative zones. This can lead to
non-monotonic behavior in the CO abundance and spectra as a function of
$T_{\rm eff}$. This is especially exciting because it provides an avenue to
observationally constrain the highly uncertain radiative zone $K_{zz}$ in
these objects.

  7. 7.

The mid-infrared photometric bands – MKO L${}^{{}^{\prime}}$ and MKO M, and
WISE 1 and 2, are quite sensitive to the $K_{zz}$ in the radiative zone. We
also find that our models provide a better fit to the infrared photometry of
these objects than chemical equilibrium models.

  8. 8.

We use our models to fit the M-band spectra of 6 brown dwarfs – Gleise 570 D,
2MASS J0415-0935, WISE 0313, UGPS 0722, WISE 2056, and WISE 1541. We determine
the best-fit atmospheric state for these objects including the quench level
$K_{zz}$ of quenched gases in these objects. Similar to Miles et al. (2020),
we also find that these objects have $K_{zz}$ much smaller than the
theoretical convective mixing limit. We also see that the best fit quench
$K_{zz}$ in these objects are very low between $T_{\rm eff}$ of 675-800 K
while it shows an increase with lowering $T_{\rm eff}$ for objects with
$T_{\rm eff}$ below 575 K.

  9. 9.

By performing some signal–to–noise calculations for JWST, we conclude that
short exposure time observations can provide us with excellent measurements of
the radiative and convective zone $K_{zz}$ in brown dwarfs.

  10. 10.

We also expect JWST to find clearly detectable disequilibrium PH3 abundances
in the atmospheres of these objects.

The non-monotonic behaviour of atmospheric chemistry, as a function of $T_{\rm
eff}$ and log(g), found in this work indicates the need to explore these
effects with even finer grid spacing than this work in parameters like $T_{\rm
eff}$ (25 K in this work), log(g) (0.25 dex in this work), and $K_{zz}$
(factor of 100 in this work). This remains the focus of a future model grid
release. The model presented in this work will be updated to include self-
consistent clouds in the near future. Additional future work will couple the
code to the chemical kinetics model VULCAN (Tsai et al., 2017, 2021) which
would enable a more rigorous self-consistent treatment of disequilibrium
chemistry in our models. This will also allow us to include effects like the
feedback of the disequilibrium abundances of gases like CH4 on the entire
chemical network of the atmosphere within our models, self-consistently. This
will also allow us to explore the empirical trends we find in this work, like
probable slower convective mixing than typically calculated assuming free
convection, more robustly. Moreover, such a coupling opens up the possibility
to explore both photochemistry and vertical mixing together self-consistently,
which is more relevant for irradiated exoplanets. This updated model will then
be used to explore effects of vertical mixing and photochemistry for warm
transiting giant planets and find additional observables to constrain $K_{zz}$
in such planets, an important comparison sample to the non-irradiated objects
explored here.

##  6\. Acknowledgments

SM thanks the UC Regents Fellowship award for supporting him for this work.
JJF acknowledges the support of NASA XRP grant 80NSSC19K0446. We acknowledge
use of the lux supercomputer at UC Santa Cruz, funded by NSF MRI grant AST
1828315. SM also thanks Caroline Morley for very useful inputs for preparing
the color-magnitude diagrams presented here. We thank the anonymous referee
for very insightful comments which helped in the betterment of this
manuscript.

Software: PICASO 3.0 (Mukherjee et al., 2022), PICASO (Batalha et al., 2019),
pandas (McKinney, 2010), NumPy (Walt et al., 2011), IPython (Pérez & Granger,
2007), Jupyter (Kluyver et al., 2016), matplotlib (Hunter, 2007), the model
grid will be formally released via Zenodo.

## References

  * Ackerman & Marley (2001a) Ackerman, A. S., & Marley, M. S. 2001a, ApJ, 556, 872
  * Ackerman & Marley (2001b) —. 2001b, ApJ, 556, 872
  * Allard et al. (2012) Allard, F., Homeier, D., & Freytag, B. 2012, Philosophical Transactions of the Royal Society of London Series A, 370, 2765
  * Allen et al. (1981) Allen, M., Yung, Y. L., & Waters, J. W. 1981, Journal of Geophysical Research: Space Physics, 86, 3617
  * Amundsen et al. (2017) Amundsen, D. S., Tremblin, P., Manners, J., Baraffe, I., & Mayne, N. J. 2017, A&A, 598, A97
  * Batalha et al. (2019) Batalha, N. E., Marley, M. S., Lewis, N. K., & Fortney, J. J. 2019, ApJ, 878, 70
  * Beer (1975) Beer, R. 1975, ApJ, 200, L167
  * Bordwell et al. (2018) Bordwell, B., Brown, B. P., & Oishi, J. S. 2018, ApJ, 854, 8
  * Borunov et al. (1995) Borunov, S., Dorofeeva, V., Khodakovsky, I., et al. 1995, Icarus, 113, 460
  * Burgasser et al. (2006) Burgasser, A. J., Geballe, T. R., Leggett, S. K., Kirkpatrick, J. D., & Golimowski, D. A. 2006, ApJ, 637, 1067
  * Burgasser et al. (2000) Burgasser, A. J., Kirkpatrick, J. D., Cutri, R. M., et al. 2000, ApJ, 531, L57
  * Burgasser et al. (2002) Burgasser, A. J., Kirkpatrick, J. D., Brown, M. E., et al. 2002, ApJ, 564, 421
  * Burrows et al. (2003) Burrows, A., Sudarsky, D., & Lunine, J. I. 2003, The Astrophysical Journal, 596, 587
  * Burrows et al. (1997) Burrows, A., Marley, M., Hubbard, W. B., et al. 1997, ApJ, 491, 856
  * Bézard et al. (2002) Bézard, B., Lellouch, E., Strobel, D., Maillard, J.-P., & Drossart, P. 2002, Icarus, 159, 95
  * Cushing et al. (2008) Cushing, M. C., Marley, M. S., Saumon, D., et al. 2008, ApJ, 678, 1372
  * Cushing et al. (2011) Cushing, M. C., Kirkpatrick, J. D., Gelino, C. R., et al. 2011, ApJ, 743, 50
  * Cushing et al. (2016) Cushing, M. C., Hardegree-Ullman, K. K., Trucks, J. L., et al. 2016, ApJ, 823, 152
  * Cutri et al. (2021) Cutri, R. M., Wright, E. L., Conrow, T., et al. 2021, VizieR Online Data Catalog, II/328
  * Dupuy & Kraus (2013) Dupuy, T. J., & Kraus, A. L. 2013, Science, 341, 1492
  * Dupuy & Liu (2012) Dupuy, T. J., & Liu, M. C. 2012, ApJS, 201, 19
  * Esplin et al. (2016) Esplin, T. L., Luhman, K. L., Cushing, M. C., et al. 2016, ApJ, 832, 58
  * Fegley & Lodders (1994) Fegley, Jr., B., & Lodders, K. 1994, Icarus, 110, 117
  * Fegley & Prinn (1985) Fegley, B., J., & Prinn, R. G. 1985, ApJ, 299, 1067
  * Fegley & Lodders (1996) Fegley, Bruce, J., & Lodders, K. 1996, ApJ, 472, L37
  * Fortney et al. (2008) Fortney, J. J., Lodders, K., Marley, M. S., & Freedman, R. S. 2008, ApJ, 678, 1419
  * Fortney et al. (2007) Fortney, J. J., Marley, M. S., & Barnes, J. W. 2007, ApJ, 659, 1661
  * Fortney et al. (2005) Fortney, J. J., Marley, M. S., Lodders, K., Saumon, D., & Freedman, R. 2005, ApJ, 627, L69
  * Fortney et al. (2020) Fortney, J. J., Visscher, C., Marley, M. S., et al. 2020, AJ, 160, 288
  * Freytag et al. (2010) Freytag, B., Allard, F., Ludwig, H. G., Homeier, D., & Steffen, M. 2010, A&A, 513, A19
  * Gandhi & Madhusudhan (2017) Gandhi, S., & Madhusudhan, N. 2017, MNRAS, 472, 2334
  * Gao et al. (2020) Gao, P., Thorngren, D. P., Lee, G. K. H., et al. 2020, Nature Astronomy, doi:10.1038/s41550-020-1114-3
  * Gardner et al. (2006) Gardner, J. P., Mather, J. C., Clampin, M., et al. 2006, Space Sci. Rev., 123, 485
  * Geballe et al. (2009) Geballe, T. R., Saumon, D., Golimowski, D. A., et al. 2009, ApJ, 695, 844
  * Geballe et al. (2001) Geballe, T. R., Saumon, D., Leggett, S. K., et al. 2001, ApJ, 556, 373
  * Gierasch & Conrath (1985) Gierasch, P. J., & Conrath, B. J. 1985, in Recent Advances in Planetary Meteorology, ed. G. E. Hunt, 121–146
  * Golimowski et al. (2004) Golimowski, D. A., Leggett, S. K., Marley, M. S., et al. 2004, AJ, 127, 3516
  * Hubeny & Burrows (2007) Hubeny, I., & Burrows, A. 2007, ApJ, 659, 1458
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science & Engineering, 9, 90
  * Karalidi et al. (2015) Karalidi, T., Apai, D., Schneider, G., Hanson, J. R., & Pasachoff, J. M. 2015, ApJ, 814, 65
  * Karalidi et al. (2021) Karalidi, T., Marley, M., Fortney, J. J., et al. 2021, ApJ, 923, 269
  * Kirkpatrick et al. (2011) Kirkpatrick, J. D., Cushing, M. C., Gelino, C. R., et al. 2011, ApJS, 197, 19
  * Kirkpatrick et al. (2012) Kirkpatrick, J. D., Gelino, C. R., Cushing, M. C., et al. 2012, ApJ, 753, 156
  * Kluyver et al. (2016) Kluyver, T., Ragan-Kelley, B., Pérez, F., et al. 2016, in ELPUB, 87–90
  * Knapp et al. (2004) Knapp, G. R., Leggett, S. K., Fan, X., et al. 2004, AJ, 127, 3553
  * Komacek et al. (2019) Komacek, T. D., Showman, A. P., & Parmentier, V. 2019, The Astrophysical Journal, 881, 152
  * Larson et al. (1978) Larson, H. P., Fink, U., & Treffers, R. C. 1978, ApJ, 219, 1084
  * Larson et al. (1977) Larson, H. P., Treffers, R. R., & Fink, U. 1977, ApJ, 211, 972
  * Lee et al. (2021) Lee, E. K. H., Parmentier, V., Hammond, M., et al. 2021, Monthly Notices of the Royal Astronomical Society, 506, 2695
  * Leggett et al. (2015) Leggett, S. K., Morley, C. V., Marley, M. S., & Saumon, D. 2015, ApJ, 799, 37
  * Leggett et al. (2013) Leggett, S. K., Morley, C. V., Marley, M. S., et al. 2013, ApJ, 763, 130
  * Leggett et al. (2017) Leggett, S. K., Tremblin, P., Esplin, T. L., Luhman, K. L., & Morley, C. V. 2017, ApJ, 842, 118
  * Leggett et al. (2012) Leggett, S. K., Saumon, D., Marley, M. S., et al. 2012, ApJ, 748, 74
  * Liu et al. (2016) Liu, M. C., Dupuy, T. J., & Allers, K. N. 2016, ApJ, 833, 96
  * Lodders (2004) Lodders, K. 2004, ApJ, 611, 587
  * Lucas et al. (2010) Lucas, P. W., Tinney, C. G., Burningham, B., et al. 2010, MNRAS, 408, L56
  * Luhman & Esplin (2016) Luhman, K. L., & Esplin, T. L. 2016, AJ, 152, 78
  * Malik et al. (2017) Malik, M., Grosheintz, L., Mendonça, J. M., et al. 2017, AJ, 153, 56
  * Mang et al. (2022) Mang, J., Gao, P., Hood, C. E., et al. 2022, arXiv e-prints, arXiv:2202.01355
  * Marley et al. (2021) Marley, M., Saumon, D., Morley, C., et al. 2021, Sonora Bobcat: cloud-free, substellar atmosphere models, spectra, photometry, evolution, and chemistry, doi:10.5281/zenodo.5063476
  * Marley & McKay (1999) Marley, M. S., & McKay, C. P. 1999, Icarus, 138, 268
  * Marley & Robinson (2015) Marley, M. S., & Robinson, T. D. 2015, ARA&A, 53, 279
  * Marley et al. (2012) Marley, M. S., Saumon, D., Cushing, M., et al. 2012, The Astrophysical Journal, 754, 135
  * Marley et al. (1996) Marley, M. S., Saumon, D., Guillot, T., et al. 1996, Science, 272, 1919
  * Marley et al. (2021) Marley, M. S., Saumon, D., Visscher, C., et al. 2021, ApJ, 920, 85
  * McKinney (2010) McKinney, W. 2010, in Proceedings of the 9th Python in Science Conference, ed. S. van der Walt & J. Millman, 51 – 56
  * Metchev et al. (2015) Metchev, S. A., Heinze, A., Apai, D., et al. 2015, ApJ, 799, 154
  * Miles et al. (2020) Miles, B. E., Skemer, A. J. I., Morley, C. V., et al. 2020, AJ, 160, 63
  * Morley et al. (2014) Morley, C. V., Marley, M. S., Fortney, J. J., et al. 2014, ApJ, 787, 78
  * Morley et al. (2018) Morley, C. V., Skemer, A. J., Allers, K. N., et al. 2018, ApJ, 858, 97
  * Moses et al. (2005) Moses, J. I., Fouchet, T., Bézard, B., et al. 2005, Journal of Geophysical Research (Planets), 110, E08001
  * Moses et al. (2021) Moses, J. I., Tremblin, P., Venot, O., & Miguel, Y. 2021, Experimental Astronomy, arXiv:2103.07023
  * Moses et al. (2011) Moses, J. I., Visscher, C., Fortney, J. J., et al. 2011, ApJ, 737, 15
  * Mukherjee et al. (2022) Mukherjee, S., Batalha, N. E., Fortney, J. J., & Marley, M. S. 2022, arXiv e-prints, arXiv:2208.07836
  * Mukherjee et al. (2021) Mukherjee, S., Fortney, J. J., Jensen-Clem, R., et al. 2021, ApJ, 923, 113
  * Noll et al. (1997) Noll, K. S., Geballe, T. R., & Marley, M. S. 1997, ApJ, 489, L87
  * Oppenheimer et al. (1998) Oppenheimer, B. R., Kulkarni, S. R., Matthews, K., & van Kerkwijk, M. H. 1998, ApJ, 502, 932
  * Parmentier et al. (2013) Parmentier, V., Showman, A. P., & Lian, Y. 2013, A&A, 558, A91
  * Patten et al. (2006) Patten, B. M., Stauffer, J. R., Burrows, A., et al. 2006, ApJ, 651, 502
  * Pérez & Granger (2007) Pérez, F., & Granger, B. E. 2007, Computing in Science & Engineering, 9
  * Phillips et al. (2020) Phillips, M. W., Tremblin, P., Baraffe, I., et al. 2020, A&A, 637, A38
  * Pontoppidan et al. (2022) Pontoppidan, K., Blome, C., Braun, H., et al. 2022, arXiv e-prints, arXiv:2207.13067
  * Pontoppidan et al. (2016) Pontoppidan, K. M., Pickering, T. E., Laidler, V. G., et al. 2016, in Observatory Operations: Strategies, Processes, and Systems VI, ed. A. B. Peck, R. L. Seaman, & C. R. Benn, Vol. 9910, International Society for Optics and Photonics (SPIE), 381 – 395
  * Prinn & Barshay (1977) Prinn, R. G., & Barshay, S. S. 1977, Science, 198, 1031
  * Saumon et al. (2003) Saumon, D., Marley, M. S., & Lodders, K. 2003, arXiv e-prints, astro
  * Schneider et al. (2016) Schneider, A. C., Cushing, M. C., Kirkpatrick, J. D., & Gelino, C. R. 2016, ApJ, 823, L35
  * Schneider et al. (2015) Schneider, A. C., Cushing, M. C., Kirkpatrick, J. D., et al. 2015, ApJ, 804, 92
  * Showman et al. (2009) Showman, A. P., Fortney, J. J., Lian, Y., et al. 2009, ApJ, 699, 564
  * Skrutskie et al. (2006) Skrutskie, M. F., Cutri, R. M., Stiening, R., et al. 2006, AJ, 131, 1163
  * Smith (1998) Smith, M. D. 1998, Icarus, 132, 176
  * Sorahana & Yamamura (2012) Sorahana, S., & Yamamura, I. 2012, The Astrophysical Journal, 760, 151
  * Suárez & Metchev (2022) Suárez, G., & Metchev, S. 2022, arXiv e-prints, arXiv:2205.00168
  * Tan (2022) Tan, X. 2022, Monthly Notices of the Royal Astronomical Society, 511, 4861
  * Tan & Showman (2021) Tan, X., & Showman, A. P. 2021, MNRAS, 502, 2198
  * Tsai et al. (2017) Tsai, S.-M., Lyons, J. R., Grosheintz, L., et al. 2017, ApJS, 228, 20
  * Tsai et al. (2021) Tsai, S.-M., Malik, M., Kitzmann, D., et al. 2021, ApJ, 923, 264
  * Visscher (2020) Visscher, C. 2020, Journal of Geophysical Research (Planets), 125, e06526
  * Visscher & Fegley (2005) Visscher, C., & Fegley, Bruce, J. 2005, ApJ, 623, 1221
  * Visscher et al. (2010a) Visscher, C., Lodders, K., & Fegley, Jr., B. 2010a, ApJ, 716, 1060
  * Visscher et al. (2006) Visscher, C., Lodders, K., & Fegley, Bruce, J. 2006, ApJ, 648, 1181
  * Visscher & Moses (2011) Visscher, C., & Moses, J. I. 2011, ApJ, 738, 72
  * Visscher et al. (2010b) Visscher, C., Moses, J. I., & Saslow, S. A. 2010b, Icarus, 209, 602
  * Vos et al. (2022) Vos, J. M., Faherty, J. K., Gagné, J., et al. 2022, ApJ, 924, 68
  * Walt et al. (2011) Walt, S. v. d., Colbert, S. C., & Varoquaux, G. 2011, Computing in Science & Engineering, 13, 22
  * Wang et al. (2015) Wang, D., Gierasch, P. J., Lunine, J. I., & Mousis, O. 2015, Icarus, 250, 154
  * Wang et al. (2016) Wang, D., Lunine, J. I., & Mousis, O. 2016, Icarus, 276, 21
  * Wright et al. (2014) Wright, E. L., Mainzer, A., Kirkpatrick, J. D., et al. 2014, AJ, 148, 82
  * Zahnle & Marley (2014) Zahnle, K. J., & Marley, M. S. 2014, ApJ, 797, 41
  * Zhang & Showman (2018) Zhang, X., & Showman, A. P. 2018, ApJ, 866, 1