© The Authors 2025

1 Introduction

The broad wavelength coverage and superb stability of JWST has enabled transformative studies of the dynamic processes in brown dwarf and exoplanet atmospheres (Biller et al. 2024; McCarthy et al. 2025). In the absence of stellar contamination present when studying their young, directly imaged exoplanet cousins, isolated brown dwarfs can be observed with exquisite spectrophotometric precision. For nearby brown dwarfs, this allows us to measure changes in their emission spectrum over time. This variability is due to changes in brightness over the 2D emitting surface of the object, which are thought to be caused by changing cloud coverage, thermal structure, chemistry, and the presence of upper-atmosphere aurorae (Radigan et al. 2014; Robinson & Marley 2014; Tremblin et al. 2020; Tan & Showman 2021; Faherty et al. 2024). Spitzer and Hubble observed light curves for dozens of brown dwarfs, finding that the variability amplitude depended on the spectral type, age, and viewing geometry (e.g. Metchev et al. 2015; Yang et al. 2016; Vos et al. 2017; Lew et al. 2020; Vos et al. 2022). This variability monitoring enabled the first steps towards three-dimensional (3D) studies of substellar atmospheres, allowing for the identification of bands and spots (Apai et al. 2017), as well as pressuredependent weather (Apai et al. 2013; Yang et al. 2016). However, it is only in the era of JWST that broad wavelength, spectroscopic light curves have been observed, allowing us to relate the dynamic atmospheric processes to changes across the entire emission spectrum (Biller et al. 2024; McCarthy et al. 2025; Chen et al. 2025).

SIMP J013656.5+093347.3 (hereafter SIMP-0136) is a nearby (6.12±0.02 pc, Gaia Collaboration 2020), well-studied exoplanet analogue, ideally suited for variability studies due to its brightness, large variability amplitude of about 2.5%, and 2.4 hour rotation period (Artigau et al. 2006, 2009). With a spectral type of T2.5±0.5 (Artigau et al. 2006) and an effective temperature of around 1100 K (Zhang et al. 2020), SIMP-0136 shares many similar atmospheric properties to directly imaged exoplanets, such as HR 8799 b or AF Lep b (Marois et al. 2008; Balmer et al. 2025). As a member of the Carina-Near Association, SIMP-0136 is estimated to be 200±50 Myr old (Gagné et al. 2017). Likewise, its relatively low surface gravity also implies a relatively youthful object, with log g = 4.46 ± 0.09 (Zhang et al. 2020). Its mass has been estimated to range from 12.7 ± 1 M_Jup (Gagné et al. 2018) to 17.8 ± 11.9 M_Jup (Vos et al. 2023). Its inclination was measured by Vos et al. (2017) to be nearly equator-on, with $i={80}_{-12}^{+10}{}^{\circ }$. In the most in-depth atmospheric analysis of SIMP-0136 to-date, Vos et al. (2023) performed atmospheric retrievals on a broad wavelength spectrum consisting of Spex PRISM (Burgasser et al. 2008), Akari (Sorahana & Yamamura 2012), and Spitzer/IRS (Filippazzo et al. 2015) data. They identified layers of patchy silicate and iron clouds and robustly inferred the atmospheric temperature structure and composition. The presence of an aurora has been identified as a mechanism to explain the observations of pulsed radio emission (Kao et al. 2016, 2018). The interaction of the aurora with the upper atmosphere has also been suggested to contribute to the observed near-infrared (NIR) variability (Morley et al. 2014; McCarthy et al. 2025).

In the NIR, the variability amplitude is mildly wavelengthdependent, with a weaker amplitude in the 1.4 µm water absorption feature (Apai et al. 2013; Lew et al. 2020). McCarthy et al. (2024) found a phase shift of ${39.9}_{-1.1}^{+3.6\circ }$ between the J and Ks band light curves. Yang et al. (2016) demonstrated that the different wavelength ranges correspond to different pressures probed by the emission spectrum. Both groups attributed the chromatic variability effects to differing cloud coverage at different levels of the atmosphere, while Apai et al. (2017) and Plummer et al. (2024) determined that planetary scale (e.g. Kelvin or Rossby) waves are likely the primary driver of the variability. Recently, McCarthy et al. (2025) (hereafter M25) published the first JWST observations of SIMP-0136, from 0.8–11 µm using NIRSpec/PRISM and MIRI/LRS. The authors proposed that the distinct variability patterns could be tied to variations in the cloud patchiness, high-altitude hotspots, and changes in the carbon chemistry.

Given this observed wavelength dependent variability, it is crucial to quantify how changes in the atmospheric state drive the observed changes in the emission spectrum. In the era of JWST, atmospheric retrievals (e.g. Madhusudhan & Seager 2009; Benneke & Seager 2012; Mollière et al. 2019; Hood et al. 2023; Grant et al. 2023; Kühnle et al. 2025; Hoch et al. 2025) have become one of the primary tools for characterising the atmospheres of exoplanets and brown dwarfs. Atmospheric retrievals typically use a 1D parametric forward model to fit the observed spectrum. They are highly dependent on the wavelength coverage of the data, as different sources of opacity only impact the spectrum at particular wavelengths. Burningham et al. (2021) and Vos et al. (2022) demonstrated the importance of broad wavelength coverage, combining NIR and MIR observations to robustly characterise silicate clouds.

To date, most retrieval studies have treated extrasolar objects as one-dimensional (1D). However, from 3D general circulation models (GCMs), it is clear that exoplanets and brown dwarfs vary across their surface and as a function of depth in the atmosphere. Blecic et al. (2017) highlighted the implications of this for atmospheric retrievals of transiting planets, showing that at best retrievals will find the average thermal structure of the object. Some attempts have been made to extend retrieval frameworks to 3D for transiting planets, which exhibit strong day-to-night contrasts due to the incident stellar radiation (e.g. MacDonald & Lewis 2022; Dobbs-Dixon & Blecic 2022). Without any incident stellar light, isolated brown dwarfs do not exhibit such diurnal contrasts and the variation in the emission spectrum across the surface is typically small, of the order of a few percent at most (Showman & Kaspi 2013; Tan & Showman 2021; Lee et al. 2024). In this work, we extended the petitRADTRANS retrieval framework (Mollière et al. 2019; Nasedkin et al. 2024a) to perform time-resolved retrievals of emission spectra, similar to Gandhi et al. (2025), to tie the atmospheric properties to the observed variability.

This study is structured as follows. In Section 2, we outline the JWST observations, as well as our novel approach to processing the time-series data. Section 3 presents the details of our retrieval framework, presenting the parameterisations used for the thermal structure, chemistry, and clouds. We present the results of this study in Section 4, highlighting the measured atmospheric properties of SIMP-0136 and how they vary over time. These results are placed in the context of existing theoretical models and JWST observations in Section 5 and the conclusions are summarised in Section 6. Additional information is included in the appendices, including a comparison of SIMP-0136 to self-consistent forward models and evolutionary models (Appendix A), an in-depth demonstration of the sensitivity of the retrievals to various mechanisms driving the variability (Appendix B) and additional results from the retrieval analysis (Appendix C).

2 Observations

We obtained time-resolved spectra of SIMP-0136 using NIR-Spec/PRISM and the MIRI/LRS as part of GO 3548 (PI: Vos). These were originally presented in M25, and the observing setup is summarised in Table 1. A total of 5726 spectra were obtained with NIRSpec/PRISM (Jakobsen et al. 2022) in bright object time-series (BOTS) mode. These observations were made with a cadence of 1.8 s spanning a total of 3.4 hours, or 1.4 rotations (P=2.4 hours, Yang et al. 2016). These spectra have a spectral resolving power ranging from R = 30-400, and cover 0.65.3 µm, and are shown in Fig. 1. The MIRI/LRS (Kendrew et al. 2015) observations were taken subsequently and 575 spectra were obtained with a cadence of 19.2 s, over a total of 3.55 hours of observations. The LRS observations have a spectral resolving power varying from R = 40-160 across the wavelength range from 5–14 µm (also presented in Fig. 1). Similar to Fig. 4 of M25 we phase-fold the LRS spectrum to align with the NIRSpec observations to produce a complete spectrum that varies with the rotation of SIMP-0136.

2.1 Data reduction

The JWST NIRSpec/PRISM data were processed using the standard stage 1 of the JWST pipeline (v1.16.1, Bushouse et al. 2025) for the bright-object time series (BOTS) mode and a customised version of the stage 2 (calwebb_detector2) pipeline, which was optimised for time-series spectral extraction and the mitigation of systematic noise sources. The reduction used CRDS version 12.0.6 with the jwst_1298.pmap context file. The output of our modified pipeline was the flux-calibrated, time-series spectra.

Unlike McCarthy et al. (2025), who used the default stage 2 reduction, we modified stage 2 to process each integration individually, resulting in better correction of systematics. By default, the BOTS reduction does not apply the nsclean step to correct for 1/f noise, which manifests as low-frequency spatial structure in the detector readout (Rauscher 2024). Since this correction was applied to fixed-slit observations, we modified the EXP_TYPE metadata from NRS_BOTS to NRS_FIXEDSLIT and applied nsclean to each integration. The standard stage 2 BOTS reduction assumes a fixed background and detector properties and applies single correction factors to each. To avoid these biases, we applied flat-field and photometric calibrations to each spectrum individually. By default, the pipeline assumes one fixed extraction aperture and spectral trace across all integrations, which does not account for time-dependent shifts due to optical distortions, telescope jitter, or variations in prism dispersion. These shifts can introduce flux losses and wavelength misalignment if a static extraction aperture is used, and so our modified pipeline automatically determines a spectral trace for each integration. We fit a first-order polynomial to the trace positions across the detector columns for each individual exposure, improving both wavelength calibration and flux accuracy. This corrects for systematics introduced by time variable detector drifts, pointing variations, and changes in pixel sensitivity. We found that these improvements resolved systematic flux variations between integrations, independent of wavelength. These are likely caused by systematic detector effects such as small gain or bias drifts, charge trapping, or variations in readout electronics. The default pipeline’s application of a static flat-field correction, as well as a uniform spectral trace, can introduce artificial time-dependent offsets if pixel sensitivity fluctuates over the course of the observation. With the spectra extracted, we removed low signal-to-noise (S/N) data from wavelengths shorter than 0.9 µm and longer than 5.3 µm. Three wavelength bins near 1.58 µm were contaminated by a hot pixel, and were also removed. We binned the time-series spectra to highlight the astrophysical variation of SIMP-0136 over its rotation. The bins were centred at every 15° of rotation, and combined 25 spectra on either side of these bin centres, for a total of 90 s of integration time into each bin. These binned integrations both improve the S/N by a factor of √50 compared to a single exposure, but are also sufficiently short to prevent the blending of the spectrum as SIMP-0136 will only rotate by 3.75° in 90 s.

The MIRI/LRS observations were reduced using v1.16.1 of the JWST pipeline using the settings described in M25. As with the PRISM data, we removed low S/N data at wavelengths shorter than 5.3 µm and longer than 10.8 µm. This range also ensures that there is no overlap with the NIRSpec data and it extends out to cover the 10 µm ammonia absorption feature. Additional wavelength bins near 6.0 µm and 6.68 µm were also removed due to poor hot-pixel correction. While we excluded data below 5.3 µm, the LRS spectrum was offset by an average of 2% below the PRISM spectrum between 5.1 and 5.3 µm. This is compatible to within the PRISM uncertainties, but the three overlapping LRS data points range between fully compatible to offset by 10σ. The data were binned similarly to the PRISM data; however, only two integrations were combined for the LRS, for a total of 38 s of integration time. This produces a S/N across the spectrum that is proportional to the emitted flux and matches the S/N at 5.3 µm in the binned PRISM spectrum. These bins were aligned such that they were in phase with the NIRSpec bins. The mean-normalised variability maps for both NIRSpec and the LRS are presented in Fig. 1.

Fig. 1 NIRSpec/PRISM (left) and MIRI/LRS (right) observations of SIMP-0136. Top: emission spectra in 15° rotation bins. Each binned NIRSpec spectrum is averaged over 90 s of integration, while each LRS spectrum is averaged over 38 s integrations to ensure that the overall signal to noise of the input spectra for the retrievals is proportional to the emitted flux. Centre: uncertainties calculated by the JWST pipeline for each spectrum in blue. Standard deviation of the 90 s (115s for the LRS) intervals in red, demonstrating that the pipeline uncertainties accurately reflect the underlying statistical distribution of the noise. Bottom: binned variability maps, also known as dynamic spectra. These variability maps are aligned in phase, as the MIRI observations were taken after the NIRSpec observations and are binned to highlight the wavelength dependence of the variability. The 24 phase bins used in these maps correspond to the binned spectra used as inputs for the retrievals.

2.2 Uncertainty estimation

