Atmospheric Waves Driving Variability and Cloud Modulation on a Planetary-Mass Object

Planetary-mass objects and brown dwarfs at the transition (T eff ∼ 1300 K) from relatively red L dwarfs to bluer mid-T dwarfs show enhanced spectrophotometric variability. Multi-epoch observations support atmospheric planetary-scale (Kelvin or Rossby) waves as the primary source of this variability; however, large spots associated with the precipitation of silicate and metal clouds have also been theorized and suggested by Doppler imaging. We applied both wave and spotted models to fit near-infrared (NIR), multi-band ( Y / J / H / K ) photometry of SIMP J013656.5+093347 (hereafter SIMP0136), collected at the Canada-France-Hawaii Telescope using the Wide-field InfraRed Camera. SIMP0136 is a planetary-mass object (12.7 ± 1 . 0 M J ) at the L/T transition (T2 ± 0 . 5) known to exhibit light curve evolution over multiple rotational periods. We measure the maximum peak-to-peak variability of 6 . 17 ± 0 . 46%, 6 . 45 ± 0 . 33%, 6 . 51 ± 0 . 42%, and 4 . 33 ± 0 . 38% in the Y , J , H , and K bands respectively, and find evidence that wave models are preferred for all four NIR bands. Furthermore, we determine the spot size necessary to reproduce the observed variations is larger than the Rossby deformation radius and Rhines scale, which is unphysical. Through the correlation between light curves produced by the waves and associated color variability, we find evidence of planetary-scale, wave-induced cloud modulation and breakup, similar to Jupiter’s atmosphere and supported by general circulation models. We also detect a 93 . 8 ◦ ± 7 . 4 ◦ (12 . 7 σ ) phase shift between the H − K and J − H color time series, providing evidence for complex vertical cloud structure in SIMP0136’s atmosphere.


INTRODUCTION
Brown dwarfs and planetary-mass objects, either isolated or orbiting host stars at wide separations, have complex atmospheric dynamics and chemistry.Due to their similar temperatures, masses, and chemical compositions (Burrows et al. 2001), they serve as analogs to lower mass or more closely orbiting gas giant exoplanets.
Spectrophotometric variability suggests active atmospheric dynamics in substellar objects.Silicate clouds have been found to form in early L dwarfs and thicken throughout the mid-L spectral class (Suárez & Metchev 2022).Brown dwarfs at effective temperatures (T eff ) of ∼1300 K (Kirkpatrick 2005, and references therein) transition from the mineral cloud-rich, late-L dwarfs (Tsuji et al. 1996;Burrows et al. 2006) to the relatively cloud-free, mid-T dwarfs (Burrows & Sharp 1999;Tsuji & Nakajima 2003;Knapp et al. 2004;Cushing & Roellig 2006).These L/T transition dwarfs exhibit elevated variability (Radigan et al. 2014;Radigan 2014;Eriksson et al. 2019;Liu et al. 2024).Younger brown dwarfs with lower surface gravity and planetary-mass objects also show higher variability rates than field brown dwarfs (Vos et al. 2022;Liu et al. 2024).Enhanced L/T transition variability has been historically associated with the precipitation and breakup of mineral cloud decks formed during the L class as the brown dwarf cools towards T class temperatures (e.g., Ackerman & Marley 2001;Burgasser et al. 2002;Reiners & Basri 2008).
However, Doppler imaging has provided tentative evidence for spot-like features in Luhman 16B through a surface map inferred (Crossfield et al. 2014) using maximum entropy principles.Re-analyzing the same data, Luger et al. (2021) identified similar features using the open source starry framework.These results must be interpreted with caution as maximum entropy-based Doppler imaging (see e.g., Vogt et al. 1987) does not deliver a unique solution but instead derives a map in which goodness-of-fit is balanced against the complexity of the image.Crossfield et al. (2014) further highlighted that as their spectra are sensitive to CO, the inferred maps may be affected by chemical abundance inhomogeneities.Notably, Crossfield et al. (2014) concluded that zonal bands would not be detectable based on the data precision; therefore, a comparison of spotted versus wave models was not possible at the time.Using Hubble Space Telescope (HST) spectrophotometry and the Aeolus code, Karalidi et al. (2016) identified two pairs of bright and dark features offset by 180 • on Luhman 16B, features which could fit either spotted or wave paradigms.
In this paper, we analyze multi-band (Y /J/H/K) photometry collected for SIMP0136 on two consecutive nights to learn about the source of spectrophotometric variability and atmospheric dynamics.We fit the photometry with both waves and spotted models, assess the performance of each model, and perform a model selection based on goodness-of-fit tests.This leads to insight into SIMP0136's horizontal and vertical atmospheric structure.
The paper is structured as follows.In §2, we explain our observational strategy, data reduction methods, and compare the periodic signals of SIMP0136 and our reference stars.We next detail the wave ( §3.1) and storm/spotted ( §3.2) models we use to fit the observed light curves and also explain the metrics we employ to assess model performance ( §3.3).In §4, we compare our model fits for the photometry of two consecutive nights ( §4.1 and §4.2).We discuss the implications for our research in terms of a preferred driver of variability ( §5.1), planetary-scale waves' physical nature ( §5.2), correlation between waves and cloud modulation ( §5.3), and a detected phase offset between color series ( §5.4).We summarize our findings and suggest future works in §6.

