# Multiple Patchy Cloud Layers in the Planetary Mass Object SIMP0136+0933

Allison M. McCarthy Department of Astronomy & The Institute for Astrophysical
Research, Boston University, 725 Commonwealth Ave., Boston, MA 02215, USA
Philip S. Muirhead Department of Astronomy & The Institute for Astrophysical
Research, Boston University, 725 Commonwealth Ave., Boston, MA 02215, USA
Patrick Tamburo Center for Astrophysics $|$ Harvard & Smithsonian, 60 Garden
Street, Cambridge, MA 02138, USA  Johanna M. Vos School of Physics, Trinity
College Dublin, The University of Dublin, Dublin 2, Ireland  Department of
Astrophysics, American Museum of Natural History, New York, NY 10024, USA
Caroline V. Morley Department of Astronomy, University of Texas at Austin,
Austin, TX 78712, USA  Jacqueline Faherty Department of Astrophysics, American
Museum of Natural History, New York, NY 10024, USA  Daniella C. Bardalez
Gagliuffi Department of Physics & Astronomy, Amherst College, 25 East Drive,
Amherst, MA 01003, USA  Department of Astrophysics, American Museum of Natural
History, New York, NY 10024, USA  Eric Agol Department of Astronomy & Virtual
Planetary Laboratory, University of Washington, Box 351580, U.W., Seattle, WA
98195, USA  Christopher Theissen Department of Astronomy & Astrophysics,
University of California, San Diego, 9500 Gilman Drive, La Jolla, CA
92093-0424, USA

(Received October 17, 2023; Accepted February 20, 2024)

###### Abstract

Multi-wavelength photometry of brown dwarfs and planetary-mass objects
provides insight into their atmospheres and cloud layers. We present near-
simultaneous J- and $K_{s}$-band multi-wavelength observations of the highly
variable T2.5 planetary-mass object, SIMP J013656.5+093347. We reanalyze
observations acquired over a single night in 2015 using a recently developed
data reduction pipeline. For the first time, we detect a phase shift between
J- and $K_{s}$-band light curves, which we measure to be
$39.9^{\circ+3.6}_{-1.1}$. Previously, phase shifts between near-infrared and
mid-infrared observations of this object were detected and attributed to
probing different depths of the atmosphere, and thus different cloud layers.
Using the Sonora Bobcat models, we expand on this idea to show that at least
two different patchy cloud layers must be present to explain the measured
phase shift. Our results are generally consistent with recent atmospheric
retrievals of this object and other similar L/T transition objects.

Brown Dwarfs (185), T dwarfs (1679), Stellar Atmospheres (1584), Exoplanet
Atmospheres (487), Exoplanet Atmospheric Variability (2020), Exoplanet
Atmospheric Structure (2310)

††journal: ApJ††software: astropy (Astropy Collaboration et al., 2022), SciPy
(Virtanen et al., 2020), Celerite2 (Foreman-Mackey, 2018), emcee (Foreman-
Mackey et al., 2013)

##  1 Introduction

Photometric variability monitoring is a useful tool for understanding the
atmospheric structure of brown dwarfs, planetary-mass objects, and exoplanets.
The observed variability is attributed to inhomogeneous cloud cover (Radigan,
2014), thermochemical instabilities (Tremblin et al., 2016), temperature
fluctuations(Robinson & Marley, 2014), and/or auroral activity (Hallinan et
al., 2015). Cloud structures typically persist for longer than one rotation
(Metchev et al., 2015; Apai et al., 2017), and as a result, we can probe the
presence of these clouds by measuring the brightness of a planetary-mass
object as it rotates. Past variability surveys in the infrared revealed trends
with spectral type (Radigan, 2014; Metchev et al., 2015), gravity (Vos et al.,
2019, 2022), and inclination (Vos et al., 2017). The objects within the L/T
transition, spanning from late-L to early-T, are known for having the largest
variability amplitude (e.g. Radigan, 2014; Metchev et al., 2015). Recently,
Ashraf et al. (2022) found that variable objects can be detected by
peculiarities in time-integrated spectra.

While single-wavelength variability monitoring observations have revealed a
large population of variable brown dwarfs and planetary-mass objects, multi-
wavelength monitoring allows us to study time-dependent effects across the
vertical structure of layers within the atmosphere. Additionally, simultaneous
multiwavelength photometry probes different depths (and thus
pressure/temperature levels) in the atmosphere providing insight into the
different cloud layers. The vertical cloud structure can reveal itself in
phase-shift and color changes across different band-passes. For example,
Buenzli et al. (2012) found that for the T6.5 dwarf 2MASS J22282889-431026,
phase lag increased for wavelengths probing decreasing pressure levels. Using
ground-based simultaneous multi-wavelength photometric observations Biller et
al. (2013) found that there were pressure-dependent phase shifts for the T0.5
object Luhman 16B (Luhman, 2013; Gillon et al., 2013). Biller et al. (2018)
detected phase offsets ranging from 200∘ to 210∘ between synthesized near-IR
light curves and a Spitzer mid-IR light curve for the L7 planetary mass object
PSO J318.5-22. Manjavacas et al. (2019) found a tentative phase shift between
J- and H-band light curves for the planetary mass T8 dwarf Ross 458 C.

In this paper we investigate the T2.5 dwarf SIMP J013656.5+093347 (SIMP0136
hereafter). SIMP0136 was discovered by Artigau et al. (2006) and is believed
to be a member of the 200-Myr-old Carina near-moving group (Gagné et al.,
2017). Gagné et al. (2017) estimated that SIMP0136 has a mass of $12.7\pm
1.0M_{\rm Jup}$, an effective temperature of $1098\pm 6$ K, and a $\log
g=4.31\pm 0.03$. In 2009, Artigau et al. found a rotation period of $\sim$2.4
hours, a peak-to-peak J-band amplitude of $\sim$ 50 mmag, and reported light
curve evolution from night-to-night. Apai et al. (2013), Metchev et al.
(2013), and Croll et al. (2016) each found a rotation period consistent with
Artigau et al. (2009)’s original $\sim$ 2.4-hour rotation period. Kao et al.
(2016) reported highly circularly polarized pulsed radio emission in the 4-8
GHz band and a magnetic field $\gtrsim 2.5$ kG, consistent with the presence
of aurorae on SIMP0136. The authors adopted a period of $2.88^{+0.34}_{-0.27}$
hr for the quasi-quiescent emission at X- band. A detailed investigation into
the inclination angle impacting observed color and variability found SIMP0136
to be close to equator-on at $80\pm 12^{\circ}$(Vos et al., 2017).

Artigau et al. (2009) also completed consecutive multi-wavelength observations
of SIMP0136 which revealed correlated J- and $K_{s}$-band light curves with a
$\frac{\Delta Ks}{\Delta J}$ = 0.48 $\pm$ 0.06. The authors showed that a two-
temperature surface could not explain the$\frac{\Delta Ks}{\Delta J}$ value,
and as a result they interpreted the variability results as clouds. Apai et
al. (2013) used the Hubble Space Telescope (HST) WFC3 grism mode to measure
spectroscopic variability and found that there were minimal J-H color changes.
The authors also determined that the surface is covered by a spatially
heterogeneous mix of two distinct regions: low-brightness, low-temperature
thick clouds, and brighter, thin, and warm clouds. Lew et al. (2020) also
reported a lack of J-H color change. No specific values are cited in the text
of the paper, but based on the linear fit shown in Figure 5, the $\Delta(J-H)$
appears to be 0, with minimal error. Yang et al. (2016) reported a phase shift
of $33.4\pm 3.9^{\circ}$ between simultaneous Spitzer 3.6 $\mu m$ and HST
light curves. Additionally, the authors report a range of phase shifts between
Spitzer 3.6 $\mu m$ and 4.5 $\mu m$ bands.