The high time resolution of the observations allowed us to validate the uncertainties produced by the JWST pipeline. As shown in Fig. 1, we grouped the time-series spectra into bins of 90 s for NIRSpec and 115 s for the LRS, with the same 15^o cadence as the binned spectra, taking the standard deviation of the spectra within these bins. This timescale is short enough to ensure only small astrophysical changes in the spectra and, thus, the dominant variation between the samples is random noise. We compared these measurements of the random noise to the uncertainties estimated by the JWST pipeline and found that this estimate closely matches the measured noise as shown in the centre panels of Fig. 1.

In contrast to known under-estimation of uncertainties with the MIRI/MRS (e.g. Miles et al. 2023; Patapis et al. 2024; Matthews et al. 2025), the uncertainties associated with NIR-Spec/PRISM and MIRI/LRS time-series data accurately reflect the underlying noise distribution. As the noise is normally distributed, we can use standard error analysis to combine N integrations, reducing the uncertainty by a factor of $1/\sqrt{N}$. In addition to the statistical uncertainties, both NIRSpec/PRISM and the LRS have systematic uncertainties in the absolute flux calibration. The PRISM spectrum is expected to have an absolute flux calibration accurate to within 5%, while the LRS spectrum calibration can vary between 3–8% (JWST User Documentation 2016; Gordon et al. 2022). We did not account for these systematic biases in our uncertainty estimates, although we note that offsets in the absolute flux will lead to biased measurements, particularly in terms of the luminosity and effective temperature.

3 Atmospheric modelling

To determine the mechanisms underpinning the observed variability, we must be able to accurately infer the atmospheric state as a function of time. We first compared the SIMP-0136 spectrum to self-consistent, radiative-convective equilibrium forward models to fit the spectrum, the results of which are discussed in Appendix A. These models did not fit the spectrum sufficiently well to resolve the observed level of variability, motivating the use of a data-driven retrieval approach to infer the atmospheric properties. In addition, we estimated the fundamental parameters of SIMP-0136 using SEDkit (Filippazzo et al. 2015), the results of which are included in Table 2. These measurements were used to provide a baseline against which to compare the results of the atmospheric retrievals.

Overall, petitRADTRANS (pRT) is a widely used, open source tool for computing emission and transmission spectra (Mollière et al. 2019; Blain et al. 2024). It combines a flexible, 1D parametric forward model with a Bayesian framework (Nasedkin et al. 2024a), which has been used to characterise the atmospheres of directly imaged exoplanets (e.g. Mollière et al. 2020; Nasedkin et al. 2024b; Landman et al. 2024; Hoch et al. 2025) and brown dwarfs (e.g. de Regt et al. 2024; Zhang et al. 2024a; de Regt et al. 2025) via atmospheric retrievals of their emission spectra. These retrievals are enabled by the use of the correlated-k method for radiative transfer (Goody et al. 1989; Lacis & Oinas 1991), allowing efficient computation of the emission spectra up to a spectral resolving power of 1000, and the Feautrier method for calculating scattering in the atmosphere in the presence of clouds (Feautrier 1964). We used pRT version 3.2.0 to perform the retrievals in this study, using exo-k (Leconte 2021) to bin the correlated-k tables to R=400 for model computation. As in previous studies, we used an adaptive pressure grid, using a total of 134 pressure levels between 10³ and 10⁻⁶ bar, with a 10× finer grid spacing at the location where the clouds are found to condense. To validate the effectiveness of pRT in identifying different potential mechanisms for variability, we produced a set of forward models in Appendix B. This allowed us to demonstrate that different processes, such as changes in the temperature structure, cloud coverage, and chemical profiles impact the spectrum in uniquely identifiable ways, reproducing the self-consistent modelling analysis of Morley et al. (2014).

The pRT retrieval module combines the forward model with Multinest (Feroz & Hobson 2008; Feroz et al. 2019) implemented in pyMultinest (Buchner et al. 2014) to sample the parameter space, estimate the model evidence, and infer posterior parameter distributions. Unless otherwise noted, we use the constant efficiency mode of Multinest with a sampling efficiency of 5% and 400 live points to reduce the computation time of the retrievals to make retrieving time-resolved spectra feasible. When combined with importance nested sampling, this reduces the accuracy of the evidence estimate, making model comparison unreliable (Feroz et al. 2019), although the parameter estimation remains minimally affected. A full listing of the priors used in the fiducial retrieval setup is available in Table 3.

Each model was convolved with a Gaussian kernel calculated from the wavelength-dependent instrumental spectral resolving power, before being binned to the wavelength bins of the data. The spectral resolution functions for NIRSpec/PRISM and MIRI/LRS are obtained from the JWST CRDS system (Greenfield & Miller 2016). The log-likelihood is then calculated to measure the goodness-of-fit between the data, D, model, F, and covariance matrix, C, expressed as $$-2\log \mathcal{L}={(\mathbf{D}-\mathbf{F})}^T{\mathbf{C}}^{-1}(\mathbf{D}-\mathbf{F})+\log (2\pi det\left(\mathbf{C}\right)).$$(1)

For the SIMP-0136 data, the covariance matrix is diagonal and consists only of the uncertainties associated with each wavelength bin. Given the broad wavelength coverage and large numbers of pixels sampled, correlations due to inter-pixel cross talk will not impact the overall slope of the spectrum, and will likely only impact a few neighbouring wavelength channels. Likewise, as the SIMP-0136 spectra were obtained using slit spectroscopy with nearly straight traces spread over many pixels, relatively few interpolations were made during the data reduction that would impart additional correlations into the data. The assumption of independent Gaussian uncertainties is therefore a reasonable approximation for this dataset, although future analyses could include fitting for a covariance matrix. We did not use any error inflation in the retrievals presented in this work, although the impacts of including error inflation are discussed in Section 5.1. We also tested the addition of an offset between the PRISM and LRS data. While a small offset was found, it did not significantly impact the retrieved parameters.

3.1 Temperature profile

Our fiducial temperature parameterisation is based on that of Zhang et al. (2023) (hereafter, Z23). This approach has become widely adopted for sub-stellar atmospheric retrievals; for example, in its application in Zhang et al. (2024b) or Nasedkin et al. (2024a). In the Z23 setup, the retrieved parameters are the i^th temperature gradients (d log T/d log P|_i) at a set of log-spaced pressure nodes, with the temperature profile calculated by integrating the temperature gradient and interpolating onto the full pressure grid. The temperature at the bottom of the atmosphere (T_bot) was also retrieved. Thus, for i increasing from the bottom of the atmosphere, we have $$\begin{array}{rl}& T_0=T_{\mathrm{B}\mathrm{o}\mathrm{t}},\end{array}$$(2) $$T_{i+1}=\exp (\log T_i+(\log P_{i+1}-\log P_i){\left(\frac{d\log T}{d\log P}\right)}_i).$$(3)

Throughout this work, we use the notation $$\mathrm{\nabla }{T}_i\equiv {\frac{d\log T}{d\log P}}_i.$$(4)

In contrast to the original implementation of this profile, we retrieved the temperature gradient from the bottom to the top of the atmosphere, without fixing an isothermal region at pressures lower than 10⁻³ bar.

In Z23, the priors on the temperature gradient parameters were determined by measuring the temperature gradients across a set of radiative-convective equilibrium models. By using relatively narrow Gaussian priors centred at these values, the retrieval was constrained to physically plausible solutions. We explored the impact of varying the prior width on the retrieved temperature profile, finding no significant variation and confirming that the posterior distributions are not determined by the prior widths at pressures between 10 bar and 0.1 mbar. To allow the fits to be driven primarily by the data, we used a modestly wider set of priors, as listed in Table 3.

To validate our choice of temperature profile, we performed test retrievals using the profiles of Mollière et al. (2020) and Guillot (2010). The more flexible Mollière profile verified that the temperature structure did not depend on the choice of parametrisation, finding qualitatively similar results. However, the Guillot profile, which was only allowed to monotonically decrease with decreasing pressure, was insufficiently flexible to fit the data, producing reduced χ² values 10 to 100 times worse than with the Z23 profile.

3.2 Chemistry

Through a small set of test retrievals, we found that a free-chemistry approach produced significantly better fits than one relying on equilibrium chemistry, with disequilibrium abundances for H₂O, CO, and CH₄, as was used in Mollière et al. (2020). In particular, the equilibrium approach was not able to reproduce the CO₂, H₂S, and PH₃ abundances. We therefore freely retrieve abundance profiles for H₂O (Polyansky et al. 2018), CO (Rothman et al. 2010), CO₂ (Yurchenko et al. 2020), CH₄ (Guest et al. 2024), C₂H₂ (Chubb et al. 2020), HCN (Barber et al. 2014), FeH (Wende et al. 2010), H₂S (Azzam et al. 2016), NH₃ (Coles et al. 2019), PH₃ (Sousa-Silva et al. 2015), Na (Allard et al. 2019), and K (Allard et al. 2016). The mass fraction abundances of each of these molecular species were retrieved as free parameters, and the remaining atmosphere was filled with a mixture of 76% H₂ and 23% He by mass, similar to the composition of Jupiter (Atreya et al. 2003; Ben-Jaffel & Abbes 2015). In addition to the impact of the molecular species on the spectra from their line absorption and emission, we also account for Rayleigh scattering from H₂ and He, as well as the collisionally induced absorption for the H₂–H₂ and H₂ –He interactions.

Previous studies have shown that retrievals are sensitive to changes in chemical abundances with pressure. For instance FeH in L-dwarfs (Rowland et al. 2023; de Regt et al. 2025) and H₂O in Y-dwarfs (Kühnle et al. 2025). Rowland et al. (2023) also suggested that non-uniform CH₄ abundances will be important for interpreting the emission spectra of L-T transition objects, such as SIMP-0136. In the presence of upper atmosphere heating from aurora or other phenomena, this effect will likely be magnified as CH₄ becomes disfavoured in equilibrium compared to CO. Therefore, we also developed a parametrisation for variable abundance profile throughout the atmosphere. The abundance at the top and bottom of the atmosphere are freely retrieved parameters. Between these boundaries, pressure-temperature pairs are also retrieved. This setup allows the retrieval of the location where changes in abundance occur, as well as the abundance at that position. The location of the pressure node is determined by retrieving the difference in log pressure between the node and the previous node, measured from the top of the atmosphere down. Between these nodes, the abundance profile can be interpolated using a linear spline, as illustrated in Fig. B.2. For the present study, we used a single pressure-temperature pair as free parameters for the CH₄ abundance to measure potential changes induced by the thermal structure of the atmosphere and to better fit the depths of CH₄ absorption features in the spectrum.

3.3 Clouds

Clouds are a critical component of the atmosphere, contributing a source of continuum opacity and reddening. We used a cloud model inspired by Vos et al. (2023); similarly to their retrievals of SIMP-0136, our fiducial model uses an amorphous Mg₂SiO₄ (Servoin & Piriou 1973) cloud together with a crystalline Fe cloud (Henning & Stognienko 1996). The cloud opacities are calculated using a distribution of hollow spheres approach (Toon & Ackerman 1981; Min et al. 2005), similarly to Mollière et al. (2020). We also explored varying the cloud composition and number of cloud layers. We tested replacing the iron cloud layer with crystalline KCl (Palik 1985) and Na₂S clouds (Morley et al. 2012), which are expected to condense at photosphere pressures in T-dwarfs, while the iron cloud should condense too deep to be observable. As McCarthy et al. (2024) found multiple layers of silicate clouds, we also tested the inclusion of patchy Mg₂SiO₄ and MgSiO₃ (Jaeger et al. 1998), as well as both silicate clouds and a full-coverage Na₂S cloud. As the clouds in SIMP-0136 are expected to occur near the base of the photosphere or deeper, we did not expect the results to depend strongly on the details of the cloud composition or particle geometry. Thus, we leave a more detailed study of the cloud properties to a future analysis.

We retrieved the mean particle radius of each layer, together with the width of a log-normal distribution about the mean. The cloud base pressure for each cloud was freely retrieved, along with the mass fraction at the base of each cloud layer. We find that this mass fraction decreases with decreasing pressure as X_i(P) = X_i,0(P/P₀)^f_sed, where f_sed is used as a parameter to determine the rate of the abundance decrease. As in Vos et al. (2023), the cloud coverage fraction of the silicate cloud, f_cloud was also allowed to vary while the non-silicate cloud is assumed to enshroud the entire object. From test retrievals, we found that the patchiness of the silicate cloud and the mass fraction were moderately degenerate; occasionally, the retrieval would find a solution where the cloud coverage approached 1, while the cloud mass fraction approached 0. On average, the silicate cloud mass fraction at the cloud base was found to be between 10⁻³ and 10⁻⁴ by mass. As we did not expect variations in the cloud abundance by any significant order of magnitude, we set the lower limit of the mass fraction prior of the silicate cloud to 10⁻⁶, thus avoiding the degenerate solution.

4 Results

4.1 Setting mass and radius