OBSERVATIONS AND DATA REDUCTION
NIR photometric observations were conducted at the Canada-France-Hawaii Telescope (CFHT) on the summit of Maunakea, Hawaii using the Wide-field InfraRed Camera (WIRCam; Puget et al. 2004).On 14 October 2012 (UT), exclusively J-band observations were col-lected, resulting in high-cadence (∆t = 0.0715 h) data.The following night, 15 October 2012 (UT), the filter was alternated between Y , J, H, and K bands, resulting in multi-color but lower cadence (∆t = 0.385 h) light curves.
We used the frames reduced by CFHT's default pipeline (iiwi version 2.1.100;Thanjavur et al. 2011) and extracted the photometry from flat-fielded and skysubtracted frames ("p" files in the cadc science archive1 ) using a fixed 2.8 ′′ aperture radius.The 50 th and 90 thpercentile seeing values were 1.18 ′′ and 1.43 ′′ .
The photometric timeseries were obtained using 12 sub-exposures with a per-sub-exposure effective exposure time of 8 s.The photometry was measured in individual sub-exposures and averaged to the values used for scientific analysis.As the sequence of sub-exposures is short (∼4 min including inter-exposure overheads), we can assume that SIMP0136 and reference stars are stable in flux within the sub-exposure sequence.We therefore use the dispersion of the 12 sub-exposure photometric measurements to determine the photometric uncertainty of their mean value.
The top two panels of Figures 1 and 2 display the raw (but normalized against each star's mean) light curves for SIMP0136 and 5 reference comparison stars for the 14 October 2012 and 15 October 2012 observations respectively.To correct for systematic effects, SIMP0136 and the comparison stars were divided by the mean normalized light curves of the comparison stars.The corrected light curves can be seen in the bottom two panels of Figures 1 and 2.
The reference comparison stars (see Table 1) were selected due to their location within the NIR detector's field-of-view, brightness, and lack of periodicity near SIMP0136's known rotational rate.As an initial test of of SIMP0136's and the comparison stars' periodicities, we performed Lomb-Scargle (L-S) periodogram analysis (Lomb 1976;Scargle 1982) using the LombScargle Python class within astropy (Astropy Collaboration et al. 2013, 2018, 2022).Figure 3 shows the results for the J-band observations on 14 October 2012.SIMP0136 exhibits two peaks above a 1% false-alarm-probability (FAP) corresponding to the rotational period (∼ 2.42 h) found in previous studies (e.g., Artigau et al. 2009) and half the rotational period (∼ 1.21 h).The reference stars do not contain signals above the 1% FAP near SIMP0136's rotational period.In the appendix, we test the reference stars for trends between airmass and seeing variations and measured brightness.We also confirm the reference stars' non-variability at periods < 1 d using external data from the Transiting Exoplanet Survey Satellite (TESS; Ricker et al. 2015).

METHODS
To model both planetary-scale atmospheric waves and storms/spots, we modify Imber (Plummer 2023(Plummer , 2024)), an open source Python code.The Imber code, the data used in this work, and a tutorial for duplicating our results and figures are openly available via GitHub2 and Zenodo3 .Imber was developed and refined in Plummer & Wang (2022, 2023).A more extensive and complete description of the package and its underlying methodology can be found in those articles.
Imber was created to analytically infer surface inhomogeneities (e.g., magnetic spots, storms, and vortices) on stars, brown dwarfs, and directly imaged exoplanets using a Doppler imaging-based technique for spectroscopic data and light curve inversion for photometry.It allows both data types to be included for an integrated multi-modal solution.Imber also includes a numerical simulation module with a full 3D grid to produce forward models of spectra, spectral line profiles, and light curves.Plummer & Wang (2022) demonstrated the numerical and analytical models produce outputs (e.g., line profiles, light curves) with residuals between the two on the order of 0.001%.
This section details the models we employ as well as the metrics by which we evaluate those models.In §3.1, we describe the wave model we use to fit the photometry.
§3.2 provides a brief description of the spotted model we implement via light curve inversion.To infer both wave and spot parameters, we employ Bayesian inference, specifically dynamic nested sampling (Skilling 2004(Skilling , 2006;;Higson et al. 2019) via Dynesty (Speagle 2020) within our framework.The inferred models are evaluated based on fit controlled for the number of free parameters (see §3.3).

Wave Model
For our wave model, we adopt the same approach as Zhou et al. (2022) based on Apai et al. (2017Apai et al. ( , 2021) ) and include a bias (C 0 ), linear term (C 1 ), and the sum of multiple (N ) sinusoidal functions: (1)  The linear term (C 1 ) accounts for a variation on time scales greater than our observation window.A i , P i , and θ i are the i th order amplitude, period, and phase.The free parameters are C 0 , C 1 , A i , P i , and θ i .Each additional wave adds three additional free parameters; therefore, the expression for the number of free parameters is thus: Uniform priors are assumed for each free parameter in the wave model.Based on SIMP0136's photometric variation of ∼5% in each NIR band, the individual wave components amplitudes (A i ) have a range of 2.5 ± 2.5%.For the period, we adjust the range of periods based on the results of the Lomb-Scargle Periodogram in §2.Phase varies 180 • ± 180 • .The bias is allowed to vary by ±0.25 and the slope by ±0.1 h −1 .

Storm/Spot Model
Imber was initially designed to both simulate spotted features' effects on observations and to infer spot parameters given a set of observational data.Here, we will briefly summarize the spotted model, but the authors refer readers to Plummer & Wang (2022, 2023) for a more rigorous explanation.
To enable computationally inexpensive Bayesian inference, the 2D stellar/substellar surface is represented by a 1D flux array.Baseline flux is modeled with a broadening kernel (Gray 2008) with linear limb darkening coefficients selected based on values derived by fitting Luhman 16B's spectral line profile (due to similar spectral type and T eff ) in Plummer & Wang (2022).Spots are modeled as Gaussian deviations to the baseline kernel to which bright spots add flux and dark spots subtract flux.The added/subtracted flux is scaled by the spot's temperature contrast and apparent size.Size is determined based on the spot radius and also latitudinal and longitudinal foreshortening due to the viewing angle (Lambert's Cosine Law).Object inclination is accounted for with rotation matrices (Euler-Rodrigues formula, Shuster 1993).The summed 1D flux at each time step is used to create the modeled light curve.
For this work, Imber was modified to allow both radius and temperature contrast to evolve in value from one rotation to the next.The code currently works by setting the rotation by which the evolution is complete, it then varies the spot parameter (radius or contrast) linearly at each time step over one full rotation, thereby, accounting for dynamic atmospheric activity.
The number of spots drives the number of free parameters for these models.Each spot has a latitude, longitude, radius, and temperature contrast.Similar to the wave models, we also include a bias to account for an unknown mean baseline flux.If spot evolution is incorporated into the model, each evolution increases the number of spot parameters.The model leads to the following expression for the number of free parameters: m = 1 + N (4 + 2E), where N is the number of spots and E is the number of spot evolutions.
Similar to the wave model, we assume uniform priors for the spotted model free parameters.Latitude and longitude priors encompass the entire sphere (±90 • and ±180 • respectively).Radius is sampled from 30 ± 30 • .Contrast varies uniformly from +1 (completely dark spot) to -1 (twice the background brightness).

Model Performance Metrics
To compare the relative merits of each model, we will use goodness-of-fit tests (χ 2 and reduced-χ 2 ) and the Bayesian Information Criterion (BIC) to provide metrics.Reduced-χ 2 is computed as where n is the number of data points (observations), m is the number of free parameters, O i and M i are the i th data points for the observation and model respectively, and σ is the photometric uncertainty.
We use the following expression to calculate BIC (Kass & Raftery 1995): For χ 2 , χ 2 ν , and BIC, a smaller value denotes a better fit.BIC and χ 2 ν account for both fit and the number of free parameters, thereby, weighting against more complex models.

RESULTS
In this section, we fit the SIMP0136 NIR photometry collected on 14 October 2012 and 15 October 2012 with both wave (see §3.1) and spotted models (see §3.2) to determine the primary driver of the planetary-mass object's spectrophotometric variability.) is shown at the top of each category; each alternate model has an associated ∆χ 2 and ∆BIC, denoting its performance with respect to the best fitting model.

High-Cadence Data
We fit the 14 October 2012, high-cadence J-band data with both wave and spotted models.For the wave models, we test N = 1 through 3 with improvement seen up to and including 3 waves.We find the 3-wave model is strongly preferred.4-wave models fail to converge, presumably due to too many free parameters (m = 14).We test 1 and 2 spot models, both with and without spot evolution.We attempt 3-Spot models but the models do not converge on a solution.
We find that SIMP0136 has a J-band peak-to-peak variability of 3.85 ± 0.14% for 14 October 2012.Here we compute the % variability based on the minimum and maximum values of the 3-wave model during the observational period.The uncertainty value is computed using the mean of the model standard deviation.The photometric variability measured for 14 October 2012 is lower than that seen the following night, suggesting a dynamic atmosphere over a relatively short time-span.
The retrieved 3-wave model contains two features that could allow for multi-rotational light curve evolution: offset phases between components and a long-term linear term (slope).The offset phases allow the superposition of each wave to create a dynamic observed light curve.The Bayesian inference also retrieves a +0.03 ± 0.01 % • h −1 linear slope term, demonstrating a gradual increase in flux and hints at dynamics on timescales longer than the period of observation.
Here we briefly highlight a few features from less preferred models.The 2-wave model (second most preferred) retrieves the same period/wavenumber as the 3-wave model's highest amplitude waves: k = 1 and 2, adding further support to the 3-wave model.The 2-spot model in which spots were allowed to evolve in size and contrast over each rotation returned both a dark spot (radius ∼20 • ) and bright spot (radius ∼10 • ) and was the third most preferred model.The spots were located at opposite polar latitudes (∼ ± 60 • ) and varied in both size and contrast over each rotation.The spotted model results will be discussed in further detail in §5.1.

Low Cadence Data
The 15 October 2012 observations provide lower cadence data with significantly fewer data points: 13, 13, 14, and 14 for the Y , J, H, and K bands respectively.The mean cadence between observations of ∆t = 0.385 h corresponds to maximum detectable wavenumbers ≲ 3 based on the Nyquist sampling criterion.For each NIR band, we are only able to constrain k = 1 and 2 wave models.3-wave models do not converge, likely due to the low cadence and small number of observations in each band.
Similar to the high cadence, 15 October 2012 data, multi-wave models outperform spotted models for each individual band and when all four NIR bands are fit simultaneously (see Table 2).The preferred 2-wave models for each NIR band are remarkably similar with each band retrieving k = 1, 2 models with periods at the rotational rate and half the rotational rate (see Table 3 and Figure 5).The phase offsets also appear similar with each component's k = 1 wave having a phase between 52  Using the 2-wave models, we compute peak-to-peak variability of 6.17±0.46%,6.45±0.33%,6.51±0.42%,and 4.33 ± 0.38% for the Y , J, H, and K bands respectively.Both % variability and uncertainty are computed as described in §4.1.The high variability (> 5%) found here is broadly consistent with previous observations of SIMP0136 (Artigau et al. 2009;Apai et al. 2013Apai et al. , 2017;;Metchev et al. 2013;Radigan et al. 2014;Wilson et al. 2014;Croll et al. 2016;McCarthy et al. 2024).
As can be seen in the bottom panel of Figure 5, each NIR band retrieval includes a negative linear slope with values of −0.14 ± 0.06 % • h −1 , −0.16 ± 0.04 % • h −1 , −0.21 ± 0.05 % • h −1 and −0.16 ± 0.05 % • h −1 for the Y , J, H, and K bands respectively.This is a steeper slope with the opposite sign as that seen in the 14 October 2012 data the night prior.The variation may lend evidence to unmodeled waves or other dynamics with timescales greater than the period of observation.
Spotted models have greater difficulty modelling the low cadence data.Evolving 2-spot models require approximately an equal number of free parameters as there are data points (and are therefore not considered) while non-evolving, 2-spot models have difficulty converging (with the exception of the composite Y + J + H + K fit which has the benefit of a higher number of observations).Both evolving and non-evolving 1-spot models retrieve polar spots (∼70 • − 80 • ) with radii of ∼30 • .Evolving 1-spot models tend to outperform non-evolving models, with the exception being the H band where spot evolution is preferred by BIC but not χ 2 ν . 5. DISCUSSION

Storms or Waves as Primary Variability Driver?
As can be seen in Table 2, wave models outperform spotted models in terms of χ 2 , χ 2 ν , and BIC.Here we seek to explore physical explanations for these results.
To further understand SIMP0136's atmospheric dynamics, we will use two planetary-scale parameters: the Rossby deformation radius (L D , e.g., Gill 1982) and the Rhines scale (L β , Rhines 1975).The Rossby deformation radius is the length at which rotational (Coriolis) effects become important (Gill 1982), and it can also be seen as the typical scale for atmospheric storms and vortices (e.g., Tan & Showman 2021b;Zhou et al. 2022).The Rossby deformation radius is computed (in km) as follows (Showman et al. 2013 where U is the flow speed, Ω is the angular rotational speed, and ϕ is latitude.At lengths greater than the Rhines scale, atmospheric structures transition from turbulent features to zonal jets as seen on Solar System planets (e.g., Cho & Polvani 1996;Showman et al. 2010Showman et al. , 2013;;Haqq-Misra et al. 2018); here the Rhines scale is computed (in km): where R is the object's radius Both inferred spots exceed the Rossby deformation radius and Rhines scale for SIMP0136, meaning the retrieved spotted models are likely unphysical.Considering the spotted model with the best χ 2 ν (see Table 2), the high-cadence J-band, 2-spot model with size and contrast evolution, we compute the Rossby deformation radius and Rhines scale at the inferred spot latitude.We conservatively assume a flow velocity (U ) of 1000 m s −1 based on the brown dwarf 2MASS J10475385+2124234's (Burgasser et al. 1999) measured wind speed of 650 m s −1 (Allers et al. 2020).We assume a planetary radius of 1.15 R J based on spectral energy distribution analysis of SIMP0136 by Vos et al. (2023).These assumptions result in L D = 0.9 • and L β = 11.8 • for latitudes of ±60  The retrieved spots' polar latitudes also argue against a spotted model explaining SIMP0136's variability.When exploring different inclinations through GCMs, objects viewed equator-on generate higher variability light curves than those viewed pole-on (Tan & Showman 2021b).For rotational periods on the scale of SIMP0136 (∼2.5 h), GCMs also demonstrate that equatorial regions have more enhanced temperature variation and cloud vertical extent than polar latitudes.This is in alignment with observations indicating that brown dwarfs' equatorial regions are more variable and redder (Vos et al. 2017(Vos et al. , 2018(Vos et al. , 2020(Vos et al. , 2022) ) and also cloudier (Suárez et al. 2023) than higher latitudes.This argument implies that polar structures similar to Saturn's polar hexagonal feature (Godfrey 1988) or Jupiter's circumpolar cyclones (Bolton et al. 2017;Adriani et al. 2018;Orton et al. 2017) are also unlikely to be the dominant driver of SIMP0136's variability.

Planetary-scale Waves in Substellar Atmospheres
The preference for atmospheric waves found in this work is in agreement with prior studies of SIMP0136 photometry (Apai et al. 2017;McCarthy et al. 2024) as well as other L/T transition dwarfs (e.g., Apai et al. 2021;Zhou et al. 2022;Fuda et al. 2024).Similar to our results in §4.1 and §4.2, peak signals corresponding to k = 1 and k = 2 waves have been previously identified for SIMP0136 (Apai et al. 2017).Based on observations over hundreds of rotations and thousands of hours, Luhman 16B has been found to also contain signals corresponding to groups of k = 1, 2 waves with three to four sine waves accounting for long-term light curve variability (Apai et al. 2021;Fuda et al. 2024).These results suggest latitudinal variation in wind speeds and atmospheric waves within a banded structure (Fuda et al. 2024).For VHS 1256 b, a 3-wave model was found to best fit observations over 2 rotations (Zhou et al. 2022).The 3-wave model was comprised of two waves (P 1 , P 2 = 18.8±0.2h, 15.1±0.2h), less than the rotational period (22.02 ± 0.04 h, Zhou et al. 2020), forming a beating pattern and a third, k = 2 wave with P 3 = 10.6±0.1 h (Zhou et al. 2022).
Brown dwarf 3D GCMs also support atmospheric waves driving variability at equatorial latitudes.For rapid rotators like SIMP0136, GCMs exhibit equatorial waves with longer zonal wavelengths and lower wavenumber values (similar to those found in §4.1 and 4.2), as well as the enhanced cloud coverage and temperature variation in their equatorial regions discussed in §5.1 (Tan & Showman 2021b).Strong evidence is also found for cloud radiative feedback-driven Kelvin waves (and more tentative evidence for Rossby waves) zonally propagating at equatorial latitudes, contributing to light curve variability (Tan & Showman 2021b).Kelvin waves move along a barrier, which in this case is formed by the equator (along which, the wave moves eastward); essentially, the forces pushing the fluid pole-ward due to a pressure gradient are balanced by Coriolis forces acting towards the equator in eastward moving fluids (e.g., Gill 1982).Rossby waves, an additional species of large-scale atmospheric wave, are driven by the conservation of potential vorticity (a fluid mechanics analog to the conservation of angular momentum) and the variation of the Coriolis parameter with latitude (Rossby 1945).
Planetary-scale waves have been observed to play an important role in the atmospheric dynamics of Jupiter.Quasi-stationary and alternating patterns of relatively  Rotations are based on a period of 2.414 h (Yang et al. 2016).(Top Row) Inferred 2-wave model fitted to observations.(Second Row) Residuals.(Third Row, First Column) Y -band, components with wave numbers (k) ≈ 1 and 2 and periods of 2.42 h and 1.27 h respectively.(Third Row, Second Column) J-band, components with wave numbers (k) ≈ 1 and 2 and periods of 2.43 h and 1.26 h respectively.(Third Row, Third Column) H-band, components with wave numbers (k) ≈ 1 and 2 and periods of 2.44 h and 1.27 h respectively.(Third Row, Fourth Column) K-band, components with wave numbers (k) ≈ 1 and 2 and periods of 2.43 h and 1.27 h respectively.(Bottom Row) Linear term data and fit.Each wave component and linear term is decomposed as described in Figure 4.
cloud-free NIR hot-spots and cooler ammonia cloudenhanced plumes have long been observed to be associated with the jet at the boundary of Jupiter's North Equatorial Belt (NEB) and Equatorial Zone (EZ) (Choi et al. 2013).These features are widely considered to be driven by a Rossby wave within Jupiter's equatorial region (Allison 1990;Showman & Dowling 2000;Friedson 2005) with the wave crests correlating to the ammonia aerosol-enhanced plumes and troughs with cloud-free regions due to condensate sublimation (de Pater et al. 2016;Fletcher et al. 2016Fletcher et al. , 2020)).

Cloud Modulation and Breakup Associated with Atmospheric Waves
Multi-band photometry offers the opportunity to explore if cloud modulation associated with planetary-scale waves (as seen in brown dwarf GCMs and observations of Jupiter) exists in SIMP0136's weather layer.If the variability from the waves is associated with silicate and metal clouds (Suárez & Metchev 2022;Vos et al. 2023;McCarthy et al. 2024), the light curve minima are where we would expect cloud coverage.In this scenario, the light curve maxima would be associated with hot spots, atmospheric areas of depleted aerosols in which light from deeper atmospheric layers is observable.We expect clouds to scatter and therefore redden light, and indeed, cloudier brown dwarfs have been found to be redder than less cloudy objects (Suárez & Metchev 2022).Light curve minima should be redder and maxima bluer if this hypothesis is true.A lack of correlation between the light curves and NIR color might support an alternate theory such as convective fingering (Tremblin et al. 2016(Tremblin et al. , 2019)).
To date, conclusive color variability has not been detected in SIMP0136 or similar L/T transition objects.Artigau et al. (2009) reported correlated J and K S band light curves ( ∆K S ∆J = 0.46 ± 0.06) for SIMP0136 which they found to best fit a scenario with dusty clouds in a clear atmosphere with two temperature components.Radigan et al. (2012) recorded comparable results ( ∆H ∆J and ∆K S ∆J ) for 2MASS J21392676+0220226 (2M2139, T1.5, Reid et al. 2008) and found the data suggested these ratios may be variable, particularly ∆K S ∆J .Using HST data, Apai et al. (2013) detected shallow J − H color variability for both SIMP0136 and 2M2139.Lew et al. (2020) found similarly negligible J − H variability across a population of L/T transition substellar objects including SIMP0136.Recently, a J − K S = 0.03 color change was reported by McCarthy et al. (2024) using NIR photometry collected with the 1.8 m Perkins Telescope Observatory, perhaps indicating that redder wavelengths should be considered for color variability studies.
Using the Y /J/H/K, 15 October 2012 photometry, we conduct a preliminary analysis to determine if the light curve minima are associated with redder NIR colors and if light curve maxima are associated with corresponding bluer colors.As a proxy for color, we use the normalized flux of the 2-wave models to create Y − J, J − H, and H − K color time series (see Figure 6).Qualitatively it can be seen that for each color series, the light curve (here we use the simultaneously-inferred composite Y + J + H + K curve) minima approximately correspond to significant dips towards redder colors while maxima correspond to bluer colors.This behaviour is more pronounced for longer wavelengths (e.g., H −K). The H −K color time series leads the composite Y + J + H + K curve to a small degree while Y − J and J − H each lag the light curve minima.We interpret these results as tentative evidence for complex vertical cloud structures with (presumably silicates and metal) clouds located at the minima and cloud-free regions existing at the maxima due to condensate sublimation.
Follow-up observations using, ideally space-based, platforms with NIR and mid-IR (MIR) spectroscopic capabilities such as the James Webb Space Telescope (JWST) would have the capability to gather broad wavelength times series (see e.g., JWST Cycle 2 GO Program 3548, PI J. Vos).With this data, broadband flux variations could be compared to variations in cloud coverage and chemical abundances, further constraining the nature of planetary-scale waves.
L/T transition objects such as VHS 1256 b, Luhman 16B (T0.5), 2M2139, and SIMP0136 have demonstrated, in general, subtler phase shifts than earlier L and later T dwarfs.Zhou et al. (2022) collected 42 h of time series observations of VHS 1256 b using the HST /Wide Field Camera 3 (WFC3) and identified no discernible phase shift between the F127M (1.270 ± 0.035 µm), F139M (1.385±0.035µm), and F153M (1.530±0.035µm) filters.Observing Luhman 16B with a 2.2 m groundbased telescope over 4 h in optical and NIR bands, Biller et al. (2013) detected significant phase offsets between K-band light curves and both H and z ′ bands.However, observing the same object over 6.5 h with HST /WFC3, Buenzli et al. (2015) found the J, H, and water bands to all be in phase.Apai et al. (2013) did not find phase shifts in either 2M2139 or SIMP0136 using HST /WFC3 G141 data over baselines of ∼ 5.9 h and ∼ 3.1 h respectively, but Yang et al. (2016) found modest phase shifts (∼30 • ) between light curves derived from HST /WFC3 G141 and Spitzer Ch.1 and Ch.2 over baselines of ∼10 h.Mc-Carthy et al. (2024) similarly found phase shifts between J and K s bands of 39.9 • +3.6 • −1.1 • using a 1.8 m ground-based telescope over a 8.5 h baseline.
Considering the color time series in Figure 6, it can be qualitatively seen that the H − K series is offset by ∼90 • from the J − H and Y − J color series, which are approximately in phase with one another.Adopting a similar approach as McCarthy et al. (2024), we use the signal function within Scipy (Virtanen et al. 2020) to determine the phase shift via cross-correlation.Performing cross-correlation on 5,000 samples from the H − K and J − H solutions provides a phase shift of 93.8 • ± 7.4 • .This 12.7σ detection provides evidence for complex vertical cloud structure between the pressure levels corresponding to these NIR bands, likely in the form of multiple cloud layers as suggested by Vos et al. (2023) and McCarthy et al. (2024).
Statistically significant phase shifts between the inferred NIR band wave components are not detected for SIMP0136 in this work (see Table 3).Within each band, the k = 1 and k = 2 components are offset by approximately 100 • − 140 • from one another.The k = 1 components each have phases ranging from ∼ 52 • − 58 • .The k = 2 components have phases ranging from ∼ 150 • to 205 • .For both k = 1 and k = 2 waves, the phase shifts between NIR bands are less than the 3σ uncertainty.Higher cadence observations in both the NIR and MIR, such as those capable by the JWST, would be able to reduce phase uncertainty and detect phase shifts at lower pressure levels (higher altitudes).Such detections would provide further collaborating evidence of cloud modulation associated with planetary-scale waves along with a more complete picture of SIMP0136's vertical architecture.

SUMMARY
To determine the driving source of spectrophotometric variability in SIMP0136, an L/T transition object, we analyze NIR photometry collected at the CFHT on the nights of 14 October 2012 and 15 October 2012.The data provides coverage for ≳5 rotations of SIMP0136.We modify the publicly available, open source Python code Imber (Plummer 2023(Plummer , 2024)), developed and honed in Plummer & Wang (2022, 2023), to fit the observed light curves with waves as well as spotted models.Here are our findings: 1.The 14 October 2012, high cadence J-band observations are best fit with a 3-wave model consisting of k = 1, 2, and 3 components with periods 2.41 h, 1.21 h, and 0.80 h (see §4.1).A linear term with a positive slope is also retrieved and indicates an increase in flux throughout observation.
2. The 15 October, low cadence Y /J/H/K light curves are each fit with 2-wave models with k = 1 and 2 components (see §4.2).Each of these components has periods approximating the rotational period and half the rotational period.In each of the inferred models, the linear term has a negative slope, demonstrating a change from the prior night's observations.
3. For the spotted models, the retrieved spot radii exceed the Rossby deformation radius (and Rhines scale for the dark spot) at the inferred latitudes and assumed mean flow speed and planetary radius, indicating such spots to likely be unphysical (see §5.1).
4. For the multi-band (Y /J/H/K) photometry, the inferred 2-wave models demonstrate correlation with shifts in color (see §5.3).The light curve minima appear to correspond to redder colors while the maxima correspond to bluer flux.This correlation is strongest for the H − K color time series.We propose this may be tentative evidence for planetary-scale waves traveling in the vertical plane of SIMP0136's atmosphere with enhanced silicate or iron cloud coverage in the wave crests and depleted aerosols in the troughs.
5. We detect a 93.8 • ± 7.4 • (12.7σ) phase offset between the H − K and J − H color series, providing evidence for complex vertical cloud structure in SIMP0136's atmosphere (see §5.4).
Moving forward, a greater understanding of the true nature of planetary-scale waves potentially driving brown dwarf and planetary-mass object variability can be achieved by multi-rotational and quasi-simultaneous observations in the NIR and MIR by platforms such as JWST.Spectroscopic modes would help to discern if the reddening associated with NIR light curve minima is tied to variations in cloud coverage or chemical abundances.Both photometric and spectroscopic modes could further constrain phase shifts between NIR and MIR bands and color indices, providing a greater understanding of SIMP0136's vertical structure throughout its atmospheric layers.
the lightkurve package (Lightkurve Collaboration et al. 2018).These observations provide insight into the variability of our reference stars.As shown in Figure 7, none of the stars display short-timescale percent-level variability, and they do not display significant periodicity above the 1% FAP level for periods shorter than 1 day.
We investigated possible correlations between observing conditions and reference star fluxes.Figures 8 and 9 illustrate the fluxes measured as a function of seeing and airmass respectively.There is an anticorrelation of flux with both quantities in both day 1 (J-only) and day 2 (interleaved Y JHK) datasets.
The purpose of reference stars is to account, at least to first order, for these extinction effects; we, therefore, showed each reference star dependency against seeing and airmass once corrected by the sum of all other reference stars (right panels of Figures 8 and 9).No trend is seen after correction by reference stars at this step and we conclude that reference stars are, for all practical purposes, stable to well within 1%.
The dependency of raw fluxes against airmass can easily be understood as an increase in atmospheric extinction at higher airmass while the dependency with seeing probably arises from a combination of aperture losses as well as the covariance of seeing with airmass (i.e., a worsening seeing correlates with higher extinction as both happen at a higher airmass).Figure 10 shows the changes of seeing, median sky level and airmass through our datasets.The expected mild correlation of seeing and sky background is seen mostly on day 1.Day 2 observations are shorter and explore a smaller airmass range.   .Each reference star flux has been normalized to its median; one sees a ∼ 6 % loss in flux between measurements with a seeing of ∼0.9 ′′ and∼1.6 ′′ .To assess if differential biases were present in our photometric measurements, we compared each of the reference stars to the average of all stars (right) against seeing.The calibrated fluxes show no dependency against seeing, which confirms that our reference stars will serve their purpose in calibrating SIMP0136 timeseries.

Figure 2 .Figure 3 .
Figure 2. Low Cadence, Y /J/H/K photometry of SIMP0136 collected at the CFHT on 15 October 2012.(Top) Raw normalized light curves of SIMP0136.(Second Row) Raw normalized light curves of comparison stars.(Third Row) Corrected light curves for SIMP0136.(Bottom) Corrected light curves for comparison stars.

Figure 4 .
Figure 4. High cadence, 14 October 2012, J-band observations and 3-wave model.Shaded regions are 1σ uncertainty.Rotations are based on a period of 2.414 h (Yang et al. 2016).(Top) The inferred 3-wave model fitted to observations.(Second Row) Residuals.(Third Row) 3-wave components with wave numbers (k) ≈ 1, 2, and 3 and periods of 2.41 h, 1.21 h, and 0.80 h respectively.To demonstrate each component wave's fit, the remaining waves and long-term linear terms are subtracted from the data.(Bottom) Linear term with wave components subtracted from the data.

Figure 5 .
Figure 5. Low cadence, 15 October 2012, Y /J/H/K-band observations and 2-wave models.Shaded regions are 1σ uncertainty.Rotations are based on a period of 2.414 h(Yang et al. 2016).(Top Row) Inferred 2-wave model fitted to observations.(Second Row) Residuals.(Third Row, First Column) Y -band, components with wave numbers (k) ≈ 1 and 2 and periods of 2.42 h and 1.27 h respectively.(Third Row, Second Column) J-band, components with wave numbers (k) ≈ 1 and 2 and periods of 2.43 h and 1.26 h respectively.(Third Row, Third Column) H-band, components with wave numbers (k) ≈ 1 and 2 and periods of 2.44 h and 1.27 h respectively.(Third Row, Fourth Column) K-band, components with wave numbers (k) ≈ 1 and 2 and periods of 2.43 h and 1.27 h respectively.(Bottom Row) Linear term data and fit.Each wave component and linear term is decomposed as described in Figure4.

Figure 6 .
Figure 6.Color time series comparison to composite Y + J + H + K light curve (gray).Shaded regions are 1σ uncertainty.Vertical dashed lines denote color series and light curve minima.Qualitatively, it can be seen that the composite light curve minima correspond to reddening (i.e.low flux in Y − J, J − H, and H − K corresponds to redder color indices).This provides evidence for cloud scattering (see §5.3).Light curve maxima correspond to bluer, presumably relatively cloud-free regions.The Y − J and J − H color series are approximately in phase but the H − K solution is ∼90 • out of phase with the remaining color series, offering evidence of SIMP0136's complex vertical cloud structure.

Figure 7 .
Figure7.TESS photometric lightcurve for the 5 reference stars (left).None of the stars display short-term (< 1 day) variability above a 1% false alarm probability and no signal is seen close to the period of SIMP0136..

Figure 8 .
Figure8.Correlation of normalized raw fluxes of reference stars against seeing (left).Each reference star flux has been normalized to its median; one sees a ∼ 6 % loss in flux between measurements with a seeing of ∼0.9 ′′ and∼1.6 ′′ .To assess if differential biases were present in our photometric measurements, we compared each of the reference stars to the average of all stars (right) against seeing.The calibrated fluxes show no dependency against seeing, which confirms that our reference stars will serve their purpose in calibrating SIMP0136 timeseries.

Figure 9 .
Figure 9. Same as Figure 8 but for the dependency of fluxes (raw and corrected) with airmass..

Figure 10 .
Figure10.Variation of the seeing, median sky background and airmass through the J-band observing night.Sky background, and to a lesser extent seeing, are correlated with the airmass through the observing sequence.

Table 2 .
Model Comparison: χ 2 , BIC, observations (n), and free parameters (m) . We can see the dark spot's inferred radius range exceeds L β and both spots' radii are far greater than L D .It is also unlikely that a group of adjacent spots of equivalent summed size are responsible for the observed variability as the same arguments above would require such spots to be in separate •