Using the Brewster retrieval framework and archival 1–15 $\mu m$ spectra, Vos
et al. (2023) completed an atmospheric retrieval of SIMP0136. Building on the
work of Apai et al. (2013), they found that a patchy, high-altitude forsterite
(Mg2SiO4) cloud above a deeper, optically thick iron cloud best describes
their spectral data. The authors determine that in order to reproduce the Apai
et al. (2013) HST variability amplitude observed in SIMP0136, the forsterite
cloud coverage changes by a small amount as the object rotates. Based on an
amplitude variability of 3%, Vos et al. (2023) find that the forsterite cloud
coverage fluctuates between 69% and 72%.

In this paper, we report the first observed phase shift between J- and
$K_{s}$-bands for SIMP0136. We reanalyze near-simultaneous J- and $K_{s}$-band
observations, originally taken in November 2015 using the 1.8-m Perkins
Telescope Observatory near Flagstaff, Arizona, and presented in Croll et al.
(2016), with the newly developed PINES Analysis Toolkit (PAT, Tamburo, 2022).
We use the Sonora Bobcat (cloud-free) models (Marley et al., 2021b) to
calculate the pressure/temperature levels probed by J- and $K_{s}$-band. We
find that J- and $K_{s}$-bands probe different regions of the atmosphere and
thus would be affected differently by the presence of different cloud levels.
We find that to explain the observed phase shifts, we must invoke the presence
of at least two different patchy cloud layers. In Section 2, we discuss our
observations and data reduction technique. In Section 3, we describe our
analysis and present our results, followed by our prediction of cloud
distribution in Section 4. In Section 5, we use a three-patch atmospheric
model to explain the observed color changes of SIMP0136. In Section 6, we
present our conclusions and future steps.

##  2 Observations and Data Analysis

The data presented in this paper were originally taken on UT 20 November 2015
as part of a larger campaign reported in Croll et al. (2016), which describes
the observations and sky conditions. Briefly, the observations were taken
using the Mimir near-infrared instrument (Clemens et al., 2007) at the 1.8-m
Perkins Telescope Observatory, located on Anderson Mesa near Flagstaff,
Arizona. The data consist of J- and $K_{s}$-band exposures taken sequentially
with 12, 30-second J-band exposures, followed by 24, 15-second $K_{s}$-band
exposures. The overhead time, including time for the readout, was $\sim$ 3
seconds for each observation, regardless of the band. This pattern continued
with minimal interruption for 8.5 hours. The guiding was performed by an off-
axis guide camera. The observations were taken in photometric conditions.
Calibration images were also obtained during the November 2015 run. Flat
fields were obtained using a lamp-illuminated dome flat screen. One hundred
fifty lamp-on and lamp-off images were taken in both J- and $K_{s}$-bands.
Twenty 30-second darks and 20 15-second darks were also taken.

We independently reduced the raw data using the newly developed PINES Analysis
Toolkit (PAT, Tamburo, 2022). PAT was originally designed to reduce data from
the Perkins INfrared Exosatellite Survey, PINES (Tamburo et al., 2022). PAT
has proven to be applicable, more broadly, to any photometric observations
taken using the Mimir Instrument at the Perkins Telescope Observatory. PAT is
explained in detail in Section 3 of Tamburo et al. (2022), but the procedure
will be explained briefly here. The following procedure is completed on the J-
and $K_{s}$-band data independently.

PAT begins by combining the calibration images and measuring a bad pixel mask,
including the hot, dark, and variable pixels on the Mimir detector. The data
are reduced using standard calibration techniques. The user then selects an
image which the pipeline then analyzes to find the target and suitable
reference stars using source detection procedures from the photutils Python
package. The user confirms this information before proceeding. This step is
repeated for each band. PAT then measures the centroid location of the target
and reference stars in each exposure. Next, the pipeline performs simple
aperture photometry with both fixed and time-variable apertures. Throughout
this process, if a bad pixel (identified by the bad pixel mask) falls on
either the target or reference star, the pipeline uses a 2D-Gaussian to
interpolate over the bad pixel. Finally, PAT produces light curves by using
the background-subtracted flux from the target and reference stars. PAT also
creates an “artificial light curve”, or ALC, a unitless weighted sum of
normalized reference star fluxes, following the approach of Murray et al.
(2020), to remove trends due to changes in the transmission in the Earth’s
atmosphere and weight the individual reference stars by their signal-to-noise
ratio and intrinsic variability. The pipeline then selects the optimal
aperture (fixed or time-variable) by selecting the light curve that minimizes
the average scatter over the duration of individual blocks of data. In this
context, the term block refers to one $\sim$6 minute set of observations
before the filter is switched. In the case of our data, the optimal aperture
for J-band was found to be fixed at a 2$\farcs$9 radius (5 pixels), and a
fixed 1$\farcs$7 radius (3 pixels) for $K_{s}$-band. The high background
present in $K_{s}$-band results in increased photon noise for larger
apertures.

In order to obtain near-simultaneous J- and $K_{s}$-band observations, the
filter was changed every $\sim$6 minutes. We measured the mean flux and error
on the mean flux for each of these 6-minute blocks of data, resulting in 37
binned points in J-band and 36 binned points in $K_{s}$-band. Figure 1 shows
the resulting light curves. The black markers represent the binned data, and
the small, light markers — orange dots for J-band, and pink x’s for
$K_{s}$-band — represent the unbinned data.

Since PAT reduces the J- and $K_{s}$-band data independently,the reference
stars used in the reduction are not necessarily the same in both J- and
$K_{s}$-band. As such, the reference stars used in the reduction only need to
have similar J- or $K_{s}$-band magnitudes as SIMP0136, not necessarily the
same $J-K_{s}$ color. The ALC is best used to measure changes in magnitudes,
not the magnitudes themselves.

We also explored the impact of precipitable water vapor on the measured flux
variations of SIMP0136 and the reference stars in each band. Figure 8 in
Tamburo et al. (2022) shows the impact of precipitable water vapor for J-band.
It shows that for a change in 1mm of precipitable water vapor, there is a
0.1%, or 1mmag, of variation, which is well within our measured variability of
SIMP0136. The $K_{s}$-band analysis was done in the same way, and found a
0.4%, or 4mmag, of variation, which is also within our measured variability of
SIMP0136. Additionally, we determined that there was no trend between
temperature of the reference star and increased variability in the object. We
concluded that the impact of precipitable water vapor is negligible.

##  3 Comparing the J- and $K_{s}$-band Light Curves

###  3.1 MCMC Modeling