To observe the variability of atmospheric parameters, it is necessary to first constrain the properties of SIMP-0136 that cannot vary as a function of phase; namely, its mass and radius. We ran a small set of retrievals at each of 0°, 90°, 180°, and 270° phase, where mass and radius were included as free parameters. This experiment was repeated using both a uniform prior distribution on the mass and radius, and with a Gaussian prior based on evolutionary expectations. For the uniform case, the radius prior was set to 𝓊(a = 0.7, b = 3.0) R_Jup and the mass to 𝓊(3,18) M_Jup, while for the evolutionary prior case, the priors were set based on the measurements of Gagné et al. (2017) 𝒩(μ = 1.1, σ = 0.1) R_Jup and 𝒩(12,1) M_Jup. We also performed an additional set of retrievals using the SEDkit measurements as priors on the mass and radius. Notably, this enforced a very strong prior on the radius of 𝒩(μ = 1.2, σ = 0.05) R_Jup. This resulted in a radius more consistent with the evolutionary models, but resulted in unphysical cloud and temperature structure parameters; thus, we discarded this set of models. Overall, this is related to the well-known ‘small radius problem’, discussed in detail in Balmer et al. (2025). There are known inconsistencies between the radii measured using atmospheric and evolutionary models, possibly due to an incomplete treatment of the clouds. As such, for the remainder of this work, we treat the mass and radius as nuisance parameters for obtaining the surface gravity, which has a more direct impact on the emission spectrum.

Both the uniform and Gagné et al. (2017) prior setups produced qualitatively similar atmospheric states. The resulting masses and radii at each phase are presented in Table 4. The evolutionary prior retrievals produce marginally better fits for 3/4 phases, as measured by the reduced χ². This set of retrievals inferred masses closer to expectations from evolutionary models, along with measured log g values consistent with literature measurements. Thus, we were able to take the average of this set of retrievals as our fixed mass and radius, excluding the outlying measurement of 7.7 M_Jup. The mass was therefore fixed at 10.381 M_Jup and the radius at 0.9329 R_Jup, resulting in a log surface gravity of 4.490. The surface gravity was kept consistent with the SEDkit results, while the mass and radius were both set to be lower, which we attribute to the small radius problem.

4.2 Atmospheric properties of SIMP-0136

Although we performed an ensemble of retrievals to understand the observed spectroscopic variability, it is insightful to consider the bulk planet properties obtained from a single retrieval. We discuss the retrievals of the 0° observations as a representative example, with the best-fit model shown in Figure 2. From the emission contribution function, shown in Fig. 3, we see that the emission spectrum probes four decades in pressure, from 10 bar where the spectrum is cut off in the NIR by clouds, to methane features at 3.3 μm and ~7 μm at pressures as low as 10⁻³ bar. This implies that the retrievals are sensitive to the thermal structure and chemical abundances across this pressure range, while there is no information imparted onto the observed spectrum from both deeper pressures and at very high altitudes.

4.2.1 Stratospheric heating

The retrieved temperature profile is presented in Fig. 4. In the troposphere, at pressures greater than 0.1 bar, the temperature profile follows an adiabat, similar to the temperature profiles predicted by radiative-convective equilibrium forward models. There is no evidence of an isothermal region, such as those produced by fingering convection (Tremblin et al. 2015, 2016), which would act to redden the spectrum without requiring clouds. Between 10⁻¹ bar and 10⁻³ bar the temperature gradient inverts, and begins increasing with increasing altitude, reaching a maximum near 2 × 10⁻³ bar. The transition from an adiabatic profile in the troposphere to an inverted temperature profile in the stratosphere at 0.1 bar is qualitatively similar to the structure of Jupiter’s atmosphere, where the tropopause is also located at 0.1 bar (Gupta et al. 2022). This is clearly in contrast with the self-consistent forward models, which are usually monotonically decreasing with altitude in the absence of external irradiation. Compared to a temperature profile that becomes isothermal at the point of the inversion, the hotspot at 0^° phase peaks at 265 K warmer at 3 × 10⁻³ bar. From the forward models of Appendix B.1, and from the emission contribution function, shown in Fig. 3, we found that this region is probed by the 3.3 µm and 7.7 µm CH₄ features. This is the strongest source of opacity in the atmosphere, and probes the lowest pressures. We performed two additional retrievals to assess the importance the CH₄ features in shaping the inversion. These retrievals were run while masking out the 3.1–3.5 µm feature, and masking out both 3.1–3.5 µm and 7–8 µm. In both cases, the strength of the inversion is greatly reduced, with the temperature profile weakly oscillating around an isothermal solution at pressures lower than 10⁻² bar. Without any strong opacity source we could not reliably infer the temperature structure at pressures lower than 1 mbar.

We calculated the derivative of the retrieved temperature profile, d log T/d log P, and compared this to the adiabatic temperature gradient to determine where the atmosphere is convectively unstable, following the Schwarzchild criterion. The adiabatic temperature gradient was calculated for an atmosphere with a metallicity of 0.2 and a C/O ratio of 0.55. It was found using the interpolated easyCHEM chemical equilibrium grid (Mollière et al. 2017). Both temperature gradients are shown in Appendix C in Fig. C.1. At 17 bar, the retrieved temperature gradient is consistent with an adiabatic gradient, but is constrained by the prior distribution. At pressures lower than 5 bar the atmosphere is stable against convection, and is therefore radiative. Based on the contribution function, this transition coincides with the top of the silicate cloud layer.

At 0° phase, we measured the effective temperature to be 1246.82 ± 0.02 K by integrating a low resolution model between 0.5 µm and 250 µm. Given the fixed radius, this is simply a function of the bolometric luminosity of SIMP-0136 at 0° phase, which is measured to high precision due to the precise spectrophotometry and wavelength coverage of the combined NIRSpec/PRISM and MIRI/LRS observations. Comparing the full wavelength range of the model to the wavelength range of the data, we found that the JWST observations measured 95.7% of the total emitted flux. We compared the T_eff found by integrating the low resolution model to that obtained by directly integrating the observed spectrum. There was significant systematic discrepancy between the methods: as several wavelength bins are excluded from the data due to hot pixels, the flux is somewhat underestimated, and the effective temperature was found to be 1197 K after applying the bolometric correction factor. This is around 100 K warmer than expectations from evolutionary models or SED fitting (see Appendix A), but is also about 80 K colder than found by Vos et al. (2023), which we attribute to their smaller retrieved radius. While the mean T_eff varied between methods, the relative variation was found to be consistent, varying by 5 K between the coldest phase and the warmest.

Fig. 2 Best-fit model from the fiducial retrieval to the observed spectrum at 0° phase, whose uncertainties are smaller than the displayed line width. This model assumed a fixed mass and radius of 10.381 MJup and 0.9329 RJup respectively. Also shown for visual reference are the log opacities at the 1 bar level of the primary absorbers in each wavelength range. The CH4 abundance was allowed to vary as a function of altitude. The reduced χ2 of the t s T. However, the model clearly matches the key absorption features of H2O, CH4, CO, and CO2, as well as additional species such as H2S and NH3.

Fig. 3 Emission contribution function for SIMP-0136 at 0° phase. The emission contribution is shown in square-root colour scale to highlight subtle variations. Significant contributions are present across four decades in pressure, bounded by a cloud layer near 10 bar and by methane absorption up to 10−3 bar.

4.2.2 Atmospheric chemistry

We were able to measure the abundance of 10 of the 12 chemical species included in the free retrieval, placing upper limits on the abundance of acetylene (C₂H₂) and phosphine (PH₃), with combined abundance profiles from each phase for each species shown in Fig. 5. Full posteriors are available in the appendix in Fig. C.3. Due to the computational cost of retrievals when not using Multinest’s constant efficiency mode, we were unable to reliably measure the Bayesian evidence and calculate detection significances for each species. Of the species with constrained abundances, H₂O, CH₄, CO, CO₂, and NH₃ have clearly visible absorption features present in the spectrum, as seen by comparing the opacities to the spectrum in Fig. 2. In particular, H₂ S is consistently a strong source of opacity in the NIR, with its opacity reaching a maxima between 2.5 µm and 3 µm, as well as between 3.5 µm and 4.1 µm, although it lacks uniquely identifiable absorption lines in the spectrum. As the retrievals fit these features well, with only small variations in the abundances with phase, the abundances of these species are reliably measured. HCN, FeH, Na, and K all contribute to the continuum opacity. The low spectral resolution of NIRSpec/PRISM prevents the identification of the alkali lines in the near-IR, and so the abundances of Na and K are largely constrained by their impact on the NIR slope, which is degenerate with the impact of clouds in this region. FeH acts similarly, and its low abundance indicates that it does not strongly impact the spectrum. Such a low abundance is consistent with all of the iron condensing out deep within the atmosphere.

The retrieved abundances of these species are mostly consistent with an equilibrium chemistry model with [M/H]=0.2 and C/O=0.55, suggesting that the retrievals have inferred plausible abundances for each of these species. Both the retrieved abundances and the equilibrium profiles are shown in the appendices in Fig. C.3. The H₂O abundance agrees well with the equilibrium abundance at a pressure of 2 bar. Then, CH₄ and CO are consistent near 2 bar, where they appear to quench. H₂ S is strongly enriched compared to the equilibrium prediction. Na is very consistent up to 1 bar. NH₃ is somewhat under abundant, but intersects with the equilibrium profile around 0.05 bar. K and HCN are both consistent up to 1 bar. Both FeH and C₂H₂ drop off to 0 at around 2 bar in equilibrium, consistent with the upper limits set by the retrievals; PH₃ is severely under abundant compared to equilibrium expectations, although this is consistent with the missing phosphine problem widely observed in brown dwarfs Visscher et al. (2006); Rowland et al. (2024).

To test the assumption of chemical equilibrium, we ran a series of retrievals allowing different levels of flexibility in the CH₄ and CO abundance profiles. As described in Section 3.2, we compared using a vertically constant abundance profile, a profile with a single freely retrieved pressure-abundance pair, and a profile with three retrieved pressure-abundance pairs. Using vertically constant abundance profiles, we found that the retrievals were unable to capture both the 1.2 µm and 3.3 µm CH₄ absorption features, and the goodness-of-fit was poor, with a reduced χ² of ~360. Adding a node to the CH₄ abundance profile provided a significantly better fit to the data, with ∆χ²∕ν ≈ 80 and a physically plausible temperature profile. While allowing additional flexibility in the CH₄ profile led to significantly better fits, adding additional nodes to the CO profile only marginally improved the reduced χ², implying that a vertically constant abundance is sufficient for CO. With a single node, we find that (as shown in Fig. 5) the CH₄ abundance decreases sharply above 2 × 10⁻² bar. Lastly, when including three nodes for each of CH₄ and CO, we find that the profiles match that of the single-node retrievals. Even with sufficient flexibility, the abundances do not reproduce the equilibrium chemistry profiles, instead suggesting that both CH₄ and CO are in chemical disequilibrium. This exercise was repeated, allowing H₂O to vary as well. Adding the additional flexibility to fit the water abundance profile resulted in a worse reduced χ²/ν, suggesting that the inclusion of these additional parameters is not justified by the data. Regardless of the degree of flexibility used, the overall interpretation of the atmosphere remains qualitatively similar, with only minor variation in the inferred atmospheric metallicity.

As both CH₄ and CO appear to be in chemical disequilibrium, we can apply the approach of Zahnle & Marley (2014) to derive the vertical mixing strength, K_zz, required to quench both species. We calculated the K_zz from the quench point of CH₄ and CO. CH₄ intersects the equilibrium profile at 17 bar, near the cloud base, while CO intersects the equilibrium profile at about 2 bar. Calculating K_zz using the temperature at the quench point, around 1500 K at 2 bar or 1850 K at 17 bar, then we find a log K_zz between 6-9 log cm² s⁻¹, which is relatively consistent with the GCM predictions of a log K_zz ~ 9 at the temperatures of SIMP-0136 (Lee et al. 2024). However, the value of K_ZZ depends strongly on the temperature, which varies throughout the atmosphere. If instead we use the effective temperature instead of the temperature at the quench point, we find a range of log K_zz = 3.3-4.6 in units of log cm² s⁻¹, depending on the choice of whether the methane or CO quench pressure is used. This relatively low value is similar to the mixing strength found for HR 8799 b (${2.9}_{-0.7}^{+0.6}$ log cm² s⁻¹, Nasedkin et al. 2024b), which has a similar spectral type to SIMP-0136 (Bonnefoy et al. 2014). This orders-of-magnitude variation in the mixing strength suggests that non-vertically-constant Kzz parameterisations are necessary to interpret the chemical state of the atmosphere.

Fig. 4 Retrieved temperature profiles for each phase bin. Temperature profiles are given as 90% confidence intervals. In grey is the emission contribution function, averaged over wavelength. We show the square root of he contribution to highlight that there is small, but significant emission from the upper atmosphere, coincident wit the location of the temperature inversion. In black lines are condensation curves for enstatite, forsterite, and iron, showing that the retrieved cloud base pressure for the silicate cloud is coincident with the expected condensation location W also show the temperature profile from self-consistent ExoRem model (log g = 4.5, [M∕H]= 0, C/O= 0.55, Teff = 1000 K) as a red dotted line.

Fig. 5 Volume mixing ratio profiles of gas-phase species for all phases. The shaded regions indicate 68% confidence intervals, taken from the set of median abundances measured at each phase. The measured precision for a single phase is typically an order of magnitude smaller than the variation between phases.

4.2.3 Patchy silicate clouds

While there is no distinct silicate feature identified in the MIR, as seen in the case of late-L dwarfs (e.g. Suárez & Metchev 2022, 2023), silicate clouds are found to contribute strongly to the shape of the emission spectrum. At 0°, the forsterite cloud base is located at 13.10 ± 0.03 bar, forming the base of the photosphere and strongly absorbing at wavelengths shorter than 2 µm. The condensation curves shown in Fig. 4, taken from Visscher et al. (2010), show that the silicate clouds are expected to condense between about 8 bar and 15 bar given the inferred temperature structure. At this depth, the silicate clouds do not significantly contribute their own MIR absorption features, as the MIR photosphere lies above the cloud layer. At the cloud base, the log mass fraction of Mg₂SiO₄ is found to be −4.33 ± 0.01. The f_sed value, which determines the rate at which the silicate abundance decreases with altitude, was found to be $f_{\mathrm{s}\mathrm{e}\mathrm{d},{\mathrm{M}\mathrm{g}}_2{\mathrm{S}\mathrm{i}\mathrm{O}}_4}=9.90\pm 0.05$, consistent with the upper limit of the prior, and creating a very compact cloud layer. Micron-sized particles were found to best fit the observed cloud properties, with a narrow log-normal distribution ($\log r_{c,{\mathrm{M}\mathrm{g}}_2{\mathrm{S}\mathrm{i}\mathrm{O}}_4}=-3.994\pm 0.003,{\sigma }_{\mathrm{L}\mathrm{N}}=1.19\pm 0.01$. The silicate clouds were found to cover 87.14 ± 0.08% of the atmosphere at this phase.

In addition to the forsterite clouds, an additional cloud is necessary to fit the observations; several different compositions were tested. In our fiducial model, we included a crystalline iron cloud, which was found to be optically thin and condensed at a slightly lower pressure (i.e. higher altitude) than the forsterite cloud. Based on the condensation curves, iron would be expected to condense at pressures greater than 30 bar, significantly deeper than the silicate clouds. When the iron cloud was replaced with an Na₂ S cloud, the condensation location was found to remain constant, without significant changes in the goodness-of-fit. As with the iron clouds however, the location of the Na₂ S cloud did not correspond with its expected condensation pressure. A KCl cloud was also tested, and resulted in a marginally worse fit than either the iron or Na₂ S clouds. The addition of an additional patchy enstatite cloud also had no significant impact on the goodness-of-fit of the spectrum. When combined, this suggests that in addition to a patchy silicate cloud, an additional source of grey opacity near the base of the photosphere near 10 bar is required to fit the data. This cloud may be compositionally different from the silicate cloud, or it may be due to the cloud parametrisation failing to accurately model the silicate clouds.

Fig. 6 Top: phase-resolved best-fit spectra divided by the mean retrieved spectrum for NIRSpec/PRISM (left) and the MIRI/LRS (right), taken from the fiducial retrieval. Bottom: difference between the observed and retrieved variability maps. The remaining difference is normally distributed, centred at 0, and with a standard deviation of 0.3% for NIRSpec and 1% for MIRI.

4.3 Time-resolved atmospheric retrievals

We ran independent retrievals on each of the 24 phase-binned spectra of SIMP-0136, using the same setup as described in Section 4.2. We refer to this setup as the fiducial, or baseline setup for the retrievals. Across both the NIRSpec and MIRI spectra the models are consistent to within 15% of the observed flux at every wavelength bin; the RMS deviation between the NIR-Spec data and the models is 2.60 ± 0.08%, while for MIRI it is 3.8 ± 0.2%. Across all phases, the largest discrepancy occurs between 5.86 µm and 6.3 µm, where the models find stronger water absorption than is present in the data, and the onset of the water absorption feature at 6.4 µm is slightly redshifted in the model compared to the data. For each retrieval, the residuals between the models and the data were dominated by systematic effects which were consistent across phase, as shown in Fig. C.4 in Appendix C. These discrepancies may be caused by missing physical processes in the models, inaccurate line-lists for the molecular species, or unknown issues in the observations or data processing. However, these systematic effects are effectively removed when dividing by the mean spectrum to examine the variability trends, resulting in typical discrepancies between the mean-normalised models and data of less than 0.3% across the NIRSpec range and <1% across the MIRI data, as shown in the lower panels of Fig. 6.

The retrievals were able to accurately fit the spectra, and were able to reproduce the morphology of the observed spectroscopic variability, as shown in Fig. 6, which can be compared to the observed variability maps of Fig. 1. For example, the brightness between 2 μm and 3.5 μm decreases between 100° and 200° of rotation, while shorter wavelengths display a higher frequency phase variation. Fig. 6 also shows the difference between the observed and modelled variability maps. In general, the residuals are normally distributed, and centred at 0, suggesting that there are few remaining systematic biases in the models. With a standard deviation of 0.3%, the residuals between the data and the models in the NIRSpec range are much smaller than the observed variability of 3%, which implies that the models are accurately capturing the varying atmospheric state.

4.3.1 Fixed parameter retrieval tests

To determine the impact of specific mechanisms driving the variability, we performed a series of tests, fixing different sets of parameters to the median retrieved values and observing if the retrievals retained sufficient flexibility to reproduce the observations.

For the first test, we fixed the temperature profile to the median profile across all phases, leaving the rest of the setup the same as in the baseline model. As shown in the top panels of Fig. 7 this setup is unable to reproduce the observed variability maps. The overall fits are poor and the structure of the variability does match that of Fig. 1. It is clear that allowing the temperature profile to vary is critical to reproducing the variability, and that the variability cannot be explained solely through changes in the cloud coverage and chemistry.

In the second test, we fixed the cloud parameters to the median retrieved values from the baseline retrieval. In contrast to fixing the temperature profile, this setup is able to reproduce the observed variability better than the baseline retrieval, as shown in the bottom panels of Fig. 7. This reinforces the findings of the fixed-temperature test, implying that changes in the thermal structure are driving the observed variability. While patchy clouds are necessary to fit the spectrum, variations in the cloud coverage or opacity are not required to explain the spectroscopic variability. This model provides a better match to the observed variability map, with a standard deviation in the residuals of 0.24% for PRISM (compared to 0.3% for the fiducial model) and 1% for the LRS (compared to 1.03%). Crucially, the remaining residuals are smooth, without the visible features between 1 to 1.5 µm and at 3.3 µm, suggesting that the fixed cloud model provides a better measure of the CH₄ abundance than the fiducial model. As such, we continue to use this model throughout this section to observe parameter variations as a function of phase, and to produce the figures in this work unless otherwise noted. For completeness, the phase variation of each parameter is presented in Appendix C in Fig. C.2.

Fig. 7 Same as Fig. 6, but with sets of parameters kept constant across phase. Top: variability map keeping the temperature profile constant, showing that the retrievals are unable to reproduce the observed variability. Bottom: variability map keeping the cloud parameters constant, which better reproduces the observed variability than the full retrieval.

Fig. 8 Fig. 8. Measurements of key atmospheric parameters as a function of phase, taken from the fiducial retrieval. Also indicated are the median and �1σ confidence intervals for each set of measurements.

4.3.2 Statistical tests

The retrieved posterior distributions, both in terms of median values and posterior widths are relatively consistent over the course of the observations, with only a few outliers. In Fig. 8, we show the inferred posterior distributions at each phase for a set of key atmospheric parameters. The inferred posterior width is significantly smaller than the observed variation in the parameters. To determine whether the measured variation in the atmospheric parameters is physical, or simply due to scatter in the measurements, we perform runs tests on each set of parameter measurements (Wald & Wolfowitz 1940). We apply the runstest_1samp function implemented in the statsmodels package, correcting for small sample statistics (Seabold & Perktold 2010). This test determines the probability 1 - p that a given series of runs (sequences of consecutive observations either greater than or less than the mean value) are drawn from a normal distribution. The sequence is considered to be distinct from random draws if p < 0.05.

Across the ensemble of retrievals, CO₂ is the only parameter that is reliably found to follow a trend with p_CO2 = 0.002. For most of the retrieval sets, H₂S is also found to be non-random, with p_H2S ≈ 0.002; however, in the fiducial case, this is reduced to p_H2S ≈ 0.14.The fiducial model is the only setup where p_H2S is not significant. Both CO₂ and H₂ S follow similar trends in each set of retrievals, increasing until around 180° before decreasing towards 360°. The effective temperature is also found to be nonrandom, with p_Teff = 2 × 10⁻⁵.

In addition to the non-parametric runs test, we also compare each parameter variability curve to the clustered light curves of M25, and to each other. Using the clustering analysis of Biller et al. (2024), the authors identified nine unique light-curve patterns in the NIRSpec data and a further two in the LRS data. These clusters were found to correspond to distinct pressures of the atmosphere. The combination of pressure and wavelength was linked to major opacity sources when compared to a Sonora Diamondback radiative-convective equilibrium forward model (Morley et al. 2024). In the NIRSpec data, cluster one was assigned to the lowest pressure, 0.12 bar, while cluster nine was assigned to the deepest pressure at 10.44 bar. The light curves curve clusters could broadly be grouped into two categories. Light curves one through five show approximately sinusoidal variation with phase, while clusters six through nine have a double-peaked shape. This suggests that these groups of light curves are sensitive to different mechanisms driving the variability, with different physical processes causing the different characteristic shapes.

The correlation between the retrieved parameters as a function of phase and the clustered light curves is illustrated in the top panels of Fig. 9. For each pair we calculate the Pearson correlation coefficient, r. An r value of -1 indicates that the pair are fully anti-correlated, while a value of 1 indicates full correlation. Given the 24 measurements for each parameter and setting a threshold of 90% confidence to determine significance, we find using a 2-sided t-test that the threshold for significant correlation is |r| > 0.271. The full set of r values between the parameters and light curves is shown in Fig. 9. Most parameters are uncorrelated with the light curves. The strongest correlation is with the effective temperature, which is strongly correlated with the light curves from 0.66 bar and deeper in the atmosphere. CO₂ is strongly anti-correlated with the light curves, with the strongest anti-correlation with light curve 6. M25 found that light curve 6 probes 0.66 bar in pressure, and was also co-located with the 4.2 µm CO₂ feature.

We also compare the correlations between the atmospheric parameters, shown in Fig. 10. Many of the parameter correlations are due to degeneracies in the model. For example, the parametrisation for the CH₄ abundance profile can set the pressure node at either the top or bottom of the photosphere to achieve the same overall abundance gradient. Likewise, the silicate cloud mass fraction and the cloud coverage fraction are anti-correlated, as a decrease in the cloud coverage can be compensated for by an increase in the cloud mass fraction.

However, some correlations seem to be astrophysical in nature. The strongest correlation between chemical parameters and temperature profile parameters is the correlation between the methane abundance parameters and the temperature gradient at the location of the temperature inversion, ∇T₆. As methane is an effective thermometer for the atmosphere of SIMP 0136, this correlation demonstrates that the observed temperature inversion is likely found because of the measured methane abundance.

Fig. 9 Top: comparison of normalised effective temperature and the ∇T4 parameter with the clustered light curves of M25. The effective temperature is strongly correlated with the light curve, while ∇T4 is anticorrelated. Bottom: matrix of Pearson correlation coefficients between atmospheric parameters and clustered light curves from M25. Nonsignificant correlations (|r| < 0.271) are not shown.

Fig. 10 Pearson correlation coefficients between the phase-variation of each pair of retrieved parameters in the fiducial retrieval. Using a two-sided t-test with 24 samples, the minimum statistically significant correlation at a 90% confidence level is 0.271, thus setting the limits on which correlations are shown. While some correlations are natural degeneracies of the model, such as the correlation between ∆ logPCH4 and the CH4 abundance at that pressure node, other pairs are likely astrophysical in nature.

4.3.3 Thermal structure evolution

The temperature profile in the troposphere (between 0.1 bar and 10 bar) is retrieved with high precision (±0.5 K per temperature point) at every phase, shown in Fig. 4. In this region the temperature changes by ±3 K between phases per pressure level. This level of temperature variation is consistent with the ~1 K temperature variations predicted by Showman & Kaspi (2013) for the top of the convective region in a rapidly rotating object. For SIMP-0136 we find the radiative-convective boundary to occur at 17 bar, as shown in Fig. C.1. At pressures between 17 bar and 0.1 bar the temperature profile is moderately sub-adiabatic, and is broadly consistent with the ExoRem temperature profile of an 1100 K object in radiative-convective equilibrium. Below the base of the photosphere at 17 bar the temperature profile is essentially unconstrained, with the width of ∇ log T₁ and ∇ log T₂ posteriors determined by the prior distributions.

Every retrieved temperature profile shows evidence of a thermal inversion in the stratosphere (Fig. 4), beginning at around 40 mbar, and reaching a maximum near 3 mbar. The amplitude of the inversion varies more strongly with phase than the deep atmosphere, with a standard deviation of 40 K between phases at the peak of the inversion. The amplitude of the inversion peak does not vary smoothly with phase, but is generally the strongest around 200° and weakest near 345°, although the actual largest amplitude hot spot occurs at 0°. The location of the inversion coincides with the top of the photosphere, as measured by the τ = 1 surface in the 3.3 μm CH₄ feature. As in the case of the single retrieval, there is no contribution at pressures lower than 1 mbar, and therefore this region is also unconstrained. At the top of the atmosphere, the uncertainties reach about 100 times larger than the uncertainties in the adiabatic region, with the standard deviation between phases reaching 70 K.

Using the same statistical tests as applied to the inferred atmospheric parameters, we test the measured temperature at each pressure for randomness using the runs test, and correlate the temperature variations with the clustered light curves of M25, which is shown in Fig. 11. None of the individual temperature variation curves were found to be statistically different from random draws. However, the effective temperature is found to be significantly different from random, with $p_{{\mathrm{T}}_{\mathrm{e}\mathrm{f}\mathrm{f}}}=2\times {10}^{-5}$. The effective temperature varies smoothly over the surface of SIMP-0136, with hemisphere-averaged temperature variations of about 5 K as shown in Fig. 12.

While the correlations between the temperature profile and the clustered light curves at any individual pressure level are generally weak, there are clear trends between different atmospheric regions. The PT profile is correlated with light curves 1 through 5 at all pressures below the stratosphere, with statistically significant anti-correlations in the inversion and near the top of the troposphere, as shown in Fig. 11. This implies that at least some of the observed variability in light curves one through five can be attributed to changes in the strength of the inversion, and its resulting impact on the opacities in this region of the atmosphere. There are transitions in the behaviour of the correlations both between the troposphere and stratosphere, and in the region of the silicate clouds. We observe a general correlation between the tropospheric temperature profile and the effective temperature at both a few hundred mbar and deep in the atmosphere at a few bar. As the bulk of the observed emission is from the troposphere, we attribute the primary variation in the luminosity with the change in temperature in this reason, consistent with the ad hoc approach of Tremblin et al. (2020). These changes may be due to cloud radiative feedback Tan & Showman (2021) driving temperature fluctuations, or by the turbulent atmosphere distributing energy inhomogeneously (Showman & Kaspi 2013).

Conversely, the stratospheric inversion is anti-correlated with the light curves associated with molecular absorption features, and is only weakly correlated with the effective temperature. This is likely due to the strength of the inversion being probed by the methane absorption at 3.3 µm and around 7 µm, acting as a thermometer for the upper atmosphere. While we cannot attribute a causal link between the changing stratospheric temperature and the methane absorption, we propose that mechanism may be changing auroral strength, which we discuss further in Section 5.2.

Fig. 11 Correlation between the phase-variation of the temperature profile and the clustered light curves from M25. Top: example temperaturevariation curves drawn from three representative pressure levels in the atmosphere (highlighted boxes) compared to correlated or anticorrelated light curve clusters. Bottom: correlation coefficient between the temperature-variation at that pressure and the clustered light curves for each pressure level, as well as the correlation to the effective temperature in the right-most column. The pressure levels for the temperature curves in the top panels are highlighted with coloured boxes. Horizontal lines and annotations indicate the location of atmospheric features.

Fig. 12 Fig. 12. Effective temperature for each observed phase. Over one rotation Teff varies by 5 K.

4.3.4 Chemical variation with phase

Of all the measured species, only CO₂ and H₂S are found to vary smoothly as a function of phase. Both of these species are anticorrelated with the effective temperature, with r_CO2 = −0.69 and r_H2S = −0.41. However, in absolute terms both species only vary by small amounts, with CO₂ varying by about ±2% about its median value, while the precision on each measurement is about 0.06%. H₂S varies by ±4% about the median value, with a precision of 0.03% on each measurement. This level of variation is comparable to the scatter seen in other species: H₂O varies by about ±1%, albeit without any clear trend identified as a function of phase. Without such a trend, we conclude that the abundances of most species, particularly H₂O, CH₄, and CO are homogeneous over the surface of SIMP-0136, to within the scatter of the measurements, or about 1%. In contrast, the abundances of CO₂ and H₂S are likely determined by the temperature of the atmosphere, increasing in abundance as the temperature falls. This pattern, particularly the sharp transition in the H₂S abundance in the fixed cloud retrieval (Fig. C.2), is qualitatively similar to the correlation between the chemical abundances and gas temperature predicted in Lee et al. (2024). This suggests that the changes in the chemical abundances may be due to a storm rotating into view, bringing localised cloud, chemistry, temperature variations.

By measuring the abundances of the primary carriers of carbon, oxygen, sulphur, and nitrogen, we are able to calculate the atmospheric metallicity, as well as the elemental abundance ratios relative to the solar value (Lodders et al. 2009). As the primary molecular species do not show trends with time, these elemental ratios can be measured by combining the abundance measurements at each phase. We find that the metallicity is slightly enriched compared to the solar value, with $\left[\mathrm{M}/\mathrm{H}\right]={0.184}_{-0.001}^{+0.006}$. The [C/H] ratio is slightly more enhanced than the oxygen ratio, with $\left[\mathrm{C}/\mathrm{H}\right]={0.218}_{-0.002}^{+0.006}$ and $\left[\mathrm{O}/\mathrm{H}\right]={0.140}_{-0.001}^{+0.003}$. [S/H] as measured by the H₂S abundance is strongly enhanced compared to solar, with $\left[\mathrm{S}/\mathrm{H}\right]={0.783}_{-0.003}^{+0.013}$, while the nitrogen abundance is very strongly depleted $\left[\mathrm{N}/\mathrm{H}\right]=-{1.85}_{-0.02}^{+0.01}$. However, the nitrogen abundance is measured only through NH₃ and HCN, and does not account for N₂, which should be the dominant carrier at the temperatures of SIMP-0136 (Zahnle & Marley 2014). From these elemental abundances, we can also calculate additional elemental abundance ratios. We find that the gas phase C/O ratio in the troposphere is modestly super-solar, with C/O = ${0.658}_{-0.002}^{+0.004}$, although this does not account for oxygen sequestration in the silicate clouds, which would reduce the bulk C/O ratio by around 20% (Lodders 2003; Calamari et al. 2024). Accounting for this, the C/O ratio is 0.54 ± 0.01, consistent with the solar value. Both values are significantly lower than the ratio found by Vos et al. (2023), where C/O= 0.79 ± 0.02. Similarly, we find the atmospheric S/O=0.12 ± 0.00 and C/S=5.60 ± 0.15. The elemental abundance ratios are shown in Fig. 13. Biazzo et al. (2012) measured metallicities for low-mass stars in the Carina-Near Association, of which SIMP-0136 is a member. While they do not measure elemental abundance ratios for volatile species, the typical metallicities for these objects is [Fe/H]=0.08 ± 0.05, suggesting that SIMP-0136 is only modestly enriched in metals relative to objects with a similar origin. Compared with its siblings, our metallicity measurement is more compatible with the expectations for a low-mass brown dwarf than previous metallicity measurements for SIMP-0136 (Line et al. 2017), as in Vos et al. (2023) who found a metal-poor atmosphere with $\left[\mathrm{M}/\mathrm{H}\right]=-{0.32}_{-0.05}^{+0.05}$.

Fig. 13 Elemental abundance ratios derived from retrieved molecular abundances combined from all phases.

4.3.5 No change in patchy clouds with phase

From the runs tests there is no evidence to support systematic variation of the cloud coverage as a function of phase p_fcloud > 0.05, reinforcing that variations in the clouds do not drive the observed variability for this early T-dwarf. Likewise, the cloud coverage is not found to be correlated with any of the clustered light curves. The cloud parameters are highly (anti)correlated with each other, and are weakly correlated with some temperature and chemical parameters. While the patchy clouds are required to reproduce the observed emission spectrum, changes in the cloud coverage do not appear to drive the observed spectroscopic variability. Finally, Fig. 7 shows that the fixed-cloud set of retrievals provided a better fit to the variability map than the fully free retrievals, highlighting that changes in the cloud parameters are not necessary to reproduce the observed variability. This is consistent with the picture of Freytag et al. (2010), where the cloud patchiness is driven by gravity waves, resulting in patchiness on small spatial scales.

5 Discussion

5.1 Error inflation

In contrast to previous studies (e.g. Vos et al. 2023), the JWST data in this work were obtained with much higher S/𝒩, and developments in atmospheric models allow us to more accurately fit the broad wavelength spectrum. However, the same S/𝒩 that allows us to measure changes in the spectra over time also lays bare the limitations of our 1D parametric models: typical reduced χ² values of retrievals in this study are ~300, with systematic residuals of the order of 5% between the models and data. The typical solution to this is to add an error-inflation term to the retrieval, and increase the uncertainties until a reduced χ² of 1 is achieved. However, we found this would increase the uncertainties, and by extension the widths of the posterior distributions, to the point where we were no longer sensitive to the observed level of variability.

Error inflation is a standard technique for inferring uncertainties on the data when they are difficult to directly measure, as is the case for the MIRI/MRS (e.g. Kühnle et al. 2025; Barrado et al. 2023; Matthews et al. 2025), as well as for accounting for the uncertainties in the models. By comparing the stochastic variations in the observed spectrum to the uncertainties calculated by the JWST pipeline in Fig. 1, it is evident that the uncertainties are well understood. However, in Fig C.4 there are strong systematic residuals present between the data and the best fit models, indicating that the models poorly describe the observations. This is consistent with the large reduced χ² values, which demonstrate poor goodness-of-fit.

To compensate for these model uncertainties, we performed an additional set of retrievals using error inflation based on the approach of Line et al. (2015), where $${\sigma }_{\mathrm{T}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l}}=\sqrt{{\sigma }^2+{10}^b}.$$(5)