Figure 1: J- and $K_{s}$-band light curves of SIMP0136. The J-band light curve
is shown above in orange with round markers. The $K_{s}$-band light curve is
shown below in purple with x markers. The small, light markers represent the
unbinned data. The large, black markers are the binned data with 1 $\sigma$
error bars, calculated using the standard deviation of the mean within a
block. The orange and purple lines show 1000 randomly selected MCMC samples
for the model light curves of J- and $K_{s}$-band, respectively. In both
plots, the black line is the maximum likelihood solution of the MCMC fits.

In order to determine the phase shift between the J- and $K_{s}$-band light
curves, we chose to cross-correlate the light curves. The observing strategy
results in the data being segmented into blocks and gaps. This segmentation
poses a challenge when calculating the phase shift between the J- and
$K_{s}$-band light curves, as it is not possible to cross-correlate data with
different x-axis timestamps. To address this issue, we chose to obtain a best-
fit model for the data and to use the model to interpolate the measurements
onto the same time stamps.

Using the reduced data from PAT, we used celerite 2 (Foreman-Mackey et al.,
2017; Foreman-Mackey, 2018), an algorithm for fast and scalable Gaussian
Processes (GP), to determine a best-fit model. The RotationTerm kernel,
designed for stellar rotation, uses a combination of two stochastically driven
damped simple harmonic oscillator (SHO) terms. The period of the second mode
is half that of the first mode’s period. The model had five free parameters
—sigma, period, Q0, dQ, and f —that are described fully in the celerite 2
documentation.

To determine best-fitting parameters and uncertainties in the model, we
performed a Markov chain Monte Carlo (MCMC) using the emcee package (Foreman-
Mackey et al., 2013). We chose a uniform prior for the period. The limits set
were between $\sim$0.9 and $\sim$3.6 hours based on previous measurements of
the system.

We ran 64 chains with 6,000 steps, with the first 1,000 discarded as burn-in.
The chains were well mixed by visual inspection. Figure 1 shows the results of
the MCMC as orange lines for J-band and purple lines for $K_{s}$-band, while
the black line shows the maximum likelihood solution of the MCMCs in both
plots. We also inspected the median posterior solution, which was nearly
identical.

###  3.2 J/$K_{s}$ Phase Shift

Figure 2: Maximum likelihood solution of the MCMC fits for the J- and
$K_{s}$-band data. The J-band data is shown in orange, with round markers, and
the $K_{s}$-band data is shown in purple, with x markers. The small, light
markers are the unbinned data.

Figure 2 shows the J and $K_{s}$ data and best-fit models on the same plot,
revealing similar features that are slightly out of phase. To quantify the
phase shift between the two light curves, we used the signal function in the
scipy Python package (Virtanen et al., 2020) to perform a cross-correlation
analysis on the maximum likelihood solution for each light curve.

To estimate the error on the cross-correlation, we randomly selected a chain
from each of the J- and $K_{s}$-band MCMC fits, and performed the cross-
correlation on those two randomly selected chains. We repeated this process
1000 times. The distribution of the phase shifts was right-skewed as shown in
Figure 3(b). Because of this, we report the 16th and 84th percentiles for the
lower and upper errors respectively. The result of this analysis is presented
in Figure 3(a), which shows a $39.9^{\circ+3.6}_{-1.1}$ phase shift between J
and $K_{s}$-band for SIMP 0136.

###  3.3 J and $K_{s}$ Period and Amplitude

From the results of the MCMC, we find a period of P=$2.55^{+0.10}_{-0.08}$
hours for J-band. For $K_{s}$-band, we find a period of P= $2.68\pm 0.28$
hours. The two period values are consistent to within 1$\sigma$.

Following the approach of Apai et al. (2013), we also measured the time
between peaks and troughs for both light curves. We used the same method to
determine the uncertainties on time between peaks as the phase shift. In both
light curves, but especially prominent in J-band, is a double peak feature.
For this study, we measured the time between the first peaks in each of the
double peak features.

For J-band, between the first and second troughs, $t_{t12}=2.55\pm 0.02$ hr,
and between the second and third troughs, $t_{t23}=2.52\pm 0.01$ hr. Between
the first and second peaks, $t_{p12}=2.49\pm 0.02$ hr, and between the second
and third peaks, $t_{p23}=2.40\pm 0.03$ hr.

For $K_{s}$-band, between the first and second troughs, $t_{t12}=2.70\pm 0.05$
hr, and between the second and third troughs, $t_{t23}=2.60\pm 0.05$ hr. For
the peaks in $K_{s}$-band, we only calculated one difference for the peaks
between $\sim 5.19$ and $\sim 2.41$ hr, $t_{p12}=2.76\pm 0.07$ hr.

In both J- and $K_{s}$-band, the amplitude of the light curve decreases over
the course of the observations. For J-band, we begin by measuring the
amplitude from the trough at $\sim 1.47$ hr to the peak at $\sim 2.37$ hr, and
find the amplitude is $3.13\pm 0.16\%$. In the next rotation, the J-band
amplitude decreases to $2.27\pm 0.10\%$. Finally, in the third rotation it
decreases again to $1.93\pm 0.12\%$. For $K_{s}$-band, beginning at the peak
at $\sim 0.75$ hr and measuring to the trough at $\sim 1.53$ hr, the amplitude
is $1.94\pm 0.27\%$. In the second rotation, the amplitude decreases to
$1.57\pm 0.19\%$. In the final rotation, the amplitude decreases again to
$1.17\pm 0.16\%$.

The change in amplitude throughout the night, as well as difference in timing
between peaks in the light curve, combined with the phase-shift between the
bands leads to color modulation throughout the course of the observation and
can be indicative of a rapidly evolving atmosphere(e.g. Metchev et al., 2015;
Apai et al., 2017).

(a)

(b)

Figure 3: Panel (a) shows the cross correlation results of the J- and
$K_{s}$-band light curve MCMC model predictions. The phase shift is found to
be 0.3 hr = $39.9^{\circ+3.6}_{-1.1}$. Panel (b) shows the distribution of
phase shifts using the cross-correlation technique between 1000 random model
fits for each J- and $K_{s}$-band. We determined the errors on the phase shift
by calculating the 16th and 84th percentile of the data, shown in dashed
lines. Since the distribution of phase shifts is right-skewed, this led to a
larger positive error than a negative error: $39.9^{\circ+3.6}_{-1.1}$.

##  4 Cloud Distribution Scenarios in the Atmosphere of SIMP0136

###  4.1 Sonora Bobcat Models

Simultaneous multiwavelength photometry probes different depths, and thus
pressure/temperature levels in the atmosphere of brown dwarfs and planetary
mass objects. These observations provide insight into the vertical cloud
structure in the atmosphere of these objects. To calculate the flux
contribution from different pressure levels, we used the cloudless Sonora
Bobcat models (Marley et al., 2021b, a; Lupu et al., 2021).