Here, the total uncertainty σ_Total is found by adding an additional inflation term to the measured uncertainties on a per-instrument basis, where b is a freely retrieved parameter. We found that the uncertainties were substantially inflated for both NIRSpec and MIRI to compensate for the systematic deviation in the models, with b_PRISM = −31.41 ± 0.02 and b_LRS = −33.765 ± 0.03, and a resulting χ²/ν of ∼1 for each retrieval performed. The typical uncertainties σ on the NIRSpec PRISM data are of the order of 10⁻¹⁶ –10⁻¹⁷ W/m²/µm. This implies that the total uncertainty is dominated by the 10^b term, as σ² + 10^b ≈ 10^b, resulting in a total uncertainty between 3–30 times greater than the measured uncertainties.

However, the inflation of the uncertainties on the posterior parameter distributions due to the increased measurement uncertainties resulted in an inability to distinguish variations as a function of phase. Key parameters such as the cloud coverage fraction and the cloud base pressure were unconstrained. The median values of these parameters were also less physically consistent than in the non-inflated retrievals: for example, in the fiducial retrievals the silicate clouds are found to condense near 10 bar, at the pressure expected from the condensation curves. In the inflated case, the silicate clouds condense at 0.01 bar, far above where they would be expected in an object with the effective temperature of SIMP-0136.

Based on the inability to reproduce the observed variability, measure variations in atmospheric parameters, and the physical implausibility of the results of the inflated retrievals, we only consider the non-inflated case in presenting results for SIMP-0136. Although there remain significant systematic discrepancies between the models and the data, the magnitude of these variations is typically small, and they do not impact the ability to reproduce the observed variability. Future model development incorporating 3D atmospheric dynamics, improved line lists, and better cloud parameterisations is likely necessary to further reduce these discrepancies.

5.2 Auroral heating in upper atmosphere

SIMP-0136 is known to emit pulsed radio signals (Kao et al. 2016, 2018), which is attributed to the large scale currents generated by electrons flowing along magnetic field lines. These electrons interact with the atmosphere, releasing energy and forming aurorae. While our retrievals find heating in the stratosphere of SIMP-0136, it is challenging to directly attribute this to heating from auroral processes: there is still no direct evidence for UV or IR auroral emission in SIMP-0136. Characteristic signatures of aurorae such as H₃+ (Gibbs & Fitzgerald 2022; Pineda et al. 2024) emission have not been observed, although higher resolution NIRSpec/G395H observations obtained with GO 5814 (PI: Fitzgerald) may place meaningful constraints on this emission. Likewise, there no observations of UV emission from any brown dwarf aurora available at present (Saur et al. 2021). Without such direct measurements, we can explore whether auroral activity could reasonably explain the observed temperature inversion and place this in context with auroral observations from within the solar system.