The Sonora Bobcat models are available for a range of atmospheres with varying
effective temperatures, gravities, and metallicities. Using an atmospheric
retrieval analysis of SIMP0136, Vos et al. (2023) found a [M/H] of $-0.32\pm
0.05$ and a C/O of $0.79\pm 0.02$. Additionally, they found a $T_{\rm
eff}=1329\pm 17$ K, which they conclude is too high, and the SED-calculated
temperature would be a better estimate. The authors’ SED temperature estimate
is $T_{\rm eff}=1150\pm 70$ K, which is in agreement with a previous result
from Gagné et al. (2017) of $1098\pm 6$ K. The estimate of the temperature
found from the retrieval is determined in conjunction with the radius. The
retrieval results also find a small radius: $0.81\pm$ 0.03 $R_{Jup}$. The
small radius, and associated large temperature, are not consistent with
evolutionary models, but are a common feature of retrieval-derived radii (e.g.
Burningham et al., 2021; Mollière et al., 2020; Barman et al., 2011). Finally,
Vos et al. (2023) find the log g = $4.5\pm 0.4$ calculated from retrievals to
be similar to the log g = $4.25\pm 0.08$ calculated from the SED fit. Using
these values as a reference, and working within the available values from the
Sonora Bobcat models, we used the atmospheric structure model with the
following parameters: [M/H]=-0.5, solar C/O ratio, $T_{eff}$=1100 K, and g=
316 m/s (log g= 4.5). We used the correlated $k$-coefficients (described fully
in Lacis & Oinas, 1991) with the following parameters: [M/H]=-0.5, and solar
C/O ratio. The solar C/O ratio which Sonora Bobcat uses is from Lodders (2010)
value for the solar C/O=0.458.

To begin, we downloaded the appropriate correlated $k$-coefficients from Lupu
et al. (2021). We then obtained the pressure-temperature (PT) profiles from
the Sonora Bobcat atmospheric structure files (Marley et al., 2021a). We
interpolated the temperature values onto a new axis that was evenly spaced in
$\Delta$ log P to create a more finely sampled profile. We then interpolated
the k-coefficient grid onto the PT profile for the model atmosphere.

Making the assumption that the atmosphere is in local thermal equilibrium, we
solved the radiative transfer equation along the PT profile to determine the
flux contribution from each layer at each wavelength. The results are shown in
the contour plot of Figure 4. The dotted lines shown in Figure 4 contain the
pressure levels at which 80% of the flux originates.

The observed data does not include H-band measurements, but we chose to
include the H-band wavelengths in Figure 4 for the sake of mapping the
complete near infrared wavelengths. The plotted dots and associated error bars
show the pressure levels at which 80% of the total flux within each band, J,
H, and $K_{s}$, originates (similar to Figure 5 in Biller et al., 2013). To
determine these values, we summed the radiative transfer equation results for
each wavelength window in each bandpass to obtain the total emergent flux.
Next, we multiplied the emergent flux and the transmission curve for the band.
We then determined the PT levels corresponding to 10% and 90% of the total
flux in each bandpass and calculated a cumulative summation of each PT level
until the PT levels from below which the flux accounted for 10% and 90% of the
total flux in the band were obtained. The location of the marker corresponds
to the weighted average wavelength and the maximum flux contribution pressure
level for each band. The plotted data points for each band do include the
transmission curve in the calculations for the location of the marker and the
size of the error bars. The pressures for which the center 80% of the received
(emitted flux multiplied by the transmission curve) flux originates for each
band are 12.6 –26.2 bar for J-band, 3.6 –20.0 bar for H-band, and 1.9 –14.1
bar for $K_{s}$-band.

The lower plot in Figure 4 displays the NIR spectrum of the Sonora Bobcat
model in black and the transmission curves for the J-, H-, and $K_{s}$-,
bandpasses in purple. The flux units are in megawatts. Similarly, Figure 5
shows the total flux per pressure for J- and $K_{s}$-band with the units in
megawatts. This plot highlights the different regions of the atmosphere that
the J- and $K_{s}$-band probe. The gray bars represent the pressure levels at
which the clouds identified by Vos et al. (2023) are situated.

Figure 4: Top: The contour plot shows the distribution of flux as a function
of wavelength and pressure for a cloudless atmosphere with effective
temperature of 1100K, log(g)=4.5, [M/H]=-0.5, and solar C/O. The color bar
distribution is logarithmic. The dotted black lines represent the pressure
levels of the atmosphere where 10% and 90% of the flux is emitted. The circle
markers are located on the weighted average wavelength and peak flux emission
pressure level for each bandpass while the bars show the range of pressure
levels where 80% of the emitted flux is weighted by the band transmission.
Bottom: This plot displays the Near Infrared spectra of the Sonora Bobcat
model in black, and the transmission curves for J-, H-, and $K_{s}$-band in
purple. Figure 5: Received flux distribution for the J- and $K_{s}$-band pass
wavelength range as a function of pressure for a cloudless atmosphere. The
solid orange line and dotted purple line represent the J- and $K_{s}$-band
flux distribution, respectively. The gray boxes in the plot depict the regions
of the atmosphere which Vos et al. (2023) found SIMP0136 contains forsterite
(Mg2SiO4) and iron (Fe) clouds.

###  4.2 Multiple Atmospheric Layers Responsible for Phase Shift

While a combination of any three unique atmospheric regions or layers could
create J- and $K_{s}$-band variability, as well as the phase shift between the
two light curves, previous studies such as Vos et al. (2023) have found strong
evidence for the presence of clouds. Therefore, we considered different
scenarios of cloud distribution in the atmosphere of SIMP0136 to determine
which would result in the modulations in the J- and $K_{s}$-band, as well as
the phase shift between the two light curves. In these scenarios, we consider
clouds to be blocking a portion of the flux from the region of the atmosphere
which the band probes. We acknowledge that, in reality, changing cloud
opacities within or above the portion of the atmosphere that the band is
probing, will alter the distribution of the emergent flux per pressure level.

We considered four simple cloud distribution scenarios to determine what could
be causing the phase shift between J- and $K_{s}$-band for SIMP0136:

  1. 1.

No clouds within or above the regions of the atmosphere which J- and
$K_{s}$-band probe. In this scenario, there would be no clouds to block the
emergent flux from the regions of the atmosphere that either the J- or
$K_{s}$-band probe. If this were the case, we would see no modulation in
either light curve.

  2. 2.

Patchy upper atmosphere clouds either within or above the region of the
atmosphere which $K_{s}$-band probes. In this scenario, the clouds in the
upper atmosphere would block the emergent flux from both the J- and
$K_{s}$-band regions of the atmosphere. This would create modulation in both
the J- and $K_{s}$- light curves and the modulation would be in-phase.

  3. 3.

Patchy lower atmosphere clouds solely within the region of the atmosphere
which J-band probes. In this scenario, the emergent flux from the portion of
the atmosphere that $K_{s}$-band probes, would be unobstructed, and there
would be no modulation in the $K_{s}$-band light curve. However, the clouds
would cause disruptions to the emergent flux from the region of the atmosphere
which J-band probes and the J-band light curve would reflect this by being
variable.

  4. 4.

Upper atmosphere either above or within the region that $K_{s}$-band probes
and lower atmosphere clouds within the region of the atmosphere which J-band
probes. In this scenario, the presence of clouds in both the upper and lower
atmosphere would produce variable light curves for both J- and $K_{s}$-band.
However, since the lower atmosphere clouds would only affect the emergent flux
from the region of the atmosphere that the J-band probes, the light curves
would likely not be in phase. The light curves could be in phase if the upper
and lower atmosphere clouds were of the same size, shape, and in-phase.

Based on these four simple atmosphere structure scenarios for SIMP0136, we
conclude that there must be at least two patchy cloud layers, one above or
within the atmospheric region that $K_{s}$-band probes (hereafter referred to
as upper atmosphere clouds), and one between the regions that $K_{s}$- and
J-band probe or within the J-band region (hereafter referred to as lower
atmosphere clouds), in order to produce the phase-shift seen in our data.