The stratospheric heating appears to be present at all phases, and while we cannot resolve latitudinal differences, must cover a significant fraction of the observed hemisphere. The maximum difference in the stratospheric temperature from an isothermal model is found to be about 265 K at 3 mbar, which is much stronger than the 37 ± 4 K heating at 1 mbar in Jupiter (Rodríguez-Ovalle et al. 2024), but comparable to the 300 K auroral heating at 1–10 mbar on WISE J1935 (Faherty et al. 2024). However, the mechanism through which the atmosphere is heated by an aurora is unclear. In Jupiter, the hotspot in the outer reaches (at pressures as low as 10⁻⁶ bar) of the atmosphere is thought to be heated by electron precipitation (Sinclair et al. 2018). The mechanism for heating Jupiter’s second hotspot, located at 1 mbar, is poorly understood (Rodríguez-Ovalle et al. 2024). Multiple mechanisms have been proposed, including adiabatic heating from local downwelling initiated by electron precipitation, and radiative heating of stratospheric aerosols produced by precipitation (Sinclair et al. 2017), with the former explanation preferred in the literature (Sinclair et al. 2023; Rodríguez-Ovalle et al. 2024). During strong auroral events, driven by interactions of the solar wind with plasma emitted by Io, Jupiter has experienced global stratospheric heating due to its aurora (Migliorini et al. 2023; O’Donoghue et al. 2021). Likewise, Neptune’s aurora is not confined to its poles, and acts to significantly heat its atmosphere Melin et al. (2025). It is therefore plausible for a local aurora on SIMP-0136 to affect the global temperature structure.

From radio observations, Kao et al. (2018) inferred the presence of a strong magnetic field, of the order of 3000 G. This is far stronger than the magnetic fields that drive the aurora on Jupiter, whose magnetic field strength is about 4 G. Such strong magnetic fields will accelerate electrons, which will follow the field lines before precipitating in the stratosphere. The depth of the atmosphere at which the electrons precipitate depends on the kinetic energy of the electron beam (Pineda et al. 2024). To precipitate between 1 and 10 mbar, electron energies between 100 and 1000 keV are required. In Jupiter, the electrons typically carry 100 keV of energy (Gérard et al. 2014; Allegrini et al. 2020), so it is plausible that the higher magnetic field strength of SIMP-0136 is capable of accelerating the electrons to high enough energies to precipitate at the pressure of the observed hotspot. Given the relatively high pressure at which a 1000 keV electron beam deposits its energy, Pineda et al. (2024) also discuss how the reaction timescales of H₃+ with water are sufficient to prevent significant H₃+ emission, which would explain the lack of observed IR auroral emission in T-dwarfs.

We can estimate the required power deposited in the stratosphere to maintain the thermal inversion by comparing the difference in emission between a model with the inversion to a model with an isothermal stratosphere, shown in Fig. 14. For this comparison we use the best fit retrieved spectrum at two different phases, and compare them to identical models with the temperature gradients above 0.1 bar set to 0. The two phases plotted were chosen based on the variability map in Fig. 6 which highlights the phases where there appears to be enhanced emission in the 3.3 µm and 7.7 µm methane absorption features. The structure of the ratio between the retrieved best-fit models and those with the isothermal stratospheres clearly resembles the forward models of Fig. B.1. By taking the difference between the models, we can calculate the difference in the bolometric luminosity, and thus the emitted power. By conservation of energy, this emitted power must be equal to the power input into the stratosphere to maintain the thermal inversion. We find that for between the minimum inversion strength and the maximum, this power varies between 3.8-4.2 × 10¹⁹ W. This is significantly stronger than the IR auroral emission of Jupiter, which emits between 1-5 × 10¹³ W (Gérard et al. 2014). However, in relative terms, Jupiter’s aurora accounts for about 0.002% to 0.01% of Jupiter’s luminosity, while for SIMP-0136 this is about 0.5% of the total emitted power. By scaling Jupiter’s auroral UV emission of 10¹² W, Saur et al. (2021) finds that brown dwarfs can emit 10¹⁹ W in the UV. Their scaling relations imply that a similar increase would be observed in the IR, which is consistent with our measurement of the additional auroral IR emission.

We can compare the Saur et al. (2021) UV luminosity prediction to an estimate of the energy deposited in SIMP-0136’s atmosphere from electron precipitation. Juno/JEDI measurements find that for characteristic electron energies of 100 keV, the auroral electron beam deposits 100 mW/m² (Clark et al. 2018). The authors show that for 1000 keV electron energies, an aurora would deposit around 10⁴ mW/m², assuming similar electron densities surrounding the object. Optimistically assuming this is deposited evenly over the surface of SIMP-0136, we calculate a total power input of 5 × 10¹⁷ W, or about 1% of the observed excess emitted power. Likewise, integrating the ionisation rate (~ 10⁸ cm⁻³ s⁻¹) of Pineda et al. (2024) for a 1000 keV beam uniformly depositing its energy in a 25 km thick shell at 1 R_Jup produces an input power of 3 × 10¹⁷ W. This inconsistency suggests that either additional heating mechanisms are necessary to maintain the stratospheric inversion, or that the mechanisms of auroral energy deposition in brown dwarfs are physically distinct from those in Jupiter.

While electron precipitation is thought to be the primary mechanism for auroral heating of the upper atmospheres of brown dwarfs and giant planets (e.g. Sandel et al. 1982; Melin et al. 2011; Müller-Wodarg et al. 2012; Pineda et al. 2024), our rough estimates of the energy deposition from electron precipitation into the atmosphere of SIMP-0136 may be insufficient to fully explain the magnitude of the stratospheric inversion. Joule heating is another mechanism that may play a key role, as the strong electric currents from the auroral electron flux will heat the atmosphere. From Cassini and HST observations of Saturn, Cowley et al. (2004) found that the Joule heating rate may be several orders of magnitude larger than from the electron precipitation. In brown dwarfs, the magnetic fields are far stronger than in Saturn, and the higher temperature creates a higher ionisation fraction than in the cold atmospheres of the Solar System planets, and so the resulting electric currents will also be correspondingly larger in magnitude. Dynamical effects such as gravity wave (Young et al. 1997) have been proposed as an additional mechanism for heating the upper atmospheres of brown dwarfs (Freytag et al. 2010), and such waves have been observed to transport energy throughout the atmospheres of Jupiter and Saturn (Brown et al. 2022; Ingersoll et al. 2021).

Fig. 14 Comparison between retrieved best-fit spectrum and models without an upper atmosphere inversion to determine the contribution from the hot spot, similar to Fig. B.1. Top: comparison of each of the spectra, with the inset panel showing the differences in the temperature profile. Centre: ratio between the two retrieved spectra with the maximum (105°) and minimum (195°) inversion strengths, as well as a model without an inversion. Bottom: comparison of the best-fit model at 105° and the non-inverted model to the observed spectrum. The data are shown with 10σ uncertainties for visibility.

5.3 Dynamical regime

A key prediction of GCMs is that the dynamical state of brown dwarf atmospheres is dependent on both the temperature, which controls the radiative timescale (e.g. Showman & Kaspi 2013), and the rotation rate (e.g. Tan & Showman 2021). Changing these parameters changes the qualitative behaviour of the atmospheric structure, determining whether the presence and scale of bands and vortices (Hammond et al. 2023). The measurement of the surface gravity, radius, effective temperature and the rotation rate, allows us to compute the thermal Rossby Number, Ro_T and the Non-dimensional Radiative Time-scale, τ̂, following Equations (7) and (8) of Hammond et al. (2023), and using the same values for the variables of γ and κ. We find dimensional quantities of τ_rad to be between 10000 and 16000 s, assuming the bulk of the flux is emitted from the 1 bar level, and depending on whether retrieved surface gravity and effective temperature are taken from the retrievals or SEDkit analysis. With a rotation rate of Ω = 7.27 × 10⁻⁴ s⁻¹, this gives a non-dimensional radiative timescale of τ̂ between 15–23. We compute an equilibrium geopotential of 10⁵−10⁶ m² s⁻², which gives a thermal Rossby number of Ro_τ between 10⁻³−10⁻⁴, placing SIMP-0136 in a regime where radiation forcing dominates and jets are not expected. Instead, the atmosphere is expected to be largely homogeneous only small-scale spatial variation, although the predicted variability amplitude of ~0.1% from the shallow water models in Hammond et al. (2023) is inconsistent with the observed variability amplitude in SIMP-0136. This inconsistency motivates the need to for GCMs to model the dynamic state of the atmosphere across this parameter space, to better capture the observed variability compared to the simplified shallow water models.

5.4 Time-dependent retrieval parametrisations

In this study, we performed a series of nearly independent atmospheric retrievals on time-series spectra of SIMP-0136. However, this is only an initial step towards full three-dimensional atmospheric characterisation of this object. To achieve this, there are several underlying assumptions that limit our approach. First, we assume that the spectrum we observe at each time is independent. In reality, we are always observing a full hemisphere of the object, and so the emission from neighbouring longitudes is correlated. This is easily captured in 3D GCM models by integrating the emission over a hemisphere, but 3D radiative transfer is still too computationally expensive to be applied to atmospheric retrievals. New approaches, such as GPU-optimised radiative transfer codes (e.g. ExoJAX, gCMCRT, Kawahara et al. 2022; Lee et al. 2022), or machine-learning accelerated approaches (e.g. FlopPITy, Ardévol Martínez et al. 2024) may enable sufficiently fast 3D calculations for future retrieval frameworks. At the same time, new sampling approaches, using gradient-descent based methods with autodifferentiable models (Kawahara et al. 2022), or machine-learning posterior estimation (Vasist et al. 2023; Lueber et al. 2025) will reduce the number of model computations required to perform parameter estimation and model comparison, allowing for the use of more complex models.

Even with sufficiently fast model computation, it remains to be seen what practical approaches exist for parameterising time-and spatially-dependent atmospheric phenomena such as winds, bands, storms, and aurora. Several approaches have been proposed for hot Jupiters, where the dynamics can be described to first order by just the equatorial jet (Blecic et al. 2017; Dobbs-Dixon & Blecic 2022; MacDonald & Lewis 2022). However, for brown dwarfs the dynamics are subtler, without displaying any strong diurnal contrasts. From our independent retrievals, we see that several parameters, such as the effective temperature, and abundances for CO₂ and H₂S vary nearly sinusoidally. Such parameters θ could be modelled as $$\theta (t)={\theta }_{0}f(\varphi (t)+\delta )$$(6)

for phase φ(t) and initial phase offset δ. Both θ₀ and δ could be freely retrieved. In the simplest form, the function f could be sinusoidal in nature, although more complex functional forms could be used in general. Such an approach could be extended to higher order variations by including additional Fourier components to approximate an arbitrary periodic function, at the cost of a large number of additional free parameters. Using such parameterisations, models could be calculated at every binned timestep t, and jointly fit to the full set of time-resolved spectra, computing a single likelihood. To speed up computation, the models could be computed at relatively coarse time steps and interpolated onto the intermediate timesteps. From our parameter measurements for SIMP-0136, we observe that most parameters would not need to be variable, at least at the precision available with current observations. The choice of parametrisation of the time-variability is difficult, and further study of atmospheric variability using both independent retrievals and through efforts to approximate 3D GCMs will be necessary to determine appropriate prescriptions.

5.5 Future observations

Even with precise JWST observations, further study of SIMP-0136 is necessary to fully unravel its atmospheric structure and the processes driving the upper atmosphere heating. High resolution observations of the 3.3 µm feature could probe an extreme dynamic range in pressure, allowing for better characterisation of the temperature structure of the upper atmosphere. Likewise, moderate resolution observations using NIRSpec/G395H and the MIRI/MRS would allow for observations of the full set of methane lines and the detection of H₃+, improving our ability to constrain the thermal structure, abundance profiles, and observe tell-tale signs of auroral emission. Deep observations with HST could test to see if SIMP-0136 is a strong UV emitter, which would conclusively demonstrate the presence of an aurora. Finally, long baseline time-series observations will allow for the observation of the same longitude across time. This would enable us to better resolve the surface structure of SIMP-0136 and determine weather patterns that vary over timescales beyond the rotational period.

6 Conclusions

Broad wavelength observations of SIMP-0136 from the NIR to the radio make it a benchmark for understanding the processes driving variability near the L–T transition. Its complex atmosphere varies with time, experiencing an evolving thermal structure, as well as stratospheric heating which we attribute to aurora. Using our novel time-resolved retrieval method applied to high-precision data from JWST, we were able to place the most precise constraints to date on its thermal structure, chemistry, and cloudiness. For the first time, we have been able to quantify how the atmospheric state evolves over time.

1. Using petitRADTRANS atmospheric retrievals applied to JWST NIRSpec/PRISM and MIRI/LRS time-series observations, we were able to precisely constrain the properties of SIMP-0136. Using SEDKit, we found the mass to be 15 ± 3M_Jup and the radius to be 1.20 ± 0.05R_jup. The retrievals suffered from the ‘small-radius’ problem and we therefore treated the mass and radius as nuisance parameters, fixing them to 10.1 M_Jup and 0.93 R_Jup, respectively. This produces a surface gravity (log g = 4.49 ± 0.01) result that is consistent with the SEDKit estimate and provided a good fit to the observations. The effective temperature was found by integrating the models calculated by the retrievals, and was found to be 1245 ± 1.4 K, varying from 1243 K at the coldest to 1248 K at the warmest;

2. We found that SIMP-0136 has a stratospheric inversion, peaking at 3 × 10⁻³ bar, and with a maximum temperature approximately 300 K warmer than if the atmosphere remained isothermal above the base of the stratosphere near 0.03 bar. We hypothesise that this inversion is due to auroral heating driven by electron precipitation. SIMP-0136 must either have a very strong aurora to sufficiently heat the stratosphere or additional heating sources are required to maintain the inversion. More detailed modelling, including ion precipitation, Joule heating, and related magnetic effects is likely required to fully understand the heating mechanisms involved;

3. The bulk atmospheric metallicity is found to be constant as a function of phase. Using precisely measured chemical abundances, we found that the atmospheric metallicity is 0.18 ± 0.01.This is only slightly enriched compared to low mass stars in Carina Near association (Biazzo et al. 2012). The atmospheric C/O ratio was found to be 0.65 ± 0.01; however, when a 20% oxygen sequestration into silicate clouds is accounted for, the bulk C/O ratio is 0.54 ± 0.01, consistent with the solar value;

4. In addition to the chemical abundance results, we demonstrated that our study is sensitive to vertical chemical gradients in the photosphere. We demonstrated that an altitudevarying chemical profile is necessary to accurately fit the CH₄ absorption features. Methane and carbon monoxide are in chemical disequilibrium. The methane abundance decreases at the location of the temperature inversions, while both CH₄ and CO are quenched between ∼10 bar and the base of the inversion;

5. We performed time-resolved atmospheric retrievals to investigate changes in atmospheric properties over time. The primary driver of the variability of SIMP-0136 is the changing thermal structure as a function of phase, similar to the temperature variations predicted by Showman & Kaspi (2013). Changes in the troposphere drive an overall change in the effective temperature of about 5 K during one rotation, while changes in the strength of the stratospheric inversion are tied to variations associated with changes at wavelengths dominated by CH₄ absorption. While patchy clouds were necessary to fit the spectrum at wavelengths shorter than 2 µm, we did not observe statistically significant variations in the cloud properties or coverage fraction as a function of phase;

6. We found that most atmospheric properties are consistent with being constant across the surface of SIMP-0136. However, we found that the abundances of CO₂ and H₂S vary as a function of phase and were correlated with the effective temperature, suggesting that we may be observing chemical changes driven by dynamics and storms.

SIMP-0136 is uniquely well-suited to studying the dynamic auroral, thermal, and chemical processes that drive variations in brown dwarf atmospheres. Further observations will reveal not only how these processes shape the atmosphere, but the timescales over which they vary. Surprisingly, we find that changes in the cloud coverage are not necessary to explain the observed variability, although this is often pointed to as the primary mechanism. Studies of a larger sample objects on both sides of the L-T transition will be necessary to determine where cloud variability becomes a relevant mechanism. Lastly, model development is necessary to fully exploit the high-precision JWST data, but also to reveal the physical processes that cause the measured changes in the thermal structure and auroral heating strength.

Acknowledgements