These conclusions build on the work of Vos et al. (2023) (which in turn builds
on the work of Apai et al., 2013) who find that the spectrum of SIMP0136 is
best described by a model containing patchy forsterite clouds lying above an
iron cloud deck. Vos et al. (2023) show that the forsterite clouds exist at a
pressure level of 1.3-1.7 bar, and the iron clouds exist at a pressure level
of and below $\sim$7 bar. We used these pressures for the basis of the
location of the gray bars shown in Figure 5. The atmospheric structure
suggested by Vos et al. (2023) matches scenario (2)—a patchy upper atmosphere
cloud over a “uniform” lower cloud deck. We argue that in order to explain the
phase shift, both cloud layers must be patchy, matching atmospheric structure
scenario (4) discussed in Section 4.2.

As shown in Figure 5 of Vos et al. (2023), the addition of an iron cloud layer
into the atmospheric model does alter the flux contribution per pressure
level, especially around 1-1.3 microns. The addition of the forsterite patch
raises the pressures probed in J-band to $\sim$0.6 bar. This discrepancy
between the clear atmospheric model, and the Brewster retrieval model from Vos
et al. (2023), while important to consider, does not alter the fact that in
order to explain the phase shift, multiple patchy cloud layers must be
present.

In order to create a phase-shift between any two curves, the addition of
another curve must be present. For the atmosphere of SIMP0136, the varying
cloud coverage in each cloud layer influences the observed light curve at each
wavelength. When forsterite clouds are present, it follows that the optical
depth of both J- and $K_{s}$-band will rise to a lower pressure level in the
atmosphere. However, in the region with patches of the forsterite cloud, the
J-band observations are still able to probe down to the level which the iron
clouds exist, while even in a hole in the forsterite cloud layer, the
$K_{s}$-band observations cannot probe down to the depth of the atmosphere
where iron-clouds exist. This is highlighted by our clear atmosphere model
shown in Figure 4. So, despite which pressure layer we are able to probe when
forsterite clouds are present, the important factor is that when forsterite
clouds are not present, the J-band can probe down to the iron cloud layer, and
observe its modulations, while the $K_{s}$-band cannot.

##  5 Implementation of Three Patch Atmospheric Model

In an attempt to model the observed phase shift, we used a combination of the
Sonora Bobcat (clear atmosphere) models (Marley et al., 2021b), the Sonora
Diamondback (cloudy atmosphere) models (Morley et al., 2024), and bespoke
Sonora Diamondback iron-only cloudy models to create a three-patch model
atmosphere that could describe the observed period, amplitude, and phase shift
when each patch was modulated. The standard Sonora Diamondback models included
five values for the sedimentation efficiency parameter —$f_{sed}$=1,
$f_{sed}$=2, $f_{sed}$=3, $f_{sed}$=4, and $f_{sed}$=8, (lower
$f_{sed}$=cloudier) —and the bespoke Sonora Diamondback iron-only models
included four $f_{sed}$ values —$f_{sed}$=2, $f_{sed}$=3, $f_{sed}$=4, and
$f_{sed}$=8. In each of the following subsections, we describe the outcomes of
each different patchy atmosphere which we created. Subsection 5.1 explores a
three-patch model with a cloudy patch (including forsterite clouds), an iron
only patch, and a clear patch all at an effective temperature of 1100K.
Subsection 5.2.1 considers an entirely clear atmosphere with three patches at
three different effective temperatures, while subsection 5.2.2 considers an
entirely cloudy atmosphere with each of the three patches having three
different effective temperatures. Subsections 5.3.1, 5.3.2, and 5.3.3, combine
clear and cloudy atmospheric patches at different effective temperatures. In
each case, our procedure was the same.

We combined three spectra, from the different models, representing the
different patches of the atmospheres, at different ratios. Representing the
three different patches spectra as $s_{1}$, $s_{2}$, and $s_{3}$, their three
ratios as $a$, $b$, and $c$ (such that $a+b+c=1$), and the total combined
spectra as $S$, we combine the three spectra as

$S=(a*s_{1})+(b*s_{2})+(c*s_{3}).$

We use the 2MASS transmission curves to determine the cumulative flux in each
of the J-, H- and $K_{s}$-bands, and use the 2MASS zero point flux
measurements (Cohen et al., 2003) to convert the measured flux in each band
from our three patch atmosphere to a magnitude. We then compared each of the
different combined three patch atmosphere colors and magnitudes to our
observed data of SIMP0136. In order to convert our SIMP0136 observations to
magnitudes, we assume that the normalized flux is equal to the currently
accepted J-, H-, and $K_{s}$-band magnitudes (Zhang et al., 2021) which is
shown in Figure 6. The J-$K_{s}$ color difference of the MCMC solutions is
shown in Figure 7.

Ideally, simultaneous observations at a wide range of wavelengths would
provide better constraints when modeling SIMP 0136’s atmosphere. Previous
studies (e.g. Lew et al., 2020) have found minimal J-H color change,
$\Delta(J-H)\approx 0$, throughout the rotation of SIMP0136. The lack of color
change could be consistent with the conclusions drawn from our contribution
plot (Figure 4), and the results of Vos et al. (2023), which place a
forsterite cloud between 1.3 and 1.7 bar, above the regions that both J- and
H-band probe, but iron clouds within the region that both J- and H-band probe.
Therefore, any patchiness in either cloud layer would affect both J- and
H-band in the same capacity. Since we do not have H-band measurements that are
simultaneous or near-simultaneous with our J- and $K_{s}$-band measurements,
we use the results of Lew et al. (2020) and inspect both the J-$K_{s}$ and J-H
color evolution to determine a model’s accuracy. To do this, we first match
our observed J-$K_{s}$ color evolution to a successful model and record the
patch ratios that align with the observed data. We do this by taking 30-second
steps along the J-band and $J-K_{s}$ evolution, and determine if there is a
point in our model within 0.001 magnitudes from both the J-band and $J-K_{s}$
magnitudes. From there, any patch coverage combination that was successful
along the J versus J-K evolution is transferred to the J-H color space. We are
then able to observe the predicted J-H behavior. An example is shown in Figure
8, where the purple line in panel (a) shows the observed data, and the orange
stars in panel (a) and panel (b) show the points of successful patch fraction.

Figure 6: Maximum likelihood solution of the MCMC fits for the J- and
$K_{s}$-band data. The J-band data is shown in orange, with round markers, and
the $K_{s}$-band data is shown in purple, with x markers. The small, light
markers are the unbinned data. The J-band magnitudes are shown on the left
axis, and the $K_{s}$-band magnitudes are shown on the right axis. Figure 7:
The J-$K_{s}$ color difference of the MCMC solutions. The purple line is
colored by time since initial observation, where the darkening of the line
signals increasing time since initial observation. The gray lines shows the
J-$K_{s}$ color difference of 1000 random MCMC samples. The pink squares on
the color magnitude diagram and pink lines on the color bar show the time of
one rotation based on the J-band period of 2.55 hr.

###  5.1 A Cloudy Patch, Iron-Only Patch, and Clear Patch at 1100K