The authors would like to thank the anonymous reviewer for their thorough and helpful feedback. We would also like to thank Elspeth Lee and Paul Mollière for their insightful discussions about this project. E.N., J.M.V., and C.O’T. acknowledge support from a Royal Society – Research Ireland University Research Fellowship (URF/1/221932, RF/ERE/221108). M.L. and M.S. acknowledge support from Trinity College Dublin via the Trinity Research Doctoral Awards. This work is based on observations made with the NASA/ESA/CSA James Webb Space Telescope. The data were obtained from the Mikulski Archive for Space Telescopes at the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-03127 for JWST. These observations are associated with program GO:3548, and can be found at 10.17909/pfnd-md36. This work was supported in part by JWST-GO-03548.004-A. Computations were performed on the HPC system Viper at the Max Planck Computing and Data Facility. B.B. acknowledges support from UK Research and Innovation Science and Technology Facilities Council [ST/X001091/1]. A.M.M. acknowledges support from the National Science Foundation Graduate Research Fellowship Program under Grant No. DGE-1840990. Approved for unlimited public release (United States Air Force, Public Affairs # USAFA-DF-2025-476). The views expressed in this article are those of the authors and do not necessarily reflect the official policy or position of the United States Air Force Academy, the Air Force, the Department of Defense, or the U.S. Government. Software used: petitRADTRANS, species, jwst-pipeline, pyMultiNest, phot_utils, Python, numpy, matplotlib, astropy.

Appendix A Benchmarking SIMP-0136

To infer the bulk properties of SIMP-0136, we used SEDkit, which reproduces the methods of Filippazzo et al. (2015). We found L_bol = −4.641 ± 0.003, T_eff = 1136K ± 22 K, Radius = 1.20 R_Jup ± 0.05R_Jup, Mass = 15 M_Jup ± 3M_Jup and log g = 4.4 ± 0.1. These measurements are largely consistent with existing literature values, (e.g. Vos et al. 2023), although the improved precision and spectral coverage of the JWST observations reduce the uncertainty on the measurements. We also used the rejection sampling method of Dupuy et al. (2023); Miles et al. (2023) to compare the properties of SIMP-0136 to the evolutionary models of Saumon & Marley (2008). These results were consistent with the SEDkit measurements, finding T_eff = 1136K ± 22 K, radius = 1.15 R_Jup ± 0.05R_Jup, Mass = 15.8 M_Jup ± 2.3M_Jup, and log g = 4.47 ± 0.1. Similarly to VHS 1256 b (Miles et al. 2023), the posterior distributions were slightly bimodal, as SIMP-0136 is both young and has a mass near the 13 M_Jup deuterium burning limit.

Fig. A.1 Fits of three self-consistent models to the SIMP-0136 spectrum at 0° phase. Statistical goodness-of-fits are poor, with Diamondback having the lowest reduced χ2 at 6327. Inferred parameters also significantly vary between each of the fits.

We performed a small set of fits using self-consistent radiative-convective equilibrium grids using species (Stolker et al. 2020) to determine a baseline for the planet properties. We used the low resolution ExoRem models (Charnay et al. 2018; Blain et al. 2021), the cloud-free ATMO models from Petrus et al. (2023), and the cloudy Sonora Diamondback models (Morley et al. 2024). The fits of each of these are included in appendix A. None of the self-consistent grids found good fits to the JWST spectrum, with reduced χ² values much greater than 1000. There was also significant variation in the inferred properties, with T_eff varying between 1000 K and 1325 K, log g between 3.5 and 5, and [M/H] between −0.2 and 0.2. While these results are broadly consistent with the literature, they lack the accuracy or goodnessof-fit to characterise the observed variability in the spectrum. This motivates the need for a retrieval analysis of SIMP-0136, whereby the precise spectra obtained with JWST can be used to infer the atmospheric parameters, which may disagree with forward model assumptions.

We additionally provided synthetic photometry for a range of filters, as well as the variability amplitude in each filter in Table A.1.

Appendix B Forward modelling of variability processes

The primary drivers of variability in SIMP-0136 have been attributed to varying cloud coverage (e.g. Artigau et al. 2009; Vos et al. 2023; McCarthy et al. 2024), auroral processes driving heating in the upper atmosphere (e.g. Pineda et al. 2017; Richey-Yowell et al. 2020), changing carbon chemistry (e.g. Tremblin et al. 2016; Tremblin et al. 2019, 2020; Biller et al. 2024; McCarthy et al. 2025). We used pRT to model how each of these processes can impact the spectrum of SIMP-0136, keeping all other parameters constant to isolate the independent effects. Across these three processes, we found that each of them imparts distinct changes in the spectrum, providing confidence that retrievals were not only be able to distinguish between the impacts, but could be used to characterise the extent to which each of these processes contributes to the observed spectroscopic variability. When compared to the typical uncertainties on each wavelength bin of 1%, we demonstrated that relatively modest changes in parameters resulted in easily detectable changes in the observed spectrum.

Fig. B.1 Variation in the spectrum due to an upper atmosphere temperature perturbation, centred at 10−4 bar. Due to the changes in the temperature structure, both the equilibrium chemical abundances and the opacities structure change, inducing changes in the outgoing emitted spectrum. The largest impact of this upper atmosphere warming is on the 3.3 µm absorption feature, due to its relatively high abundance in the upper atmosphere and large cross section at 3.3 µm, which causes the atmosphere to become optically thick at lower pressures.

B.1 Upper atmosphere heating

SIMP-0136 is known to be auroral from radio observations (Kao et al. 2016, 2018). Auroral processes can induce heating of the upper regions of the atmosphere due to the ionisation of atoms by high-energy electrons, leading to an increase in the temperature, as has been measured on Jupiter (e.g. Gérard et al. 2014; O’Donoghue et al. 2021; Sinclair et al. 2023; Rodríguez-Ovalle et al. 2024), Saturn (Gezari et al. 1989), Neptune (Melin et al. 2025) and on Y-dwarf WISE-J1935 (Faherty et al. 2024). To test the impact of the upper atmosphere thermal structure on the observed spectrum, we generated an emission spectrum using parameters derived from test retrievals of SIMP-0136. We used an equilibrium chemistry model (easyCHEM, Mollière et al. 2017), with solar metallicity and C/O (0.55, Lodders et al. (2009). This was repeated using vertically constant abundances to ensure variations in the spectrum were attributable to changes in the relative absorption and emission in each layer, as opposed to changes in the chemical abundances. As a baseline, we use an Eddington temperature profile (Eddington 1930), where $$T^4(P)=\frac{3T_{\mathrm{i}\mathrm{n}\mathrm{t}}^4}{4}(\frac{2}{3}+\frac{{\kappa }_{\mathrm{I}\mathrm{R}}P}{g}),$$(B.1)

setting T_int = 773 K, κ_IR = 0.18 and log g = 4.49 based on a test retrieval of SIMP-0136 using the Guillot profile. We applied a Gaussian perturbation to the temperature profile, centred at 10⁻⁴ bar and with a width of 1 bar. The peak of the perturbation was set to 50 K and 100 K from the nominal temperature profile. This is a similar setup to Morley et al. (2014), who identified such upper atmosphere heating as a source of variability in brown dwarfs.

The temperature profiles, together with the spectra computed for each model are shown in Fig. B.1. The increase in the upper atmosphere temperature primarily impacts the 3.3 µm methane feature. At the temperatures of SIMP-0136, this is the strongest source of opacity in the atmosphere, and thus probes the lowest pressures, and is the only molecule which probes the region affected by the temperature perturbation. With sufficient increase in temperature (around 400 K at the effective temperature of SIMP-0136), this methane feature transitions from an absorption feature to an emission feature, as observed in a Y-dwarf by Faherty et al. (2024). Given the clear absorption feature in the SIMP-0136 spectrum, we therefore expect to be sensitive to the temperature structure at pressures as low as 10⁻³ bar. This temperature sensitivity is limited by the low resolution of NIR-Spec/PRISM, and higher resolution observations would probe a larger dynamic range in pressure.

These results were similar between the equilibrium models and for the constant abundance models, which saw an even larger change in the 3.3 µm methane feature, as well as larger changes in the 2.3 µm CO feature. We also explored varying the location of the inversion, as well as how the inversion impacts objects with different effective temperatures. At higher pressures, the inversion acts similar to an increase in the effective temperature, broadly increasing the flux emitted throughout the spectrum. This results in a closer match to the flux ratios of Morley et al. (2014), who placed the hotspot at pressures of 0.1 bar or deeper. At lower pressures, the impact becomes more localised to the strongest opacity sources, primarily the 3.3 µm and 7.7 µm CH₄ features. At lower effective temperatures, the dominant source of opacity is ammonia, and so the strongest impact of an inversion is on the 10 µm absorption feature. At higher temperatures, as CO becomes dominant over CH₄ the relative impact on the methane features decreases, while the impact on the CO features increases.

B.2 Carbon chemistry

Changes in the carbon chemistry have been identified as a potential mechanism for the observed variability in the absorption features of carbon-bearing molecules, and thus precise characterisation of the abundances of these species is necessary to explore their evolution over time (Tremblin et al. 2016; Tremblin et al. 2019, 2020). Equilibrium chemistry calculations show that varying the temperature at a given pressure will change the chemical composition. While many L-T transition objects are thought to experience strong vertical mixing (Mukherjee et al. 2022), which homogenises the chemical abundances, this has not yet been tested using a data-driven approach. Furthermore, the vertical mixing strength is expected to decrease with altitude (Tan & Showman 2021; Mukherjee et al. 2022), allowing the chemical abundances to equilibrate, albeit at very low temperatures; then, the reaction rates become low enough to remain out of equilibrium. The resulting abundances in the upper atmosphere would therefore be sensitive to upper atmosphere heating.

To determine the sensitivity to changes in the abundance profile, we calculated an emission spectrum for each of the profiles of, keeping all species other than CH₄ vertically constant. For the equilibrium case, we only allow the CH₄ abundance to vary under equilibrium assumptions ([M/H]=0, C/O=0.55), again keeping all other abundances vertically constant. The abundances at the top and bottom of the atmosphere, as well as at the pressure node were chosen to approximate the equilibrium case. Model 1 is vertically constant, with an abundance equal to the equilibrium abundance in the upper regions of the photosphere (between 1-10⁻³ bar). Model 2 reaches the same abundance at 10⁻³ bar, and approximates the slope in the very upper atmosphere. Model 3 reaches the same abundance deeper in the atmosphere, near 1 bar, and a comparison with Model 4 allows us to determine the sensitivity to the chemical gradient in the photosphere. Each of these profiles is shown in the top panel of Fig. B.2.

Fig. B.2 Top: Set of mass-fraction CH4 profiles, demonstrating the linear spline parameterisations. The black dashed line shows the equilibrium profile for CH4 for [M/H]=0 and C/O = 0.55, for the temperature profile illustrated by the grey dotted line. Bottom: Equilibrium chemistry models have large variations in abundances throughout the atmosphere. In this model, we fixed all of the abundances to the median retrieved values of Table 3, except for that of methane. The different colours correspond to the same abundance profiles as in the top panel. Even though these parameterisations only affect the methane abundance, the impact on the spectrum is clearly distinct from the variations caused by the hotspot.

The resulting spectra are shown in the bottom panel of Fig. B.2. The relative strengths of the 1.6 µm and 2.2 µm CH₄ features make each case clearly identifiable relative to each other and to the equilibrium case. The structure of the variation is clearly distinct from the variations induced by changing the thermal structure in the upper atmosphere, and primarily changes the strength of the 3.3 µm CH₄ feature. We expect this sensitivity to chemical gradients to extend to other species as well. In practice, the methane abundance could not decrease without a sink for the additional carbon atoms in another carbon bearing species, and from equilibrium chemistry we would expect the CO abundance to increase in proportion to the decrease in methane.

Fig. B.3 Cloud patchiness is thought be the primary driver of near IR variability. Here we show that changes in the silicate cloud patchiness of only a few percent have large impacts on the spectrum, with larger changes at short wavelengths. This is again clearly distinct from the spectral changes due to the thermal structure or composition.

B.3 Patchy clouds

Patchy cloud coverage is thought to be the primary driver of variability in L-T transition objects, and has been identified as a primary mechanism for the variability in SIMP 0136 (McCarthy et al. 2024). As the silicate clouds condense deeper in the atmosphere, the cloud deck occurs near the photosphere, leading to relatively substantial changes in the cloud coverage fraction over the surface. We use the same equilibrium model as with the baseline for the temperature perturbation. Test retrievals indicated that SIMP-0136 is around 90% covered in silicate clouds, somewhat more than the 70% found in Vos et al. (2023). We therefore varied the cloud patchiness by 2% about this nominal value, as shown in Fig. B.3. This change in the cloud coverage induces a 20% change in the NIR flux, where the continuum absorption of the silicate clouds increases. This degree of the variation demonstrates that we are sensitive to small changes in the cloudiness of the object, of the order of a percent. There are also smaller changes induced in the K and M bands, where H₂O and CO absorption deep in the atmosphere are impacted by the presence of the silicate cloud layer. The overall shape of this variation reproduces the patchy models of Morley et al. (2014), although we find a greater impact on the spectrum for a smaller change in the cloud coverage, at a similar effective temperature. While this patchiness is clearly idealised and does not account for degeneracies within the cloud model which may compensate for decreased cloud coverage with, for example, an increased mass fraction of the cloud, it highlights that the silicate clouds impart a unique spectral variation which cannot be replicated by varying the temperature structure or chemistry.

Appendix C Additional retrieval results

Fig. C.1 Temperature gradient as a function of pressure. At pressures higher than 17 bar, the retrieved temperature gradient is steeper than an adiabatic profile, and therefore convectively unstable, while at pressures lower than 17 bar it is shallower than an adiabatic profile and is convectively stable.

Fig. C.2 Measurements of all atmospheric parameters as a function of phase for the fixed cloud retrieval, which best reproduced the observed variability. Also indicated are the median and ±1σ confidence intervals for each set of measurements.

Fig. C.3 Combined volume mixing ratios of gas-phase species from all phases for the fixed-cloud retrieval setup. Shown in dashed lines are the abundance profiles in equilibrium for [M/H] = 0.18 and C/O=0.55. The grey shading indicates pressures at which the retrievals are insensitive to the chemical abundances.

Fig. C.4 Residuals at all phases between the JWST data and the best-fit models for the fixed-cloud retrieval setup. Uncertainties on data too small to see at this scale, typically about 0.1%. The structure present across phase indicates that the primary source of the residual discrepancy is systematic in nature, due to either incomplete modelling or artefacts from the data processing.

References

All Tables

All Figures

Fig. 1 NIRSpec/PRISM (left) and MIRI/LRS (right) observations of SIMP-0136. Top: emission spectra in 15° rotation bins. Each binned NIRSpec spectrum is averaged over 90 s of integration, while each LRS spectrum is averaged over 38 s integrations to ensure that the overall signal to noise of the input spectra for the retrievals is proportional to the emitted flux. Centre: uncertainties calculated by the JWST pipeline for each spectrum in blue. Standard deviation of the 90 s (115s for the LRS) intervals in red, demonstrating that the pipeline uncertainties accurately reflect the underlying statistical distribution of the noise. Bottom: binned variability maps, also known as dynamic spectra. These variability maps are aligned in phase, as the MIRI observations were taken after the NIRSpec observations and are binned to highlight the wavelength dependence of the variability. The 24 phase bins used in these maps correspond to the binned spectra used as inputs for the retrievals.

In the text

Fig. 2 Best-fit model from the fiducial retrieval to the observed spectrum at 0° phase, whose uncertainties are smaller than the displayed line width. This model assumed a fixed mass and radius of 10.381 MJup and 0.9329 RJup respectively. Also shown for visual reference are the log opacities at the 1 bar level of the primary absorbers in each wavelength range. The CH4 abundance was allowed to vary as a function of altitude. The reduced χ2 of the t s T. However, the model clearly matches the key absorption features of H2O, CH4, CO, and CO2, as well as additional species such as H2S and NH3.

In the text

	Fig. 3 Emission contribution function for SIMP-0136 at 0° phase. The emission contribution is shown in square-root colour scale to highlight subtle variations. Significant contributions are present across four decades in pressure, bounded by a cloud layer near 10 bar and by methane absorption up to 10−3 bar.
In the text

Fig. 4 Retrieved temperature profiles for each phase bin. Temperature profiles are given as 90% confidence intervals. In grey is the emission contribution function, averaged over wavelength. We show the square root of he contribution to highlight that there is small, but significant emission from the upper atmosphere, coincident wit the location of the temperature inversion. In black lines are condensation curves for enstatite, forsterite, and iron, showing that the retrieved cloud base pressure for the silicate cloud is coincident with the expected condensation location W also show the temperature profile from self-consistent ExoRem model (log g = 4.5, [M∕H]= 0, C/O= 0.55, Teff = 1000 K) as a red dotted line.

In the text

	Fig. 5 Volume mixing ratio profiles of gas-phase species for all phases. The shaded regions indicate 68% confidence intervals, taken from the set of median abundances measured at each phase. The measured precision for a single phase is typically an order of magnitude smaller than the variation between phases.
In the text

	Fig. 6 Top: phase-resolved best-fit spectra divided by the mean retrieved spectrum for NIRSpec/PRISM (left) and the MIRI/LRS (right), taken from the fiducial retrieval. Bottom: difference between the observed and retrieved variability maps. The remaining difference is normally distributed, centred at 0, and with a standard deviation of 0.3% for NIRSpec and 1% for MIRI.
In the text

	Fig. 7 Same as Fig. 6, but with sets of parameters kept constant across phase. Top: variability map keeping the temperature profile constant, showing that the retrievals are unable to reproduce the observed variability. Bottom: variability map keeping the cloud parameters constant, which better reproduces the observed variability than the full retrieval.
In the text

	Fig. 8 Fig. 8. Measurements of key atmospheric parameters as a function of phase, taken from the fiducial retrieval. Also indicated are the median and �1σ confidence intervals for each set of measurements.
In the text

	Fig. 9 Top: comparison of normalised effective temperature and the ∇T4 parameter with the clustered light curves of M25. The effective temperature is strongly correlated with the light curve, while ∇T4 is anticorrelated. Bottom: matrix of Pearson correlation coefficients between atmospheric parameters and clustered light curves from M25. Nonsignificant correlations (|r| < 0.271) are not shown.
In the text

Fig. 10 Pearson correlation coefficients between the phase-variation of each pair of retrieved parameters in the fiducial retrieval. Using a two-sided t-test with 24 samples, the minimum statistically significant correlation at a 90% confidence level is 0.271, thus setting the limits on which correlations are shown. While some correlations are natural degeneracies of the model, such as the correlation between ∆ logPCH4 and the CH4 abundance at that pressure node, other pairs are likely astrophysical in nature.

In the text

Fig. 11 Correlation between the phase-variation of the temperature profile and the clustered light curves from M25. Top: example temperaturevariation curves drawn from three representative pressure levels in the atmosphere (highlighted boxes) compared to correlated or anticorrelated light curve clusters. Bottom: correlation coefficient between the temperature-variation at that pressure and the clustered light curves for each pressure level, as well as the correlation to the effective temperature in the right-most column. The pressure levels for the temperature curves in the top panels are highlighted with coloured boxes. Horizontal lines and annotations indicate the location of atmospheric features.

In the text

	Fig. 12 Fig. 12. Effective temperature for each observed phase. Over one rotation Teff varies by 5 K.
In the text

	Fig. 13 Elemental abundance ratios derived from retrieved molecular abundances combined from all phases.
In the text

Fig. 14 Comparison between retrieved best-fit spectrum and models without an upper atmosphere inversion to determine the contribution from the hot spot, similar to Fig. B.1. Top: comparison of each of the spectra, with the inset panel showing the differences in the temperature profile. Centre: ratio between the two retrieved spectra with the maximum (105°) and minimum (195°) inversion strengths, as well as a model without an inversion. Bottom: comparison of the best-fit model at 105° and the non-inverted model to the observed spectrum. The data are shown with 10σ uncertainties for visibility.

In the text

	Fig. A.1 Fits of three self-consistent models to the SIMP-0136 spectrum at 0° phase. Statistical goodness-of-fits are poor, with Diamondback having the lowest reduced χ2 at 6327. Inferred parameters also significantly vary between each of the fits.
In the text

Fig. B.1 Variation in the spectrum due to an upper atmosphere temperature perturbation, centred at 10−4 bar. Due to the changes in the temperature structure, both the equilibrium chemical abundances and the opacities structure change, inducing changes in the outgoing emitted spectrum. The largest impact of this upper atmosphere warming is on the 3.3 µm absorption feature, due to its relatively high abundance in the upper atmosphere and large cross section at 3.3 µm, which causes the atmosphere to become optically thick at lower pressures.

In the text

Fig. B.2 Top: Set of mass-fraction CH4 profiles, demonstrating the linear spline parameterisations. The black dashed line shows the equilibrium profile for CH4 for [M/H]=0 and C/O = 0.55, for the temperature profile illustrated by the grey dotted line. Bottom: Equilibrium chemistry models have large variations in abundances throughout the atmosphere. In this model, we fixed all of the abundances to the median retrieved values of Table 3, except for that of methane. The different colours correspond to the same abundance profiles as in the top panel. Even though these parameterisations only affect the methane abundance, the impact on the spectrum is clearly distinct from the variations caused by the hotspot.

In the text

	Fig. B.3 Cloud patchiness is thought be the primary driver of near IR variability. Here we show that changes in the silicate cloud patchiness of only a few percent have large impacts on the spectrum, with larger changes at short wavelengths. This is again clearly distinct from the spectral changes due to the thermal structure or composition.
In the text

	Fig. C.1 Temperature gradient as a function of pressure. At pressures higher than 17 bar, the retrieved temperature gradient is steeper than an adiabatic profile, and therefore convectively unstable, while at pressures lower than 17 bar it is shallower than an adiabatic profile and is convectively stable.
In the text

	Fig. C.2 Measurements of all atmospheric parameters as a function of phase for the fixed cloud retrieval, which best reproduced the observed variability. Also indicated are the median and ±1σ confidence intervals for each set of measurements.
In the text

	Fig. C.3 Combined volume mixing ratios of gas-phase species from all phases for the fixed-cloud retrieval setup. Shown in dashed lines are the abundance profiles in equilibrium for [M/H] = 0.18 and C/O=0.55. The grey shading indicates pressures at which the retrievals are insensitive to the chemical abundances.
In the text

	Fig. C.4 Residuals at all phases between the JWST data and the best-fit models for the fixed-cloud retrieval setup. Uncertainties on data too small to see at this scale, typically about 0.1%. The structure present across phase indicates that the primary source of the residual discrepancy is systematic in nature, due to either incomplete modelling or artefacts from the data processing.
In the text