The final scenario discussed in Section 4.2, in which there are three viewing
possibilities —observing a forsterite cloud, observing through a patch in the
forsterite cloud layer down to an iron cloud, or observing a patch in both the
forsterite and iron cloud layers (clear atmosphere)—is described in this
section. Therefore, we created atmospheres with three patches, all at
effective temperatures of 1100K. In each atmosphere, the first patch used one
$f_{sed}$ value from the Sonora Diamondback (cloudy) model which includes both
forsterite and iron cloud layers. The second patch used one $f_{sed}$ value
from the bespoke Sonora Diamondback (iron-only) model, and the third patch
used the Sonora Bobcat (clear) atmosphere.

In this scenario, multiple model combinations, meaning different combinations
of $f_{sed}$ values for patch 1 and patch 2, were able to explain our observed
J-$K_{s}$ color change. However, when these same models were translated into
J-H color space, all of the models predicted a significant J-H color change
throughout the rotation period, conflicting with previous observations. In one
such model, shown in Figure 8, the $\Delta(J-H)=0.12$.

The model shown in Figure 8 uses an $f_{sed}=3$ for both the Sonora
Diamondback cloudy atmosphere models, and the bespoke Sonora Diamondback iron-
only models. Combining our observations of SIMP0136 with these models, shows
that the total fractional cloud coverage would need to vary between 31-81%,
with the forsterite cloud patch varying between 15-30%, and the iron cloud
patch varying between 1-66%. The forsterite and iron cloud patches are out of
phase. This is expected since a patch without forsterite clouds would allow
for more iron clouds to be visible, and contribute more to the overall flux.

While this model can explain 93% of the observed data, the model cannot
explain the entirety of the observation (namely the section from 2.2-2.6 hr),
and we attribute this to a combination of several factors. First, there may be
inaccuracies in the models. Second, an error in the radius measurements of
SIMP0136, which are needed to convert the model flux values to magnitudes, and
third, an error in the magnitudes arising from using the measured 2Mass values
as a reference point. At the peak variability of the J-band light curve, from
the trough at $\sim$ 1.47 hr to the peak at $\sim$ 2.37 hr, the amplitude is
$3.13\pm 0.16$%. The cloud coverage at the trough is $\sim$ 15%
forsterite+iron and $\sim$ 68% iron only, corresponding to $\sim$83% complete
cloud cover. At the time of the peak, our data is outside of the model. Our
best approximation would be $\sim$ 31% forsterite+iron clouds, and $\sim$ 0%
iron only cloud coverage. This corresponds to an increase in the
forsterite+iron cloud coverage of $\sim$ 16% decrease in iron cloud coverage
of $\sim$ 68%, and a decrease in the total cloud coverage of $\sim$ 52%.

As mentioned, other model combinations were successful at describing our
observed J-$K_{s}$ color change. Interestingly, there is little to no change
when changing either patch 1 or patch 2’s $f_{sed}$ from $f_{sed}=3$ to
$f_{sed}=4$. However, when changing either patch to $f_{sed}=1$ or
$f_{sed}=2$, while the other patch remains a higher $f_{sed}$, the models can
no longer match the data. Once both patches have an $f_{sed}=1$ or
$f_{sed}=2$, the model can match the data again. When the patches have a lower
$f_{sed}$, the overall cloud coverage drops to a maximum of 35%. The
forsterite and iron clouds are still out of phase, and both species coverage
fraction fluctuates $\sim$30%.

(a)

(b)

(c)

Figure 8: A three-patch model where patch 1 uses the Sonora Diamondback models
including a layer of forsterite clouds, patch 2 uses the bespoke Sonora
Diamondback iron-only models, and patch 3 uses the Sonora Bobcat clear
atmospheric models. All three patches have a temperature of 1100K. In this
model, both cloudy atmospheres had an $f_{sed}=3$. Panel (a) shows $M_{J}$ and
$J-K_{s}$ values from the highest posterior results of the MCMC of our data
(Figure 1) in purple. The orange stars represent the points in our model that
most closely match our observed data. The fractional coverage ratios that the
orange stars represent are then translated into the J-H color space, shown in
panel (b). In panel (b) the orange line traces the path that the three patch
model predicts. In both panel (a) and panel (b),the gray line represents the
boundary at which there is 100% cloud cover (any combination of forsterite and
iron), and 0% contribution by the Sonora Bobcat clear model. The pink line
represents the boundary where there is 0% contribution by the iron only model,
so the entirety of the flux is created by the Sonora Bobcat clear model and
Sonora Diamondback model which includes forsterite. The corresponding line for
when there is 0% contribution from the Sonora Diamondback which contains
forsterite is outside the plotted color-magnitude space, and is thus not
included. Panel (c) shows the variation of each patch coverage over time. The
purple stars represent patch 1, the pink stars represent patch 2, and the gray
stars show (patch 1 + patch 2) or (1 - patch 3). Panel (c) uses a polynomial
fit to trace the change in patch coverage over time. The line is only there as
a data visualization tool. From 2.0-2.6 hours, the model cannot reproduce the
measured magnitudes, and in panel (c) this is highlighted by the polynomial
fit. There is no physical meaning to the fit within this time range.

###  5.2 Inhomogenous Surface Temperature Patches for Both Entirely Clear and
Entirely Cloudy Atmospheres

For the sake of completeness, in this section we consider whether a three-
temperature model can also explain the observed phase shift in J- and
$K_{s}$-band, and lack of phase shift between J\- and H-band.

####  5.2.1 Clear Atmosphere

Robinson & Marley (2014) found that atmospheric temperature fluctuations could
be responsible for observed photometric variability in brown dwarfs.

Using the Sonora Bobcat (clear) atmospheric models, we found that in order to
create a red enough atmosphere to match our observed range in J-$K_{s}$, we
required three patches at effective temperatures of 1500K, 900K, and 800K. In
this scenario, our $\Delta(J-H)=0.03$, which is still different than the
previous J-H observations of SIMP0136.

####  5.2.2 Cloudy Atmosphere

Using different effective temperature patches of the standard Sonora
Diamondback (cloudy) models, we are unable to replicate the atmosphere of
SIMP0136. The spectra are too red to match the observed range in J-$K_{s}$. It
should be noted that as described in Section 5, our J- and $K_{s}$-band
magnitudes were based on the currently accepted values for SIMP0136. It is
possible that at the time our observations were completed, the J- and
$K_{s}$-band magnitudes were different than at the time the currently accepted
values were measured.

###  5.3 A Combination of Clear and Cloudy Patches at Different Temperatures

Finally, we used patches of different species of clouds, at different
$f_{sed}$ values, at different effective temperatures, to attempt to describe
the observed J-$K_{s}$ color change and lack of J-H color change.

####  5.3.1 Two Clear Patches at Different Temperatures, with a Third Patch
which is Cloudy

Multiple models of combining two clear patches at different effective
temperatures with a third cloudy patch were able to explain the observed
J-$K_{s}$ color change. However, none of them are able to explain the
previously observed J-H color change. For one such case, where patch 1 was a
clear patch at $T_{eff}=$1300K, patch 2 was a clear patch at $T_{eff}=$1200K,
and patch 3 was a cloudy patch at $T_{eff}=$900K and $f_{sed}=4$, the
$\Delta(J-H)=0.02$.

####  5.3.2 Two Cloudy Patches at Different Temperatures, with a Third Patch
which is Clear

Many different combinations of two cloudy patches with a third clear patch at
different effective temperatures were able to explain the observed J-$K_{s}$
color change but not the previously observed J-H color change. For a clear
patch at $T_{eff}=$1200K, and cloudy patches at $T_{eff}=$900K and
$T_{eff}=$1100K, both with $f_{sed}=3$, the resulting $\Delta(J-H)=0.02$.

####  5.3.3 A Cloudy Patch, Iron-Only Patch, and Clear Patch, all at Different
Temperatures

Several different combinations of effective temperatures and $f_{sed}$ values
were able to explain the observed J-$K_{s}$ color change, but unable to
explain the lack of J-H color change. For a model atmosphere with a clear
patch at $T_{eff}=$1000K, an iron-only patch with an effective temperature of
$T_{eff}=$1100K and $f_{sed}=3$, and a cloudy patch with an effective
temperature of $T_{eff}=$1200K and $f_{sed}=4$, the $\Delta(J-H)=0.01$. This
scenario represents our most successful attempt at matching both the observed
J-$K_{s}$ color change, as well as the previously observed lack of J-H color
change. As mentioned in the introduction, Lew et al. (2020) does not mention
specific values for the $\Delta(J-H)$ color change, but based on the
appearance of the linear fit in Figure 5, we can only approximate it to
$\Delta(J-H)=0$ with zero uncertainty. Without having simultaneous J-H
measurements of our own, we are unable to say with certainty that this model
accurately represents the atmosphere of SIMP0136.

##  6 Conclusion

In this work, we present the first detected J- and $K_{s}$-band phase shift
measurement for SIMP0136, a planetary-mass T2.5 dwarf. We determine the phase
shift to be $39.9^{\circ+3.6}_{-1.1}$. The period of the light curve is found
to be $2.55^{+0.10}_{-0.08}$ and $2.68\pm 0.28$ hours, for J- and $K_{s}$-band
respectively. The J- and $K_{s}$-band periods are consistent to within
1$\sigma$. Both the J- and $K_{s}$-band light curves have amplitudes which
decrease over time. The J-band amplitude decreases from $3.13\pm 0.16\%$ to
$2.27\pm 0.10\%$ to $1.93\pm 0.12\%$. For $K_{s}$-band, the amplitude
decreases from $1.94\pm 0.27\%$ to $1.57\pm 0.19\%$ to $1.17\pm 0.16\%$.

Additionally, using the Sonora Bobcat models (Marley et al., 2021b), we
calculated the flux distribution of a cloudless model atmosphere for SIMP0136
as a function of pressure and wavelength (Figures 4, 5), and find that the J-,
H-, and $K_{s}$-band primarily probe the pressures 12.6–26.2, 3.6–20.0, and
1.9–14.1 bar, respectively.

We find that among the explored various non-homogeneous temperature and cloudy
scenarios tested in this paper multiple patchy cloud layers is the preferred
model in order to describe the phase shift detected between the J- and
$K_{s}$-band light curves. The two patchy layers include a patchy upper layer,
which sits within or above the region of the atmosphere that the $K_{s}$-band
probes, and a patchy lower cloud layer that sits either between the regions of
the atmosphere which J- and $K_{s}$-band probe or within the region of the
atmosphere that the J-band probes (final scenario in Section 4.2). This
expands on the results of Vos et al. (2023), who find that the spectrum of
SIMP0136 is best described by a model containing patchy forsterite clouds at a
pressure range of 1.3–1.7 bar above an iron cloud deck sitting at and below
$\sim$7 bar.

We attempted to model the atmosphere of SIMP0136 by combining multiple spectra
as different patches representing fractional coverage in the atmosphere. We
used the Sonora Bobcat (clear) atmospheric models, the Sonora Diamondback
(cloudy) atmospheric models containing forsterite and iron clouds, and bespoke
Sonora Diamondback models which contain only iron clouds. In every combination
of model we could create, we were unable to match both the observed
J-$K_{s}$=0.03 color change and previously measured lack of J-H ($\Delta
J-H\approx 0$) color change from Lew et al. (2020).

The combination of different spectra is an approximation of the potential
atmospheric structure of SIMP0136. In every combination, the models are not
able to replicate the previously observed lack of J-H color change from Lew et
al. (2020). However, these observations are not simultaneous or near-
simultaneous with the J- and $K_{s}$-band observations presented in this work,
and it is unclear how much of a problem this poses. If the minimal J-H color
change is a constant feature, this suggests that the atmosphere of these
objects is much more complex, and our understanding of the cloud layers
between 1 and 20 bar is incomplete. Moving forward, simultaneous or near-
simultaneous observations of the J-, H-, and $K_{s}$-band will be necessary to
gather information about the atmosphere of these objects, and the atmospheric
layers that the near infrared bands probe.

Forthcoming JWST observations will reveal far more about the atmospheres of
substellar objects given that multiple programs will be monitoring objects
over extended wavelength coverage.Long term monitoring is necessary to
determine if phase shifts, amplitudes, and color-changes are constant
overtime. Simultaneous multiwavelength photometry combined with the power of
atmospheric models, such as the Sonora Bobcat and Diamondback models, but also
3D atmospheric modeling, provide an opportunity to create a more complete map
of the atmospheric structure of exoplanet analogs, such as SIMP0136.

We thank the reviewer for their helpful comments, which improved the quality
of this publication. We also thank Michael Line for his insight into the use
of atmospheric models, such as Sonora Bobcat. This material is based upon work
supported by the National Aeronautics and Space Administration under grant No.
80NSSC20K0256 issued through the Science Mission Directorate. This material is
based upon work supported by the National Science Foundation Graduate Research
Fellowship under Grant No. DGE-1840990. J. M. V. acknowledges support from a
Royal Society - Science Foundation Ireland University Research Fellowship
(URF$\backslash$1$\backslash$221932).

## References

  * Apai et al. (2013) Apai, D., Radigan, J., Buenzli, E., et al. 2013, ApJ, 768, 121, doi: 10.1088/0004-637X/768/2/121
  * Apai et al. (2017) Apai, D., Karalidi, T., Marley, M. S., et al. 2017, Science, 357, 683, doi: 10.1126/science.aam9848
  * Artigau et al. (2009) Artigau, É., Bouchard, S., Doyon, R., & Lafrenière, D. 2009, ApJ, 701, 1534, doi: 10.1088/0004-637X/701/2/1534
  * Artigau et al. (2006) Artigau, É., Doyon, R., Lafrenière, D., et al. 2006, ApJ, 651, L57, doi: 10.1086/509146
  * Ashraf et al. (2022) Ashraf, A., Bardalez Gagliuffi, D. C., Manjavacas, E., et al. 2022, ApJ, 934, 178, doi: 10.3847/1538-4357/ac7aab
  * Astropy Collaboration et al. (2022) Astropy Collaboration, Price-Whelan, A. M., Lim, P. L., et al. 2022, apj, 935, 167, doi: 10.3847/1538-4357/ac7c74
  * Barman et al. (2011) Barman, T. S., Macintosh, B., Konopacky, Q. M., & Marois, C. 2011, ApJ, 733, 65, doi: 10.1088/0004-637X/733/1/65
  * Biller et al. (2013) Biller, B. A., Crossfield, I. J. M., Mancini, L., et al. 2013, ApJ, 778, L10, doi: 10.1088/2041-8205/778/1/L10
  * Biller et al. (2018) Biller, B. A., Vos, J., Buenzli, E., et al. 2018, AJ, 155, 95, doi: 10.3847/1538-3881/aaa5a6
  * Buenzli et al. (2012) Buenzli, E., Apai, D., Morley, C. V., et al. 2012, ApJ, 760, L31, doi: 10.1088/2041-8205/760/2/L31
  * Burningham et al. (2021) Burningham, B., Faherty, J. K., Gonzales, E. C., et al. 2021, MNRAS, 506, 1944, doi: 10.1093/mnras/stab1361
  * Clemens et al. (2007) Clemens, D. P., Sarcia, D., Grabau, A., et al. 2007, PASP, 119, 1385, doi: 10.1086/524775
  * Cohen et al. (2003) Cohen, M., Wheaton, W. A., & Megeath, S. T. 2003, AJ, 126, 1090, doi: 10.1086/376474
  * Croll et al. (2016) Croll, B., Muirhead, P. S., Lichtman, J., et al. 2016, arXiv e-prints, arXiv:1609.03587.  https://arxiv.org/abs/1609.03587
  * Foreman-Mackey (2018) Foreman-Mackey, D. 2018, Research Notes of the American Astronomical Society, 2, 31, doi: 10.3847/2515-5172/aaaf6c
  * Foreman-Mackey et al. (2017) Foreman-Mackey, D., Agol, E., Ambikasaran, S., & Angus, R. 2017, AJ, 154, 220, doi: 10.3847/1538-3881/aa9332
  * Foreman-Mackey et al. (2013) Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306, doi: 10.1086/670067
  * Gagné et al. (2017) Gagné, J., Faherty, J. K., Burgasser, A. J., et al. 2017, ApJ, 841, L1, doi: 10.3847/2041-8213/aa70e2
  * Gillon et al. (2013) Gillon, M., Triaud, A. H. M. J., Jehin, E., et al. 2013, A&A, 555, L5, doi: 10.1051/0004-6361/201321620
  * Hallinan et al. (2015) Hallinan, G., Littlefair, S. P., Cotter, G., et al. 2015, Nature, 523, 568, doi: 10.1038/nature14619
  * Kao et al. (2016) Kao, M. M., Hallinan, G., Pineda, J. S., et al. 2016, ApJ, 818, 24, doi: 10.3847/0004-637X/818/1/24
  * Lacis & Oinas (1991) Lacis, A. A., & Oinas, V. 1991, J. Geophys. Res., 96, 9027, doi: 10.1029/90JD01945
  * Lew et al. (2020) Lew, B. W. P., Apai, D., Zhou, Y., et al. 2020, AJ, 159, 125, doi: 10.3847/1538-3881/ab5f59
  * Lodders (2010) Lodders, K. 2010, in Astrophysics and Space Science Proceedings, Vol. 16, Principles and Perspectives in Cosmochemistry, 379, doi: 10.1007/978-3-642-10352-0_8
  * Luhman (2013) Luhman, K. L. 2013, ApJ, 767, L1, doi: 10.1088/2041-8205/767/1/L1
  * Lupu et al. (2021) Lupu, R., Freedman, R., & Visscher, C. 2021, Correlated k coefficients for H2-He atmospheres; 180 spectral windows and 1060 pressure-temperature points, v1.0.0, Zenodo, doi: 10.5281/zenodo.4755012
  * Manjavacas et al. (2019) Manjavacas, E., Apai, D., Lew, B. W. P., et al. 2019, ApJ, 875, L15, doi: 10.3847/2041-8213/ab13b9
  * Marley et al. (2021a) Marley, M., Saumon, D., Morley, C., et al. 2021a, Sonora Bobcat: cloud-free, substellar atmosphere models, spectra, photometry, evolution, and chemistry, v1.0.0, Zenodo, doi: 10.5281/zenodo.5063476
  * Marley et al. (2021b) Marley, M. S., Saumon, D., Visscher, C., et al. 2021b, ApJ, 920, 85, doi: 10.3847/1538-4357/ac141d
  * Metchev et al. (2013) Metchev, S., Apai, D., Radigan, J., et al. 2013, Astronomische Nachrichten, 334, 40, doi: 10.1002/asna.201211776
  * Metchev et al. (2015) Metchev, S. A., Heinze, A., Apai, D., et al. 2015, ApJ, 799, 154, doi: 10.1088/0004-637X/799/2/154
  * Mollière et al. (2020) Mollière, P., Stolker, T., Lacour, S., et al. 2020, A&A, 640, A131, doi: 10.1051/0004-6361/202038325
  * Morley et al. (2024) Morley, C. V., Mukherjee, S., Marley, M. S., et al. 2024, arXiv e-prints, arXiv:2402.00758, doi: 10.48550/arXiv.2402.00758
  * Murray et al. (2020) Murray, C. A., Delrez, L., Pedersen, P. P., et al. 2020, MNRAS, 495, 2446, doi: 10.1093/mnras/staa1283
  * Radigan (2014) Radigan, J. 2014, ApJ, 797, 120, doi: 10.1088/0004-637X/797/2/120
  * Robinson & Marley (2014) Robinson, T. D., & Marley, M. S. 2014, ApJ, 785, 158, doi: 10.1088/0004-637X/785/2/158
  * Tamburo (2022) Tamburo, P. 2022, patricktamburo/pines_analysis_toolkit: PINES Analysis Toolkit Release v1.0.0, v1.0.0, Zenodo, doi: 10.5281/zenodo.6200054
  * Tamburo et al. (2022) Tamburo, P., Muirhead, P. S., McCarthy, A. M., et al. 2022, AJ, 163, 253, doi: 10.3847/1538-3881/ac64aa
  * Tremblin et al. (2016) Tremblin, P., Amundsen, D. S., Chabrier, G., et al. 2016, ApJ, 817, L19, doi: 10.3847/2041-8205/817/2/L19
  * Virtanen et al. (2020) Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261, doi: 10.1038/s41592-019-0686-2
  * Vos et al. (2017) Vos, J. M., Allers, K. N., & Biller, B. A. 2017, ApJ, 842, 78, doi: 10.3847/1538-4357/aa73cf
  * Vos et al. (2022) Vos, J. M., Faherty, J. K., Gagné, J., et al. 2022, ApJ, 924, 68, doi: 10.3847/1538-4357/ac4502
  * Vos et al. (2019) Vos, J. M., Allers, K., Apai, D., et al. 2019, arXiv e-prints, arXiv:1903.06691.  https://arxiv.org/abs/1903.06691
  * Vos et al. (2023) Vos, J. M., Burningham, B., Faherty, J. K., et al. 2023, ApJ, 944, 138, doi: 10.3847/1538-4357/acab58
  * Yang et al. (2016) Yang, H., Apai, D., Marley, M. S., et al. 2016, ApJ, 826, 8, doi: 10.3847/0004-637X/826/1/8
  * Zhang et al. (2021) Zhang, Z., Liu, M. C., Best, W. M. J., Dupuy, T. J., & Siverd, R. J. 2021, ApJ, 911, 7, doi: 10.3847/1538-4357/abe3fa