© The Authors 2026

1 Introduction

Brown dwarfs are substellar objects that cannot sustain stable fusion and steadily cool over time (Burrows et al. 1997; Baraffe et al. 2003; Saumon & Marley 2008). Y dwarfs represent the coolest end of the brown dwarf sequence, with currently known effective temperatures of 250 K ≤ T_eff ≤ 500 K (Leggett et al. 2021) and rapid rotation periods of a few hours (Tannock et al. 2021), occupying a physical regime that overlaps with the atmospheric conditions of giant planets (Showman et al. 2020). More than ten years after Wide-field Infrared Survey Explorer (WISE) discovered the first Y dwarfs (Cushing et al. 2011; Kirkpatrick et al. 2011; Luhman 2014), the James Webb Space Telescope (JWST) has both expanded that sample and is poised to detect and characterize young sub-Jovian companions. Consequently, an expanding sample of these cool brown dwarfs has become a primary focus of recent studies (e.g., Leggett & Tremblin 2023; Beiler et al. 2024; Lew et al. 2024). Their atmospheric conditions make them valuable laboratories for testing dynamical and cloud physics, which could deliver insights into colder planetary atmospheres.

The spectra of Y dwarfs indicate the presence of water (Burrows et al. 2003; Morley et al. 2014b), sulfide, and salt clouds (Morley et al. 2012, 2014a), and point to vigorous vertical mixing (Robinson & Marley 2014), signaling that their atmospheres depart from simple equilibrium expectations (Lacy & Burrows 2023). Using the Ackerman & Marley (2001) model, Morley et al. (2012) suggested that, compared to the silicate clouds of L dwarfs, cooler brown dwarf atmospheres are expected to contain different condensates (such as Cr, Na₂S, MnS, and KCl) and ran a suite of models taking these species into account, finding that cloudy models explain the spectra of these objects better than cloud-free models. Morley et al. (2014a) highlight that, although photometric variability has been a focus of L/T transition studies, late-T to early Y dwarfs (T_eff ≥ 375 K) also display a degree of variability (Cushing et al. 2016). They suggest that this variability could be a result of patchy cloud cover or the global circulation creating an irregular temperature structure, and try to mimic these effects by combining multiple 1D models. They conclude that both of these effects induce variability in different spectral ranges.

Analytical studies have provided a first guide to the relevant dynamics for these substellar atmospheres and a growing body of theory links brown dwarf variability to interacting thermal, dynamical, cloud, and chemical processes. Showman & Kaspi (2013) presented a simple analytical framework, estimating the horizontal temperature variations, together with zonal-and vertical wind magnitudes, connecting them to a dimensionless efficiency parameter η, quantifying how much of the present heat flux is used to drive the circulation in these atmospheres. Using a two-layer shallow-water model that is forced from the interior by randomly varying perturbations, Zhang & Showman (2014) demonstrated that brown dwarf atmospheres can settle into either of two broad circulation regimes: large-scale zonal jets or fields of local vortices. Their outcomes depended on the balance between convective forcing and radiative damping. Hammond et al. (2023) used a similar one-layer shallow-water model and by exploring a large parameter space covering possible brown dwarf configurations, identified four dynamical regimes. Both works conclude that, given the radiative forcing expected of Y dwarfs, they should lie in a regime dominated by zonal jets.

In addition to these analytical works, recent years saw multiple works that ran 3D general circulation models (GCMs) of field brown dwarf atmospheres. The first such work by Showman & Kaspi (2013) ran an idealized brown dwarf case with the MIT-gcm (Adcroft et al. 2004) and found brown dwarf atmospheres to be rotation-dominated. By comparing a fast and a slow rota-tor case, along with analytical considerations, they concluded that, on large scales, the flow should align in columns parallel to the rotation axis and have latitudinal variations resulting in pole-to-equator temperature differences. Showman et al. (2019) performed idealized GCM simulations, forcing the atmosphere through a physically motivated parametric model of interior thermal perturbations, and they found that their models produced zonally-banded structures ubiquitously, along with oscillations in the predicted structures that they tie to the observed variability.

Tan & Showman (2021a) explored the effects of cloud thermal feedback and its potential to drive atmospheric circulation in brown dwarfs using a model assuming a Cartesian geometry. They showed that cloud and rotation coupling can produce significant temperature contrasts, strong winds, and persistent, time-variable structures under suitable conditions. By coupling the MITgcm to a cloud model, Tan & Showman (2021b) conducted a high-resolution study of a multitude of brown dwarf scenarios of varying rotation rates. They confirmed that cloud radiative feedback can be the main driver of dynamical activity in brown dwarf atmospheres and explored possible reasons behind the observed light curve variability (Cushing et al. 2016; Artigau 2018). They concluded that variability may arise both from wave activity triggered by cloud feedback and from changes in cloud coverage with viewing geometry. Tan (2022) extended the model to explore jet formation and tracer transport in brown dwarf and isolated giant planet atmospheres. The study found that thermal forcing can generate strong zonal jets, vertical shear, and long-term oscillations in equatorial winds, depending on the strength of radiative and frictional damping. It also showed that these circulations can mix cloud particles and chemical tracers vertically, with the extent of mixing varying by particle size and chemical properties.

In a series of works, Lee et al. (2023) and Lee et al. (2024), using ExoFMS (Lin 2004; Lee et al. 2020), presented GCM simulations of brown dwarf scenarios coupling their radiative transfer, chemistry and clouds with increasingly more sophisticated treatments of these physical processes. They found that these processes feed back onto each other in a complex fashion and are crucial in explaining the observed dynamical and thermal structure, along with the variability of a given brown dwarf atmosphere.

In a related work, Lee & Ohno (2025) developed a two-moment cloud microphysics module, coupled it to the ExoFMS GCM and applied it to an isolated Y dwarf case. They found that the parameter regime they tested produced weak dynamical activity and emphasized the need for a wider sweep of parameter space, an issue we address with the current work. This work extends GCM studies of cold brown dwarf atmospheres with a targeted survey of Y dwarf parameter space using a unified physical treatment. We span rotation rates and effective temperatures characteristic of early Y dwarfs and exclude the cold tail where H₂O clouds form, while variations in R_p, log g, metallicity, and disequilibrium chemistry are left for future studies.

The structure of the paper is organized as follows. In Section 2.1, we introduce diagnostic quantities that have been used in previous theoretical works on brown dwarfs and give an analytical order-of-magnitude estimate of the dynamical regime we expect Y dwarfs to fall into. Section 2.2 gives an account of our chosen sample of models, along with our reasoning for the presented modeling choices. Sections 2.3–2.4.6 detail our GCM setup along with the additional physics modules that were developed for this study. Section 3 goes over our main results, detailing identified trends in our sample of 12 Y dwarf cases, their corresponding flow regimes and the influence of mineral clouds in their atmospheres. Subsequently, in Section 4, we discuss how our results fit into the existing literature, offer insights from the first large-scale GCM study on Y dwarfs, and lastly, highlight possible improvements for future studies.

2 Methodology

2.1 Characteristic flow quantities

To ground our global circulation modeling in first-principles expectations and as a prelude to the numerical simulations, we derive order-of-magnitude estimates for the thermal and dynamical regimes utilizing prior shallow water theoretical work done on brown dwarf atmospheres. The analytical framework of Hammond et al. (2023) uses a dimensionless two-parameter space to describe the circulation regimes of a brown dwarf atmosphere. The first characteristic quantity is the thermal Rossby number $$R{o}_{\text{T}}=\frac{{\Phi }_0}{{\left(2\Omega R_{\text{p}}\right)}^2},$$(1)

where Φ₀ is the equilibrium geopotential. This can be connected to the standard Rossby number (Holton 1992; Showman & Guillot 2002), defined as $$Ro=\frac{U}{fL},$$(2)

through thermal wind balance Φ₀ = 2ΩUR_p (Mitchell & Vallis 2010). Here, U [m s⁻¹] is the characteristic wind speed, L [m] the characteristic length scale, and f = 2Ω sin θ [s⁻¹] the Corio-lis parameter at latitude θ for a planet rotating at rate Ω [rad s⁻¹]. While, for the purposes of this section, we can estimate Ro_T via an idealized thermal wind calculation, within the GCM we calculate the geopotential as (Zhang & Showman 2014) $${\Phi }_0=\gamma \frac{R_d^2}{c_p}T,$$(3)

where R_d [J K⁻¹ kg⁻¹] is the specific gas constant, c_p [J K⁻¹ kg⁻¹] the specific heat at constant pressure, T [K] the temperature and γ a measure of adiabacity given by $$\gamma =\left[1+\left(\frac{c_{\text{p}}}{g}\frac{dT}{dz}\right)\right].$$(4)

Both definitions of the Rossby number gauge the relative importance of inertial to Coriolis forces, where small (Ro≪ 1) values imply rotation-dominated dynamics. The second coordinate is the dimensionless radiative timescale $${\widehat{\tau }}_{\text{rad}}=2\Omega {\tau }_{\text{rad}},$$(5)

which incorporates the physical radiative relaxation time (Showman & Guillot 2002) $${\tau }_{rad}~\frac{p}{g}\frac{c_p}{{\sigma }_{\text{SB}}{T}^3},$$(6)

where p [Pa] is pressure, g [m s⁻²] surface gravity and σ_SB [W m⁻² K⁴] the Stefan–Boltzmann constant. Although different physical conditions can lead to the same pair of parameter values $(R{o}_{\text{T}},{\widehat{t}}_{\text{rad}})$, Hammond et al. (2023) state that every point in this plane corresponds to a unique qualitative circulation state, providing a useful description of the underlying dynamical state of the atmosphere.

Unfortunately, to estimate the Rossby number and dimensionless radiative timescale requires specifying the values of quantities that are the outcomes of GCM simulations. We thus use analytical scaling relations to make educated guesses for the values of these quantities. To estimate the horizontal wind velocity U we adopt the wave-driven scaling relation of Showman & Kaspi (2013). By combining geostrophic thermal-wind balance with a simple radiative-damping closure, they derived an expression of the characteristic horizontal wind speed U in a fast-rotating, nonirradiated atmosphere as $$U\simeq \sqrt{\frac{\eta c_{p}T}{4}}\frac{\ell \Delta \tilde{z}HN}{f},$$(7)

where H [m] is the pressure scale height, defined as $$H=\frac{k_{B}T}{mg},$$(8)

and N [s⁻¹] is the Brunt–Väisälä frequency $$N=\sqrt{\frac{g}{T}\left(\frac{dT}{dz}+\frac{g}{c_p}\right)},$$(9)

ℓ [m⁻¹] the characteristic wavenumber, $\Delta \tilde{z}[H]$ the vertical thickness of the overturning cell in scale heights, and η the fraction of the interior heat flux that goes into driving atmospheric waves, called the dimensionless wave driving efficiency.

Although the Y dwarf parameter range spans 250 K ≤ T_eff ≤ 500 K, we assume 400 K ≤ T_eff ≤ 600 K for our calculations presented in this section, matching our modeled cases detailed in Section 2.2. The radius is taken to be equal to the radius of Jupiter R_p = R_J, since the radii of brown dwarfs are not expected to change significantly past their initial evolution and typically fall very close within one Jupiter radius (Burrows et al. 1997; Leggett et al. 2017). For the surface gravity, we choose a representative value of log g = 4.5 and later discuss the qualitative impact of changing this parameter. Rotation periods in the range P_rot = 2.5–20 h encompass most of the current observations (Leggett et al. 2017; Tannock et al. 2021) and imply angular velocities Ω = 9 × 10⁻⁵–7 × 10⁻⁴ rad s⁻¹. We evaluate the Corio-lis parameter at a representative mid–latitude θ = 45^◦. Assuming an approximate solar–composition mixture, with volume mixing ratios (VMRs) ${\chi }_{{\text{H}}_2}=0.84$ and ${\chi }_{\text{He}}=0.16$, yields a mean molecular weight µ = 2.33, provides the thermodynamic background, with a specific gas constant $R_d=\frac{k_{\text{B}}}{\mu m_u}=3569\text{J}{\text{K}}^{-1}{\text{kg}}^{-1}$ and specific heat at constant pressure $c_p=\frac{R_d}{{\kappa }_{\text{ad}}}=12039\text{J}{\text{K}}^{-1}{\text{kg}}^{-1}$, where m_u is the atomic mass unit and κ_ad is the adiabatic gradient.

Following Showman & Kaspi (2013) and aiming for simplicity, we assume an isothermal atmosphere. This allows us to calculate $NH=R_d\sqrt{\frac{T}{c_p}}$. For wave driving, we adopt the characteristic horizontal wavenumber of the same work ℓ = 3 × 10⁻⁷ m⁻¹, which corresponds to a characteristic length scale of ≈10⁴ km, matching observational values (Simon-Miller et al. 2006) for Jupiter. A quick comparison to the Rossby deformation radius, as it is defined in Showman & Kaspi (2013) $$L_D=\frac{\Delta \tilde{z}HN}{f},$$(10)

yields a similar average value of ≈5 × 10³ km when calculated over all the ranges of values given above. The vertical extent of the overturning cell is assumed to range between $\Delta \tilde{z}=1-3$ (Showman & Kaspi 2013), meaning that the circulation spans one to three scale heights. Following the same study, we define the wave-driving efficiency as the fraction of the interior heat flux that goes into accelerating the mean flow $$\eta \equiv \frac{\mathcal{A}\Delta up}{F_{\text{OLR}}g},$$(11)

where 𝒜 is the eddy acceleration, ∆u a characteristic wind speed and F_OLR the emitted flux. Using observationally inferred 𝒜, ∆u, p, F_OLR for Earth and Jupiter, Showman & Kaspi (2013) infer η ∼ 10⁻³ for both cases and suggest a range of η 10⁻⁴–10⁻² for brown dwarf studies. In our calculations, we explore values over η = 10⁻⁴ to 10⁻¹, a range that encompasses very weak to vigorous wave forcing.

We consider two photospheric pressure values p = [0.1, 1] bar in the radiative timescale expression. Combining these values with the rotation periods above gives a dimensionless radiative time ${\widehat{\tau }}_{\text{rad\hspace{0.17em}}}$ that lies between 5 and 10³. The resulting pair $(R{o}_{\text{T}},{\widehat{\tau }}_{\text{rad\hspace{0.17em}}})$ therefore situates Y–dwarf atmospheres in the circulation regime characterized by radiative forcing with no jets expected to form (Hammond et al. 2023). The full range of possible values is illustrated in Figure 1, with the extent of the axes corresponding to the ones presented in Hammond et al. (2023) to enable a direct comparison.

Adjusting each parameter causes a corresponding shift in the position on the plot. Lower effective gravities of log(g) = 3.5–4 shift the entire plot to the right and higher effective gravities shift it to the left. Although surface gravity affects both the predicted thermal Rossby numbers and radiative timescales, in our idealized framework, its effect is more dominant on the radiative timescales. Selecting a deeper photosphere, for example, p = 1 bar (also illustrated with fainter points in Figure 1), has the same effect of shifting the points to the right and different wave-driving efficiencies shift the plot upward or downward, as expected.

Fig. 1 Summary of our calculations, presented in Section 2.1. Highlighted are the cases with gravity log g = 4.5, a wave-driving efficiency of η = 10−3 and a vertical scale of $\Delta \tilde{z}=2$ at the approximate photosphere level pphotosphere = 0.1 bar. Other possible combinations are overplotted in fainter colors to illustrate the complete space spanned by the range of values considered in this study. Dotted lines mark the edges of minimum and maximum values. The red star corresponds to (Teff = 400 K; Prot = 10 h; log g = 4.46) and is from a recent GCM study, included as a reference point.

2.2 Design of simulation cases

To capture the center of the temperature distribution of the early Y dwarf regime, while avoiding the parameter space where H₂O clouds would form, we adopt 400 K ≤ T_eff ≤ 600 K, where the warm end (T_eff = 600 K) overlaps with The late-T regime. We chose our representative Y dwarf scenarios to have a radius of R_jup, together with an effective surface gravity of log g = 4.5. Although observed Y dwarf radii vary about 25%, they are mostly around one Jupiter radius (Burrows et al. 1997; Leggett et al. 2017). For the surface gravity, while retrieval (e.g., Zalesky et al. 2019) and modeling (e.g., Lacy & Burrows 2023) studies suggest a range of 3.5 ≤ log(g) ≤ 5, we select a single representative value of log g = 4.5. As discussed in Section 2.1, changing g mainly produces a horizontal shift in ${\widehat{\tau }}_{\text{rad\hspace{0.17em}}}$ on Figure 1, whereas the shallow-water regimes identified by Hammond et al. (2023) vary most strongly in the vertical (Ro_T) direction. We, therefore, keep the surface gravity fixed and concentrate our computational efforts on exploring a wider range of effective temperatures, rotation rates and cloud species. To decide which cloud species to include in our simulations, we look at previous studies (Morley et al. 2012; Gao & Benneke 2018) and remark that for the temperature range we find ourselves in, the main relevant cloud species are KCl, MnS and Na₂S. Figure 2 illustrates the condensation curves for these species overplotted with our initial T–p profiles. Finally, we assume solar metallicity for all of our runs. Current observations put most Y dwarfs close to solar metallicities (Leggett et al. 2017) and since most of these objects are observed nearby our Solar neighborhood, this provides a reasonable assumption. Our constructed suite of run cases that are briefly summarized in Table 1.

We start our GCM runs from rest with a gray analytical T–p profile (Milne 1921; Eddington 1926) derived for a self-luminous object in the gray opacity limit $$T^4=\frac{3T_{\text{int}}^4}{4}\left[\frac{2}{3}+{\tau }_{\text{lw}}\right],$$(12)

with T describing the temperature structure with optical depth, T_int the bottom boundary internal temperature of the atmosphere, and τ_lw the long-wave optical depth. The optical depth term τ_lw(p) controls the location of the photosphere and is solely used to construct the initial T–p profile. Following the radiative equilibrium solutions calculated in previous studies for brown dwarf atmospheres (e.g., Burrows et al. 2006), we aim for an adiabatic interior, an isothermal upper atmosphere and a profile with the approximate photosphere (p ≈ 0.1–1 bar) corresponding to the effective temperature for each case. Therefore, we choose a long-wave gray opacity of κ_lw = 2.5 × 10⁻³ m² kg⁻¹ to generate our initial T–p profile using the relation (Deitrick et al. 2020) $${\tau }_{\text{lw}}={\tau }_0\frac{p}{p_{\text{ref}}}+\left(\frac{1}{f_{\text{lw}}}-1\right){\tau }_0{\left(\frac{p}{p_{\text{ref}}}\right)}^2$$(13)

where f_lw is the strength of the uniformly mixed absorbers set to 0.5, $0.5,{\tau }_0={\kappa }_{\text{lw}}\frac{p_{\text{ref\hspace{0.17em}}}}{g}$ the reference optical depth, p_ref the reference pressure of the bottom boundary set to 100 bar for our current set of simulations. It should be noted that the quadratic dependence τ_lw ∝ p² makes this solution convectively unstable at initialization. To remove this initial instability without altering the target profile, we apply a brief initialization (10 computational steps of 100 s) with dry convective adjustment and then continue with our mixing-length theory (see Section 2.4.3) routine. The T–p profile after this step is essentially unchanged and the simulation proceeds stably. Figure 2 illustrates various examples for initial T–p profiles, and we can see that the κ_lw = 2.5 × 10⁻³ m² kg⁻¹ case fits the characteristics we were aiming for.

In our tests, we also initialized simulations using an isothermal temperature profile, T (p) = T_int, letting the dynamics dictate the evolution of the T–p profile. This approach ultimately converged to the same steady-state T–p structure as with the gray profile given by Eq. (12). However, the latter profile notably accelerates convergence to the steady state, as it closely approximates the final equilibrium solution. This approach also circumvents possible numerical issues where starting from an isothermal profile can generate a confined cloud condensate layer at the bottom of the atmosphere that results in a physically unrealistic scenario.

Recent studies by Tan & Showman (2021b), Lee et al. (2023) and Lee et al. (2024) suggest that these atmospheres reach a statistical steady state after a couple of hundred days. Our simulations confirm this fact and we choose a uniform simulation time of 1000 Earth days for all of our presented results, aiming for a reasonable yet computationally feasible value. The last critical consideration for the general structure of runs is the average horizontal resolution. Since Y dwarfs are rapid rotators (Leggett et al. 2017; Tannock et al. 2021), we require a high enough resolution to capture the length scale of the dominant dynamical features. To estimate the needed length scale, we take the minimum possible Rossby deformation radius L_D for our setup, calculated for a range of plausible parameters (as they are presented in Section 2.1). The minimum, achieved only for the fastest rotators (P_rot = 2.5 h), is L_D ≈ 6.5 × 10⁵ m, corresponding to a refinement level g_level ≈ 7 as described in Mendonça et al. (2016); Deitrick et al. (2022). At the opposite end, the slowest rotators (P_rot = 20 h) have a larger minimum L_D ≈ 5.271 × 10⁶ m, which would only require g_level ≈ 4. To conservatively resolve small-scale structure and to ease cross-comparison, while also aiming to remain computationally feasible, we adopt g_level = 7 for all but the slowest cases. For the P_rot = 20 h runs, we use g_level = 6, which still comfortably resolves the deformation scale at those rotation rates. These considerations land us at a chosen average angular resolution of $\overline{\theta }\approx 1^{°}$ for the slowest rotators and $\overline{\theta }\approx {0.5}^{°}$ for the rest of the runs and are comparable to the high-resolution simulations performed by Tan & Showman (2021a) and Tan (2022).

Fig. 2 Various T–p profiles overplotted to illustrate the initial conditions from which the simulations were spun up. The solid lines correspond to the κlw = 2.5 × 10−3 m2 kg−1 curves with the lighter shades of each color illustration curves corresponding to the same internal temperature Tint with varying long-wave opacities of κlw = [2.5 × 10−5, 2.5 × 10−4, 2.5 × 10−2] m2 kg−1. Also overplotted are condensation curves for the cloud condensate species KCl, MnS, Na2S. The shaded horizontal region (p = 0.1 – 1 bar) corresponds to the approximate photosphere.

2.3 The dynamical core

In our study, we use the THOR¹ dynamical core (Mendonça et al. 2016; Mendonça et al. 2018a,b; Deitrick et al. 2020, 2022; Noti et al. 2023) to simulate the atmospheric circulation of Y dwarf atmospheres. THOR is an open-source, fully three-dimensional dynamical core that uses a horizontally explicit, vertically implicit (Tomita et al. 2002; Wicker & Skamarock 2002) time integrator on an icosahedral grid (Tomita et al. 2001) to solve the non-hydrostatic deep Euler equations. This split-time integration scheme allows THOR to maintain numerical stability with vertically propagating waves in its solutions while still taking a reasonable time step, which would otherwise be constrained by the Courant–Friedrichs–Lewy condition. Our selected grid structure avoids the usual limitations of latitude–longitude grids (Staniforth & Thuburn 2012), which experience singularities and resolution difficulties at the poles. Such problems arise in all nonuniform grid structures characterized by fixed points and substantial variations in cell sizes or resolution crowding. The icosahedral grid proposed by Tomita et al. (2001) and Tomita & Satoh (2004) provides a quasi-uniform horizontal grid layout. In the vertical direction, the THOR model employs an altitude grid and is implemented in two possible configurations. The first configuration, described in Mendonça et al. (2016) and Deitrick et al. (2020), uses a uniformly spaced altitude grid. The bottom-layer altitude is obtained from the ideal gas law, and the user defines the top of the atmosphere as a free parameter. An alternative introduced in Noti et al. (2023) provides more detailed control over the vertical spacing, enabling a nonuniform vertical structure. Although the default THOR model primarily solves the non-hydrostatic deep Euler equations, as outlined in Deitrick et al. (2020), Deitrick et al. (2022), and Noti et al. (2023), it also provides the option to implement various simplifications to the governing equations. These alternatives include quasi-hydrostatic deep or hydrostatic shallow formulations, each with its own assumptions. In this work, we use the non-hydrostatic deep equations with a uniformly spaced vertical grid structure.

To dissipate energy that accumulates at the grid scale, every run employs two explicit schemes, horizontal and vertical hyperdiffusion, and a divergence-damping operator, together with a Rayleigh sponge near the model top. The hyperdiffusion terms remove small-scale enstrophy, the divergence operator acts directly on the flow’s divergence field, and the sponge relaxes all wind components toward their zonal means above a prescribed height, thereby absorbing waves that would otherwise reflect off the top of the atmosphere. No basal drag is applied at the lower boundary. The numerical orders, coefficients, and sponge parameters match those used in our works using the THOR GCM (e.g., Akın et al. 2025) and are listed in Table 2.

2.4 Additional physics modules for brown dwarfs

As introduced in Mendonça et al. (2016) and Deitrick et al. (2020), the THOR GCM was built from the ground up for exo-planetary atmospheres and has, at the time of writing, been mainly applied to model hot Jupiters (e.g., Mendonça et al. 2018a; Deitrick et al. 2022; Noti et al. 2023; Akın et al. 2025). Brown dwarfs, meanwhile, occupy an intermediate region between giant planets and low-mass stars. Although there are significant similarities between brown dwarfs and giant planets, such as similar surface gravities, molecular compositions, and cloud chemistry (Showman et al. 2020), there are also distinct differences that require the development of additional physics modules to accurately model their atmospheres.

The main differences between the circulation regimes of brown dwarfs and hot Jupiters arise from the fact that hot Jupiter atmospheres are externally forced by their host star. Together with tidal locking and a lack of significant internal heating, the dynamics of these atmospheres can be character-ized by a pronounced horizontal (day–night) heating gradient and the interplay of this gradient with planetary rotation. In field brown dwarfs, the situation is inverted. These objects lack any external forcing and slowly radiate away the residual heat from their formation, and in more massive cases, the negligible output of halted nuclear fusion. Their atmospheres are therefore driven primarily by the interior. The resulting thermal forcing is vertical, extending from the deep convective interior up to the observable photosphere, and, when combined with the typically rapid rotation of field brown dwarfs, sets up a circulation in which convective transport and vertical mixing play a central part. Accordingly, we implement a routine that introduces thermal perturbations at the approximate radiative-convective (RCB) boundary and we replace the present convective adjustment scheme (Deitrick et al. 2020) with a mixing-length-theory formulation that treats sub-grid convection more realistically. Many widely used exoplanetary GCMs implement a convective adjustment scheme (e.g., Heng et al. 2011; Leconte et al. 2013; Deitrick et al. 2020) that mixes the entropy in a convectively unstable zone, approximating convection as an instantaneous process compared to the dynamics. While Lee et al. (2023) model a simplified brown dwarf case using convective adjustment, they also subsequently advocate for the improvements offered by their mixing-length theory implementation in their follow-up work (Lee et al. 2024). Fully convection-resolving global models (e.g., Lefèvre et al. 2022) remain prohibitively expensive for the parameter space explored here.

Throughout this study, we confine ourselves to field brown dwarfs and therefore neglect any external instellation. This choice has a direct bearing on the radiative transfer strategy. THOR presently contains three schemes: the original two-stream double-gray module (Mendonça et al. 2016; Mendonça et al. 2018b; Deitrick et al. 2020); a two-stream, picket-fence non-gray module modeled on the work of Lee et al. (2021) and described by Noti et al. (2023); and a multiwavelength “improved two-stream” method (dubbed THOR+HELIOS in Deitrick et al. 2022), that couples the THOR dynamical core to the HELIOS radiative transfer solver and provides a rigorous treatment of scattering by medium-sized and large aerosols (Heng & Kitzmann 2017; Heng et al. 2018; Malik et al. 2017, 2019). For the cloud-forming, nonirradiated objects considered here, the existing options are not ideal for our goals. The picket-fence module devotes bands to stellar shortwave and offers little benefit without instellation, while the multiwavelength THOR+HELIOS coupling is too costly for long integrations with time-varying, cell-wise cloud properties. We require an efficient radiative transfer scheme for the planetary emission that supplies a locally evaluated opacity in every grid cell so that evolving cloud opacities can be coupled to the thermal structure calculated by the solver. We therefore adopt a single thermal band with locally evaluated Rosseland-mean gas and cloud opacities, mirroring the approach of many previous brown dwarf studies (e.g., Tan & Showman 2021b; Tan 2022; Lee et al. 2023).

In summary, we have developed four complementary additions to THOR: an interior thermal perturbation model, a mixing-length sub-grid convection scheme, a gray radiative transfer module using a Rosseland-mean gray opacity scheme, and a simplified cloud tracer module that includes thermal feedback and scattering due to clouds. Our implementation follows the approaches set out by Showman et al. (2019), Tan & Showman (2021a), Tan (2022) and Lee et al. (2023), and adapts them for the THOR GCM, as detailed in the following subsections. A quick overview of our parameters for these additional physics modules can be found in Table 2.

2.4.1 Radiative transfer

We extend THOR with a computationally efficient, gray radiative transfer module tailored to field brown dwarfs, adding the ability to couple cloud-opacity feedback while eliminating costly stellar-insolation calculations that are irrelevant for our nonirradiated targets. To achieve this, we prune the original picket-fence scheme (Noti et al. 2023), retaining a single IR band for the planetary emission. This makes the scheme well-tested and efficient while satisfying our requirements for a temperature- and pressure-dependent local opacity structure that can couple to the cloud module we detail in Section 2.4.4.

Our method follows a classic two-stream, short-characteristics formalism with linear interpolants (Olson & Kunasz 1987) to solve the radiative transfer equation. At every timestep, the local gas opacity is interpolated from the Freedman et al. (2014) grid, which provides analytical fits to the Rosseland-mean opacity κ_R, defined as $$\frac{1}{{\kappa }_{\text{R}}(T,p)}={\left({\displaystyle {\int }_0^{\infty }\frac{\partial B_v}{\partial T}}dv\right)}^{-1}\left({\displaystyle {\int }_0^{\infty }\frac{1}{{\kappa }_v}}\frac{\partial B_v}{\partial T}dv\right).$$(14)

The Rosseland mean weighs the frequency-dependent opacity κ_ν by the temperature gradient of the Planck function B_ν and as a harmonic mean it emphasizes the contributions from the least opaque (smallest opacity) windows. It therefore reproduces the correct flux in the optically thick layers, exactly where brown dwarf interiors inject their heat. Replacing the constant, that is, the gray opacity of the double-gray method and using a single mean opacity rather than tens of thousands of wavelength bins used for multiwavelength methods, we both preserve the essential T–p dependence and achieve this at a reasonable computational cost. It should be noted that while the Rosseland mean delivers an accurate depiction of the diffusive and optically thick limit, in the optically thin parts of the upper atmosphere, the Planck mean should be preferred (Marley & Robinson 2015). In our single-band gray scheme, we do not switch means; following prior work by Tan & Showman (2021b), Komacek et al. (2022) and Lee et al. (2023), we impose a gray lower limit of κ_R,min = 10⁻³ m² kg⁻¹.

We expect the gas-phase thermal scattering to be negligible (Lee et al. 2023), but cloud particles can contribute appreciably. For this reason, we implement a further modification (detailed in Section 2.4.6) to account for the scattering induced by the presence of clouds in our models and only activate it when clouds are present.

2.4.2 Thermal forcing of the atmosphere

In our curated sample, with no instellation present, the dynamical response of the atmosphere is triggered through variations in the internal thermal flux coupling to the brown dwarfs’ rotation. Conceptually, the physical scenario we aim to model corresponds to convective plumes that rise and sink at the RCB. Several authors (e.g., Zhang & Showman 2014) have modeled this mechanism and Showman et al. (2019) have outlined their approach that works through stochastic thermal perturbations introduced at the RCB, which we follow in our current work.

Briefly summarized, thermal perturbations are characterized through functional forms S (λ, ϕ, p, t) that are separated into the vertical S_v(p) and horizontal S_h(ϕ, p, t) terms, where λ is longitude, ϕ is latitude, p is pressure, and t is time. The vertical structure is designed to decay over two scale heights in order to keep the perturbations localized near the RCB. In our study, this level is set to a constant p_rcb = 10 bar, matching other recent modeling efforts (Lee et al. 2024).

Physically, the above-mentioned buoyant plumes are present at a much smaller scale than our grid scale (Freytag et al. 2010) and are introduced at the smallest resolved scale in our GCM setup. This implementation implicitly assumes that the smaller-scale plumes reorganize themselves into a continuous larger-scale structure corresponding to our grid scale (Showman et al. 2019). Finally, using their physically motivated description, we estimate the amplitude of the thermal perturbations as $$T_{\text{amp}}\approx \frac{\Delta z}{H}\frac{T_{\text{eff}}}{{\tau }_{\text{rad}}\sqrt{n_f}},$$(15)

where ∆z is the maximum amplitude of the displacement, H the pressure scale height, T_eff the effective temperature, τ_rad the radiative timescale and n_f the forcing wavenumber. This completes the physical picture as we have chosen the relevant length scale and the amplitude of the perturbations we would like to introduce. The horizontal structure, then, is represented as a Markov process, that is, a sequence of stochastic events that depend on the previous state $$S(\lambda ,\varphi ,p,t+\delta t)=rS(\lambda ,\varphi ,p,t)+\sqrt{1-r^2}F(\lambda ,\varphi ,t),$$(16)

which has its roots in turbulence studies (Scott & Polvani 2007; Showman et al. 2019). This parametrization provides the simplest stochastic representation, in the absence of a detailed physical model, consistent with the idea that convective events at the RCB persist for a short duration before being replaced by new plumes.

To give a conceptual explanation, Eq. (16) updates the horizontal forcing pattern by combining two contributions. The first term, rS (λ, ϕ, p, t) with $r=1-\frac{\delta t}{t_{\text{storm\hspace{0.17em}}}}$, retains part of the previous pattern, such that tstorm sets the timescale over which the forcing loses memory of its past. Since δt ≪ t_storm, we have 0 < r < 1, so the pattern evolves gradually from one step to the next. The second term introduces new variability through F(λ, ϕ, t) constructed as a sum of spherical harmonics with amplitude set by T_amp and a characteristic horizontal length scale set by the forcing wavenumber n_f. The prefactor $\sqrt{1-r^2}$ scales the strength of this stochastic contribution so that the forcing neither decays nor grows without bound over many updates. Given a sufficient number of iterations, this yields a continually refreshed pattern with the prescribed horizontal scale, intended to represent the aggregate effect of unresolved convective plumes on the resolved flow. Our chosen parameters match the Lee et al. (2024) study and are summarized in Table 2.

2.4.3 Mixing length theory

To model the convective interior of brown dwarfs, we introduce a mixing-length scheme after Marley & Robinson (2015), Joyce & Tayar (2023) and Lee et al. (2024). Mixing length theory (MLT) (Prandtl 1925; Böhm-Vitense 1958) provides us with a parametric approach for representing the convective mixing that happens on a sub-grid scale.

Physically, the modeled scenario can be described as a parcel of fluid, rising buoyantly in the vertical direction due to its excess temperature, in regions where the atmosphere is unstable to convection. The stability of a region of the atmosphere is assessed through the Schwarzschild (Schwarzschild 1906) criterion $$\Gamma \ge {\Gamma }_{\text{ad}},$$(17)

where $\Gamma =-\frac{dT}{dz}$ is called the lapse rate and ${\Gamma }_{\text{ad}}=\frac{g}{c_p}$ is the adiabatic lapse rate, in other words, the rate that the temperature would decrease with altitude without exchanging heat with its surroundings. The core assumption of mixing-length theory is that a displaced parcel exchanges heat with its surroundings after traveling a characteristic distance L, called the mixing length. The mixing length is approximated as a multiple of the scale height H as L = αH, where α is an order-unity constant. We adopt α = 1, following Lee et al. (2024), while noting that the methodology of choosing of this parameter is still actively debated (Smith 1998; Joyce & Tayar 2023). For a detailed discussion of how adjusting this free parameter influences the results, we refer the reader to Robinson & Marley (2014).

MLT is inherently a one-dimensional approximation, where the mixing is assumed to occur only in the vertical direction. The mixing between the rising fluid parcel and its new surroundings results in a convective flux $$F_{\text{conv\hspace{0.17em}}}=\frac{1}{2}\rho c_{p}wL\left(\Gamma -{\Gamma }_{\text{ad\hspace{0.17em}}}\right),$$(18)

where ρ is bulk density of the atmosphere, c_p the specific heat capacity under constant pressure and w the vertical mixing velocity (Marley & Robinson 2015) $$w=L{\left[\frac{g}{T}\left(\Gamma -{\Gamma }_{\text{ad}}\right)\right]}^{1/2}.$$(19)

Eq. (18) is sometimes written in terms of the convective eddy diffusivity, approximated as, K_zz = wL (Marley & Robinson 2015). The convective heat flux F_conv is used to update the local temperature tendencies in every computational cell according to $${\left(\frac{\partial T}{\partial t}\right)}_{\text{conv}}=\frac{-1}{\rho c_P}\frac{\partial F_{\text{conv}}}{\partial z},$$(20)

using a sub-timestepping approach. We adopt an MLT timestep of ∆t_MLT = 0.5 s. Once the cumulative MLT updates reach the main dynamical timestep, the resulting temperature correction is passed to the dynamical solver as an update to the heating rate Q_heat.

As described above, this approach yields convective heat fluxes and, in turn, a local and time-varying thermal eddy diffusivity coefficient $$K_{zz}=L^2{\left[\frac{g}{T}\left(\Gamma -{\Gamma }_{\text{ad}}\right)\right]}^{1/2},$$(21)

evaluated in each layer that meets the instability criterion. This usage of K_zz is local and time-dependent and differs from the globally averaged eddy diffusivity coefficients sometimes adopted in theoretical or modeling studies (e.g., Steinrueck et al. 2021; Tan 2022). Our approach yields a self-consistent mixing profile that can distinguish vertically separated convective cells. While not as accurate as convective-resolving models (such as Lefèvre et al. 2022), the scheme captures the key features of interior mixing that we can couple to our cloud tracer scheme detailed in Section 2.4.4. At the same time, the MLT routine can only capture mixing where the convective heat flux is nonzero (F_conv > 0), while the vertical velocity above a convective zone decays rather than dropping to zero. This “overshoot” effect has been noted by multiple studies relevant for brown dwarf atmospheres (Freytag et al. 2010; Lefèvre et al. 2022; Lee et al. 2024). We therefore extrapolate the velocity above the top (i.e., to pressures lower than p_top) of each convective zone with the prescription of Woitke et al. (2020) $$\ln w_{\text{ov}}=\ln w_{\text{top}}-\beta ·\max \left(0,\ln p_{\text{top}}-\ln p\right),$$(22)

where β is a free parameter that relates to Freytag et al. (2010) and adopted as β = 2.2, following the work of Lee et al. (2024) in our study, and w_top represents the vertical velocity at the top of a given convective cell, as determined using Equation (19). This overshoot prescription feeds directly into convective K_zz structure via K_zz,tot = K_zz + K_zz,ov where K_zz,ov = w_ovL, giving a smooth decline of K_zz in the stably stratified region above each convective boundary. The overshoot term K_zz,ov is applied solely as a correction for tracer mixing (detailed in Section 2.4.4) and does not affect the MLT temperature tendencies (through Q_heat), since F_conv = 0 in a stable region.

2.4.4 Modeling of cloud tracers

To probe the effect of clouds on the thermal structure and the ensuing feedback, we implement a simple tracer-based equilibrium cloud scheme. Our implementation closely follows the works of Tan & Showman (2019), Tan & Showman (2021a), Komacek et al. (2022) and Lee et al. (2024). The model evolves the following equations: $$\begin{array}{l}\frac{D{q}_v}{dt}=(1-s)\frac{\min \left(q_s-q_v,q_c\right)}{{\tau }_c}-s\frac{q_v-q_s}{{\tau }_c}\\ -\frac{q_v-q_{\text{deep\hspace{0.17em}}}}{{\tau }_{\text{deep\hspace{0.17em}}}}+\frac{1}{\rho }\frac{\partial }{\partial z}\left(K_{zz}\rho \frac{\partial q_v}{\partial z}\right),\end{array}$$(23) $$\begin{array}{l}\frac{D{q}_c}{dt}=s\frac{q_v-q_s}{{\tau }_c}-(1-s)\frac{\min \left(q_s-q_v,q_c\right)}{{\tau }_c}\\ -\frac{1}{\rho }\frac{\partial }{\partial z}\left(\rho \langle q_c{V}_s\rangle \right)+\frac{1}{\rho }\frac{\partial }{\partial z}\left(K_{zz}\rho \frac{\partial q_c}{\partial z}\right),\end{array}$$(24)

where q_v is the vapor mass-mixing ratio (MMR), q_c the condensate MMR, q_s the equilibrium saturation ratio, τ_c the condensation timescale, q_deep the deep reservoir MMR, τ_deep the deep mixing timescale and s the saturation coefficient. Conceptually, this is an equilibrium condensation model with a prescribed relaxation timescale applied to the tracers, while advection of tracers follows the bulk flow of the atmosphere passively. For the relaxation timescale, we follow Tan & Showman (2019), who prescribe τ_deep = 10³ s and demonstrate that varying τ_deep by orders of magnitude has little effect. The condensation timescale τ_c follows Lee et al. (2024) (see Table 2), who adopt τ_c = 120 s. This is also consistent with Lefèvre et al. (2022), who use a smaller value (τ_c = 10 s) and find that varying τ_c from 1 s to 100 s produces only small changes in the cloud layer depth. For the stability criterion, we adopt the convention of Lee et al. (2024), defining it in terms of the supersaturation ratio S, given by $$S=\frac{{\chi }_v}{{\chi }_s},$$(25)

Where χ_s, the equilibrium VMR, is computed as χ_s = min (p_vap(T)/p, 1), with p_vap(T) the saturation vapor pressure as a function of temperature. The saturation coefficient s is then defined as $$s=\{\begin{array}{ll}1,\hfill & \text{\hspace{0.17em}if\hspace{0.17em}}S>1\hfill \\ 0,\hfill & \text{\hspace{0.17em}if\hspace{0.17em}}S<1.\hfill \\ \frac{d{q}_v}{dt}=\frac{d{q}_c}{dt}=0,\hfill & \text{\hspace{0.17em}if\hspace{0.17em}}S=1\hfill \end{array}$$(26)

This ensures that phase changes only occur when the system is either supersaturated (condensation) or subsaturated (evaporation), and that no mass exchange occurs in equilibrium (S = 1). Inside the GCM, this set of equations is represented through a standard continuity equation with a source term in the following form: $$\frac{\partial \rho q}{\partial t}+\nabla ·(\rho vq)=\mathcal{S},$$(27)

where ρ is the bulk density of the background atmosphere, q represents the MMR of either the condensable vapor or the condensate of a given cloud species, v the wind field, and $\mathcal{S}$ the source and sink terms for each equation.

The left-hand side of Equation (27) has the same structural form as the entropy equation evolved by the THOR dynamical core and is solved using the same finite-volume scheme introduced in Mendonça et al. (2016) and Deitrick et al. (2020). A more detailed derivation of how Equations (23) and (24) relate to Equation (27) is given in Appendix A. The source terms, represented by the term $\mathcal{S}$ on the right-hand side, are integrated separately with a sub-timestepping approach. Specifically, the treatment of source terms is split into three components. The first part, the mass exchange between vapor and condensate phases, along with vapor replenishment from the deep reservoir is integrated using a third-order adaptive Bogacki–Shampine Runge–Kutta scheme (Bogacki & Shampine 1989). Similarly, the vertical mixing terms calculated with the K_zz values derived from the MLT routine are integrated using an adaptive Bogacki– Shampine routine. We use the third-order Bogacki–Shampine method for both cases because lower-order integrators were not stable at reasonable time step sizes. Additionally, its embedded second-order solution gives a cheap error estimate for adaptive time stepping, which avoids a conservative fixed time step that would bottleneck performance and saves the time otherwise spent searching for suitable step values. The treatment of the gravitational settling of condensate particles requires a more detailed approach and is discussed in the following subsection.

The physical scenario we model assumes that there is a deep (p > 100 bar) vapor reservoir that has a solar elemental composition. For each condensate, we set q_deep (see Table 2) using the solar number density ratios (n (El) /n (H)) of Asplund et al. (2021), taking the value for the molecule’s least abundant constituent, which serves as the bottleneck in forming the condensate. Our assumption is guided by the fact that most Y dwarfs are observed in the Solar neighborhood, and have metallicities close to solar Leggett et al. (2017). We adopt a solar-abundance deep reservoir when no stronger observational or first-principles constraints are available, following the conventions of Gao & Benneke (2018), Komacek et al. (2022) and Lee et al. (2024). At the beginning of our simulations, we set the condensate MMRs to a very small value, q_c = 10⁻³⁰, while the corresponding vapor MMRs follow their saturation curves but never exceed their respective q_deep. In other words, we assume the deep reservoir has already come into equilibrium with the overlying atmosphere before cloud formation begins.

2.4.5 Gravitational settling of clouds

Equations (23) and (24) provide a simplified framework that permits mass exchange between vapor and condensate. In equilibrium, the system maintains a delicate balance between the upward mixing from the deep vapor reservoir and the gravitational settling of cloud condensates. The resulting large-scale cloud structure then depends on the efficiency of vertical mixing, both advective and convective, throughout the atmosphere.

Because we do not use an explicit microphysical model, assumptions about the cloud particle size distribution must be made to estimate settling velocities. In the absence of a first-principles theory, we follow convention and adopt a prescribed log-normal size distribution (Ackerman & Marley 2001) $$\frac{dn}{dr}=\frac{N_0}{r\sqrt{2\pi }\ln {\sigma }_g}\exp \left(-\frac{\ln {\left(r/r_{\text{med}}\right)}^2}{2\ln {\sigma }_g^2}\right),$$(28)

defined by a median particle radius r_med and a fixed geometric standard deviation σ_g, to parameterize the bulk settling behavior of condensates. This can be integrated to yield the total number density $$N_0=\frac{3ϵq_c\rho }{4\pi {\rho }_c{r}_{\text{med}}^3}\exp \left(-\frac{9}{2}{\ln }^2\sigma \right),$$(29)

where ρ_c is the density of the condensate particle and ϵ = µ_c/µ_gas is the molecular weight ratio between the condensate and the background gas. In this framework, the effective cloud particle radius can be expressed as (Lee et al. 2024) $$r_c=r_{\text{med}}\exp \left(\frac{7}{2}{\ln }^2\sigma \right).$$(30)

This corresponds to the Ackerman & Marley (2001) expression (Eq. (13) in their work) with f_rain = 1 and α = 1. Having defined the size distribution of cloud condensates, we estimate the settling velocity through (Ohno & Okuzumi 2018; Lee & Ohno 2025) $$V_s=\frac{2\beta g{r}_c^2}{9\eta }{\left[1+{\left(\frac{0.45g{r}_c^3{\rho }_c\rho }{54\eta }\right)}^{2/5}\right]}^{-5/4},$$(31)

where g is surface gravity, η the dynamic viscosity, ρ_c the density of the condensate and ρ the bulk density of the background atmosphere. This expression is an extension of the standard expression for terminal velocity (Rossow 1978), derived assuming a Stokes flow, for a particle falling in a viscous flow. The empirical factor β, called the Cunningham slip factor (Cunningham 1910; Davies 1945) $$\beta =1+Kn\left(1.257+0.4e^{-1.1/Kn}\right),$$(32)

accounts for deviations from continuum behavior in the flow and is a function (Li & Wang 2003) of the Knudsen number $$Kn=\frac{\lambda }{r_c},$$(33)

where λ is the mean free path of the fluid, calculated as (Jacobson 2005) $$\lambda =\frac{2\eta }{\rho }\sqrt{\frac{\pi \mu }{8RT}},$$(34)

and η is the dynamical viscosity according to Rosner (2000) $$\eta =\frac{5}{16}\frac{\sqrt{\pi \mu m_u{k}_{\text{B}}T}}{\pi d_{{\text{H}}_2}}\frac{{\left(\frac{k_{\text{B}}T}{{ϵ}_{\text{LJ}}{\text{H}}_2}\right)}^{0.16}}{1.22},$$(35)

where m_u is the atomic mass unit, d_H2 the molecular diameter of H₂ and ϵ_LJ,H2 the Lennard-Jones potential of H₂. Lee & Ohno (2025) adopt the same expression, while they consider contributions H₂, H and He. In this study, we limit ourselves to H₂ for the sake of simplicity. Numerically, we calculate the vertical transport using a second-order explicit MacCormack (MacCormack 2001) integrator with a Koren (Koren & van der Maarel 1993) flux limiter.

2.4.6 Thermal feedback and scattering due to clouds

So far, we have described how cloud tracers are distributed throughout the atmosphere, generally following the bulk motion of the background gas, while accounting for mass exchange due to evaporation and condensation, mixing from convective motions, and gravitational settling of condensates. As mentioned in Section 1, previous studies (Morley et al. 2012) indicate that for Y dwarf atmospheres, clouds are expected to drive dynamical activity. This necessitates that we allow the cloud distribution to thermally feedback on the background atmosphere.

Physically, to capture the radiative impact of clouds, we must specify the optical depth and scattering properties of each cloud species. The optical depth is computed from specifying the extinction opacity of clouds κ_ext (m² kg⁻¹). The scattering properties of the cloud particles are described by the single scattering albedo, ω [–], which quantifies the fraction of incident radiation that is scattered, rather than absorbed by cloud particles, and the asymmetry factor, g [–], which characterizes the directional bias of scattering, that is, the ratio of forward to backward scattering. These quantities, however, are wavelength-dependent and require knowledge of the size distributions of cloud particles. To integrate them into our radiative transfer scheme, which uses a single thermal band for brown dwarf emission, we follow the procedure described in Lee et al. (2024), adopting the size distributions and parameters detailed in Section 2.4.5. We use the LX-MIE (Kitzmann & Heng 2018) code to pre-calculate the wavelength-dependent cloud properties and, assuming a lognormal size distribution, integrate over size to obtain the normalized attenuation coefficient α (m⁻²), the single-scattering albedo ω (–) and the asymmetry factor g (–). We then form Rosseland-mean values by weighting with ∂B_ν/∂T on a temperature grid and tabulate the resulting normalized attenuation coefficient α_R(T), single-scattering albedo ω_R(T), and asymmetry factor g_R(T), which we interpolate to local background gas temperature at runtime. This allows us to calculate, on the fly, the Rosseland mean cloud opacity of a single cloud species as $${\kappa }_{\text{R\hspace{0.17em}},\text{\hspace{0.17em}cloud\hspace{0.17em}}}=\frac{{\alpha }_{\text{R}}{N}_0}{\rho },$$(36)

where the number density, N₀ (m⁻³) is given by Eq. (29). We approximate the contribution of our cloud species’ thermal contribution by adding it to the background gas opacity, interpolated from the tabulated Freedman et al. (2014) opacities (as described in Sect. 2.4.1) $${\kappa }_{\text{R},\text{\hspace{0.17em}tot\hspace{0.17em}}}={\kappa }_{\text{R},\text{\hspace{0.17em}gas\hspace{0.17em}}}+{\displaystyle \sum _i^{n_{\text{cloud\hspace{0.17em}}}}{\kappa }_{\text{R},\text{\hspace{0.17em}cloud\hspace{0.17em}},i}}.$$(37)

It should be noted that, although used in several other studies (Tan & Showman 2021b; Komacek et al. 2022; Lee et al. 2023), this is an approximation of total opacity. We apply a further correction to account for scattering introduced by the cloud particles present at a given layer of the atmosphere, through the absorption approximation of Lee (2024). The idea behind this approximation is to account for scattering through a modification of the transmission function during our radiative transfer calculations. The total single scattering albedo ω_R,tot and the total asymmetry factor g_R,tot are calculated as, $${\omega }_{\text{R},\text{tot\hspace{0.17em}}}=\frac{{\omega }_{\text{R},\text{\hspace{0.17em}gas\hspace{0.17em}}}{\kappa }_{\text{R},\text{\hspace{0.17em}gas\hspace{0.17em}}}+{\displaystyle {\sum }_i^{n_{\text{cloud\hspace{0.17em}}}}{\omega }_{\text{R},\text{\hspace{0.17em}cloud\hspace{0.17em}},i}}{\kappa }_{\text{R},\text{\hspace{0.17em}cloud\hspace{0.17em}},i}}{{\kappa }_{\text{R},\text{\hspace{0.17em}tot\hspace{0.17em}}}},$$(38) $$g_{\text{R},\text{\hspace{0.17em}tot\hspace{0.17em}}}=\frac{g_{\text{R},\text{\hspace{0.17em}gas\hspace{0.17em}}}{\omega }_{\text{R},\text{\hspace{0.17em}gas\hspace{0.17em}}}{\kappa }_{\text{R},\text{\hspace{0.17em}gas\hspace{0.17em}}}+{\displaystyle {\sum }_i^{n_{\text{cloud\hspace{0.17em}}}}{g}_{\text{R},\text{\hspace{0.17em}cloud\hspace{0.17em}},i}}{\omega }_{\text{R},\text{\hspace{0.17em}cloud\hspace{0.17em}},i}{\kappa }_{\text{R},\text{\hspace{0.17em}cloud\hspace{0.17em}},i}}{{\omega }_{\text{R},\text{\hspace{0.17em}tot\hspace{0.17em}}}{\kappa }_{\text{R},\text{\hspace{0.17em}tot\hspace{0.17em}}}}.$$(39)

For each layer, we first apply the δ-M+ scaling of Lin et al. (2018), to remove the strong forward-scattering peak and keep the two-stream calculation accurate. The resulting modified single-scattering albedo and optical depth are denoted as ω_R,tot^⋆ and τ^⋆. The modified co-albedo is then $${ϵ}^{*}=\sqrt{\left(1-{\omega }_{\text{R},\text{\hspace{0.17em}tot\hspace{0.17em}}}^{\star }\right)\left(1-g_{\text{R},\text{\hspace{0.17em}tot\hspace{0.17em}}}{\omega }_{\text{R},\text{\hspace{0.17em}tot\hspace{0.17em}}}^{\star }\right)},$$(40)

which multiplies the layer optical depth to give the absorption optical depth τ_a = ϵ^∗τ^⋆. Our radiative transfer solver (see Section 2.4.1) then evaluates the transmission function $$\mathcal{T}(\mu )=\exp \left[-{\tau }_{\text{a}}/\mu \right],$$(41)

for every Gaussian angle µ. Lee (2024) conclude that this method provides a good approximation of the scattering behavior while coming at negligible computational cost.

With these modifications in place, initially, we turn cloud tracers on passively, allowing them to evolve without contributing to the thermal budget. After 250 Earth days, sufficient time for the dynamics to redistribute the tracers and for condensates to settle into statistically steady cloud layers, we switch on their thermal feedback. Starting this feedback later avoids spurious heating at the lower boundary. The reasoning behind this is that the initial state is saturated yet condensation has not yet occurred, tracer concentrations are initially greatest near the surface, and immediate coupling would artificially amplify their effect there.

3 Results

3.1 Atmospheric structure of the Y dwarf sample

Having constructed a broad suite of Y dwarf simulation cases, we now examine the steady-state behavior of each run after reaching statistical equilibrium. Figure 3 presents the full sequence of simulations, with each panel showing the effective radiating temperature T_eff, rad for a given combination of effective temperature T_eff and rotation period P_rot. The effective radiating temperature is obtained from the outgoing longwave radiation (OLR) flux via $F_{\text{OLR\hspace{0.17em}}}={\sigma }_{\text{SB\hspace{0.17em}}}{T}_{\text{eff},\text{\hspace{0.17em}rad\hspace{0.17em}}}^4$, providing a temperature-equivalent measure of the gray outgoing flux emerging from the photosphere. Across cases of the same effective temperature, it varies only slightly, typically by less than 1%. While equatorial features are prominent in all simulations, their detailed morphology differs. For T_eff = 400 K, faster rotators exhibit slightly enhanced equatorial flux, with patchy structures appearing at higher latitudes in all cases. For the slower rotating P_rot = 10 h, 20 h cases at T_eff = 400 K, 500 K, an oscillatory equatorial pattern emerges, with small-scale, high-latitude variability. The T_eff = 600 K cases appear more uniform from the equator to mid-latitudes, while the polar regions show increased patchiness. As expected, higher resolution and faster rotation generally lead to finer-scale atmospheric features, although the most significant jump is from the P_rot = 20 h to the faster rotators, which display similar characteristic length scales in the structures that end up forming. Although Figure 3 may suggest that many T_eff = 500 K, 600 K cases are uniform away from the equator, this impression arises because equatorial variability is much stronger than at higher latitudes, as the T_eff = 500 K; P_rot = 2.5 h case indicates, high-latitude vortices are present in all runs but are too weak to stand out against the equatorial signal. Both the T_eff = 400 K and 500 K cases develop oscillatory equatorial patterns, with the 500 K runs showing the largest contrast and strongest winds, whereas at T_eff = 600 K the equatorial variability weakens and the flow is dominated by broad zonal bands, with variability shifted toward the polar regions.

Overall, the effective radiating temperature maps in Figure 3 are close to horizontally uniform, with only very small peak-to-peak contrasts, while still exhibiting coherent banded and tilted structures indicative of three-dimensional circulation. The small horizontal variations in these maps are consistent with a weak dynamical regime, which we discuss further in Section 3.4.

Figure 4 illustrates the evolution of the temperature–pressure (T–p) profiles across the entire set of simulations. As discussed in Section 2.2, the final T–p profiles remain remarkably close to the initial conditions shown in Figure 2. A weak inversion forms early in the simulation (not shown), but the atmosphere settles into a structure featuring an adiabatic interior and a nearly isothermal upper atmosphere extending to pressures as low as p ≈ 10⁻⁴ bar. The profiles shown are global averages, but individual columns reveal negligible horizontal temperature variation, consistent with the absence of any imposed longitudinal or latitudinal asymmetry in thermal forcing.

3.2 Vertical mixing

We expect vertical mixing to be one of the main drivers of dynamics in Y dwarf atmospheres. In Figure 4, we overplot the T–p structure with the globally averaged eddy diffusivities K_zz, calculated according to Eq. (21) by the MLT routine, giving a measure of the convective mixing that is going on in the atmosphere, and the “advective” mixing that is calculated from the GCM’s own vertical velocity w value as $H(z){\sqrt{\langle w^2(z)\rangle }}_h$.

One of the first trends that can be pointed out is the fact that convective mixing seems to be consistent across every sample in a given effective temperature T_eff, independent of rotation period P_rot of a given case. This can be seen by observing each row separately. This matches the design in which each atmosphere is forced from below, as outlined in Section 2.4.2. The thermal forcing is governed by the amplitude of the imposed thermal perturbations, T_amp, which is fixed by the effective temperature (see Section 2.4.2) and does not vary with rotation period.

Further trends stand out across our sample, namely all of the T_eff = 400 K cases present a uniform, deep, convective column of K_zz = 10⁵ m² s⁻¹ extending from the base of the atmosphere to the p = 1 bar level, tapering off linearly due to convective overshoots, arriving at the minimum level K_zz,min = 10 m² s⁻¹ around p ≈ 0.1 bar. For the T_eff = 500 K and T_eff = 600 K cases, the convective mixing structure appears to be split into two regions, one extending from the base of the pressure domain to about p = 40 bar, followed by another disjoint convective zone that sits between 1 bar ≤ p ≤ 10 bar for T_eff = 500 K and 1 bar ≤ p ≤ 3 bar for T_eff = 600 K. In a similar study, Lefèvre et al. (2022) used a convection-resolving model and found that their K_zz estimates peak at ≈10⁶ m² s⁻¹ within convective zones but drop by several orders of magnitude above them. They caution that these very small post-convective values are not appropriate for GCMs given the lack of wind shear. Our overshoot parameter (β = 2.2) follows this suggestion and produces a more gradual decline of K_zz, so the profiles shown here should be read as an upper limit on mixing.

While the individual structure of the advective vertical velocity changes across our sample, all of them fall in a range of 10⁻¹–10 m² s⁻¹, indicating a regime characterized by weak advective mixing. This is in line with our expectations from these substellar atmospheres, where convective mixing is expected to dominate. In all of our examples, the RCB falls within the p = 1–5 bar level and sits below the pressure level of the photosphere, which is at about p_photosphere ≈ 0.7 bar for all the presented runs.

Fig. 3 Effective radiating temperature maps of individual runs at trun = 1000 d. The shown results are time-averaged over the last 10 days of the given run time. Run cases vary in effective temperature from Teff = 400 K to Teff = 600 K from top to bottom and correspond to a rotation period of Prot = 2.5–20 h from left to right.

Fig. 4 Vertical mixing profiles of all simulations. The T–p profile is overplotted with the convective mixing Kzz and the advective mixing contribution $H(z)\sqrt{{\langle w^2(z)\rangle }_h}$. The pressure extents of convective regions are shown as gray horizontal bands in each panel.

3.3 Cloud properties

Figure 5 shows, for each case, the thermal band photosphere, the globally averaged VMRs of vapor and condensate for every cloud species together with the equilibrium saturation VMR χ_s. The T–p profile is superimposed to allow direct comparison with the vertical-mixing plots in Figure 4. Different rows correspond to the models with T_eff = 400, 500, and 600 K and a fixed rotation period of P_rot = 20 h. Although different rotation periods, for the same effective temperature, result in varying horizontal structures within a given layer, the globally averaged vertical profile of clouds and their associated number density distributions remain nearly identical. That is why a single rotation period, P_rot = 20 h in Figure 5, suffices to outline the main behavior. In every panel, the vapor VMR curve, shown in red, follows the saturation curve in green up to the deep-reservoir value listed as χ_deep in Table 2, showing that the simulated atmospheres stay saturated with vapor throughout. The condensate VMR curve, plotted in blue, marks the altitudes where cloud condensates settle, forming a cloud layer. As discussed in Section 2.4.5, the vertical cloud profile results from a balance among convective mixing, gravitational settling, and the supply from the deep reservoir. The deep reservoir sets an upper limit, yet whether that limit is reached depends mainly on the strength of vertical mixing, which, in our simulations, is dominated by convection. This link becomes clear when one compares the mixing intensities in Figure 4 with the shapes of the condensate curves in Figure 5. Gravitational settling confines the cloud layer vertically, and thermal coupling has a similar regulating effect, though by a different route. When cloud thermal coupling is active, part of the internally deposited heat is absorbed and scattered within the cloudy regions, giving a slight warming that evaporates a small fraction of the local condensate and lowers its VMR accordingly.

For T_eff = 400 K, KCl and Na₂S clouds appear near p = 10 bar and extend upward to about p = 0.3 bar, to the level of the brown dwarf’s photosphere. MnS clouds occupy a similar pressure but extend deeper and have significantly weaker abundances. At this effective temperature, the saturation curve lies below the deep-reservoir mixing ratio even at the lower boundary, so the vapor mixed upward becomes supersaturated immediately. Condensate, therefore, collects at the base of the atmosphere, where our model removes it, and its mixing ratio can rise only as far as convective transport can redistribute it. For the T_eff = 500 K and 600 K models, the behavior of clouds is qualitatively similar, though the extent of the cloud layers becomes narrower as the effective temperature increases. While the cloud layers show some overlap with the photosphere level, their optical depths peak around the base of the illustrated condensate curve and remain well below unity (τ_cloud ≪ 1) across the whole atmosphere.

Figure 6 shows the globally averaged vertical profiles of dia-batic heating rate produced by the dynamical core. In each panel, we plot the volumetric heating rate Q_heat (W m⁻³) along with the cloud layers where condensate mixing ratios peak, as indicated in Figure 5, and the convective zones identified by a positive convective heat flux F_conv > 0 from the MLT routine. Convective overshoot enhances the mixing of cloud particles but does not contribute directly to Q_heat. The Q_heat values shown here are the source terms that enter the model’s energy equation and are therefore expressed per unit volume. They can be converted to a temperature tendency through dT/dt = Q_heat/(ρc_p), resulting in peak magnitudes comparable to the thermal perturbations imposed at p = 10 bar. In the statistically steady state, Q_heat departs from zero primarily in the deep layers where internal forcing and convection act, while above the thermal-band photosphere (p ≈ 0.7 bar) it remains close to zero in the global mean. Although cloud layers extend above the convective zones, they remain optically thin in our gray thermal band and do not coincide with comparable enhancements in Q_heat in these profiles. It should be noted that the cloud layers that are shown in Figure 6, do not all correspond to the same VMR, they simply indicate where cloud condensates settle.

Fig. 5 Globally averaged vertical profiles of cloud condensate (blue) and vapor (red) VMRs, together with the saturation vapor VMR (green), for three cloud species (columns, left to right: KCl, Na2S, MnS) and three effective temperatures (rows, top to bottom: 400–600 K). Each panel shows results at a simulation time of trun = 1000 d, averaged over the final 10 days, for a fixed rotation period of Prot = 20 h. The horizontal gray shaded area marks the photosphere for the single-band thermal emission.

Fig. 6 Vertical profiles of diabatic heating rates for the run cases with Teff = [400, 500, 600] K and Prot = 20 h at trun = 1000 d. Overlaid are convective zones and cloud condensate layers, intended to give an overview of the physical processes that contribute to the radiative balance of the atmospheres under study. The smaller panels on the right-hand side shows the extent of each overplotted layer.

3.4 Flow regime

Figure 7 shows the thermal Rossby number Ro_T plotted against the nondimensional radiative timescale τˆ_rad, following the conventions of Hammond et al. (2023). The axis ranges match that study, and the gray line marks the boundary separating two of their flow regimes, the upper one with small vortices and an equatorial jet, the lower dominated by radiative forcing and lacking jets (see their Figure 8). Each symbol apart from the red star represents one of the GCM simulations from this study, with Ro_T and ${\widehat{\tau }}_{\text{rad\hspace{0.17em}}}$ calculated from our model outputs using the definitions in Section 2.1 (see Eqs. (3) and (5)). We also include an additional GCM run with T_eff = 600 K, P_rot = 20h, log g = 3.5, designed to test whether reduced gravity enhances variability. This case is plotted in a lighter green to distinguish it from the standard T_eff = 600 K; P_rot = 20 h run and its variability implications are discussed in Section 4.3. All of our simulated cases lie below the plotted regime line, inside the radiative-forcing-dominated regime, and the plotted values also reveal an even more rotation-dominated regime than the theoretical estimates in Figure 1 had implied.

All points belonging to the same effective temperature (same color) fall on a straight line with a negative slope. This is a simple result of the normalizations adopted in Section 2.1 (also seen in Figure 1), as the rotation period increases (i.e., the angular velocity Ω decreases), the thermal Rossby number Ro_T ∝ Ω grows, whereas the dimensionless radiative timescale ${\widehat{\tau }}_{\text{rad\hspace{0.17em}}}\propto \Omega$ shrinks. This implies, for a constant temperature T, $R{o}_{\text{T}}\propto {\widehat{\tau }}_{\text{rad\hspace{0.17em}}}^{-2}$ and our results reproduce this scaling. Conceptually, small Ro_T and large ${\widehat{\tau }}_{\text{rad\hspace{0.17em}}}$ signal a flow that is strongly constrained by Coriolis forces and has ample time (in units of a rotation) to redistribute heat before radiation can act. The point corresponding to the Lee & Ohno (2025) study falls within our sample, albeit at a different point than our (T_eff = 400 K; P_rot = 10 h) case. We would like to note that the chosen pressure level corresponding to the photosphere affects the radiative timescale expression significantly, although not the thermal Rossby number. A higher photosphere (lower pressure) would push all of the plotted points to the left and account for the difference, since most of the discrepancy is on the x-axis.

Finally, we would like to connect these results to the ones presented in Figure 1. In the absence of the vertical profiles calculated by our GCM, we generated the values in Figure 1 using a thermal wind approximation, estimating the thermal Rossby number through thermal wind velocities. In that scenario, the angle of the slope is a function of the magnitude of horizontal winds in the atmosphere, acting as a proxy for dynamical activity. In contrast, the GCM analysis evaluates the geopotential directly (Eq. (3)) and reports Ro_T. If we were to plot, instead, the standard Rossby number from the GCM winds (U ∼ 1–10 m s⁻¹), the points would lie at Ro ∼ 10⁻⁷–10⁻⁶ and the corresponding slopes of points sharing the same effective temperature would be smaller.

Fig. 7 Thermal Rossby number RoT plotted against the non-dimensional radiative timescale ${\widehat{\tau }}_{\text{rad\hspace{0.17em}}}$ for the complete sample of planets discussed in this study. The plotted quantities are calculated according to the definitions given in Sect. 2.1 and at the approximate photosphere corresponding to p = 0.7 bar. Chosen examples are highlighted together with their OLR profiles overplotted in the figure. The dashed lines correspond to a boundary between circulation regimes identified by Hammond et al. (2023). The red star is included as a reference point, the same as in Figure 1, corresponding to a run case of (Teff = 400 K; Prot = 10 h; log g = 4.46) from a recent GCM study.

4 Discussion

4.1 Dynamical activity

Across our suite of Y dwarf simulations, the results place all cases firmly in the radiative-forcing-dominated regime of the Hammond et al. (2023) study. While this is consistent with the weak dynamical signatures seen in the effective radiating temperature maps and with the near-uniform horizontal temperatures in the T–p profiles, compared to our theoretical estimations in Section 2.1, we end up producing less dynamical activity in the modeled atmospheres. Within each T_eff sequence, decreasing P_rot produces more pronounced structures but does not organize flow into coherent jets. Flux variations, while very small, seem to be concentrated around the polar regions with a more uniform equatorial region. Over the course of our simulations, a quasi-oscillatory equatorial pattern develops (visible in the T_eff = 400 K, 500 K; P_rot = 10 h, 20 h cases), persists for a few hundred Earth days, and then transitions into a uniform zonal band structure. The associated flux variations, likely reflecting the imposed internal thermal perturbations, remain negligible. The general lack of strong circulation in our results can be attributed to the weak wave activity. All models exhibit characteristics of mid-to-high-latitude Rossby waves, produced by rapid rotation interacting with the latitudinal gradient of the Coriolis force. The primary indicator of the presence of Rossby waves is the hemispherically antisymmetric tilts, visible in variations of the photospheric temperature (Figure 3), seeded by our internal thermal perturbations. These tilted features can be observed in higher latitudes, indicating westward traveling waves, matching the characteristics of Rossby waves. In principle, the interaction of Rossby waves with the mean flow could excite and maintain zonal jets, as suggested by Tan & Showman (2021a), but the lack of such features suggests that the strength of the excited waves compared to radiative processes is not efficient enough to drive such activity.

4.2 Cloud radiative feedback

Our vertical mixing analysis indicates that convection dominates transport over the weak advective contribution throughout the photosphere and below. This dominance is reflected in the cloud distributions (Figure 5), condensate layers form where the saturation curve and deep reservoirs permit, but their amplitudes and vertical extents are shaped primarily by convective mixing and gravitational settling. If the convective diffusion terms (in Equations (23) and (24)), whose strength is determined by the K_zz values from our mixing-length scheme, is omitted and the dynamical core alone advects cloud tracers upward, the condensate VMR peaks at roughly the same pressures (Figure C.2) but attains values several orders of magnitude lower than those plotted in Figure 5 and is more vertically confined. This implies that, despite clouds not having a strong effect on our presented results, their coupling to convective mixing is crucial for the resulting dynamical structure, if clouds are expected to contribute significantly.

Our results show that clouds are replenished efficiently but do not seed larger-scale dynamical structures. Thermal coupling within cloudy layers slightly warms the local environment, yet these effects remain secondary in the global energy balance. Previous studies (Tan & Showman 2021b; Lee et al. 2024) had indicated that cloud heating could drive a self-sustaining evaporate–cool–recondense cycle and act as a source of dynamical activity. In our models, this effect is too weak to significantly influence the circulation. Notably, Tan & Showman (2021b) explicitly state that their simulation setup was intended to study a cloud-driven dynamical scenario rather than to represent a realistic parameter set.

We note that our results are obtained with a gray radiative transfer scheme, an approach that is widely used in brown dwarf and exoplanet GCMs (e.g., Tan & Showman 2021b; Malsky et al. 2024; Kennedy et al. 2025). This method does not explicitly resolve spectral windows that sample deep, cloud-bearing layers and may therefore underestimate cloud radiative effects in those bands. Given the low condensate abundances and cloud optical depths in our simulations, we expect this bias to be modest for the parameter space explored here. This implies that our inference of weak cloud feedback should be viewed as conservative rather than invalidated by the gray approximation. Moreover, building on Lee et al. (2024) with a more realistic cloud model, Lee & Ohno (2025) apply the framework to a Y dwarf and obtain the same dynamical regime as our models, as noted in our results. Taken together, our results do not contradict previous studies; rather, they extend that framework to Y dwarf conditions.

In all models, our convective regions coincide with the pressure ranges where the base of the cloud layers forms (Figure 6) and where we impose thermal perturbations. This likely explains the insignificant cloud thermal feedback, as convective heat fluxes dominate the energy transport in these layers. A further possibility is that cloud number densities remain too low to yield a significant contribution to the Rosseland-mean opacity. This interpretation is consistent with Figure 6, where cloud layers extend above the convective zones but the diabatic heating rates show no corresponding increase.

4.3 Variability

Together, these results suggest clear directions for a more dynamically active test case. Following the scalings in our methodology, the most promising configuration would combine a slow rotation rate with the hotter end of our sample, and lower surface gravity (such as log g = 3.5). Such a setup should result in longer radiative timescales, while the choice of higher effective temperatures more strongly seeds variability, thereby allowing dynamical activity to organize into larger-scale structures and enhance variability. In Figure 7, we explore such a case and observe that it shifts in the intended direction in the $(R{o}_{\text{T}}-{\widehat{\tau }}_{\text{rad\hspace{0.17em}}})$ space. While this indicates a more dynamically active regime (further discussed in the Appendix), we still remain in the radiative-forcing-dominated regime. Overall, this extension, while not exhaustive, supports the trend indicated by our main GCM runs: within the physical processes considered here, Y-dwarf variability naturally lies in this regime.

Another direction to enhance variability is to explore alternative cloud-particle size distributions, which could strengthen cloud–radiation coupling. Finally, incorporating disequilibrium chemistry, in future works, will be important to assess whether compositional feedbacks alter the thermal structure and the variability observed here.

5 Conclusions

We have presented a curated sample of twelve GCM runs spanning the Y dwarf regime (400 K ≤ T_eff ≤ 600 K), to map the dynamical behavior of these atmospheres across a first, parameter-space-wide study. The model setup is observationally motivated and couples a parametrized convection scheme with a simple saturation-adjustment cloud model, enabling a self-consistent assessment of how convective mixing and cloud formation interact to shape the flow, an effect studied in previous studies (e.g., Tan & Showman 2021b) on brown dwarf variability.

Our results can be summarized as follows: Salt and sulfide species (KCl, Na₂S, MnS) form patchy (Appendix C) and optically thin cloud layers close to the photosphere, but their radiative feedback, and thus, effect on the T–p profile is minimal. Across the grid, the circulation remains in the radiative-forcing-dominated regime (Hammond et al. 2023), with small-scale, off-equatorial features that do not organize into jets. The T–p profiles and thermal emission equilibrate within a few hundred days, and top-of-atmosphere temperature variations remain small, with the observed temporal evolution primarily reflecting the imposed thermal perturbations rather than strong feedbacks among included processes. In the presence of the sluggish dynamical activity that we observe, variability is primarily controlled by rotational modulation from long-lived vortices. Overall, these findings align with the theoretical predictions of Hammond et al. (2023), in that no large-scale dynamical features are present in this circulation regime, with the atmospheres primarily controlled by radiative processes.

Our conclusions are also in line with the recent results of Lee & Ohno (2025), despite the inclusion of more cloud species and a wider parameter range coverage, supporting the view that models of Y dwarf atmospheres produce little dynamical activity in general. In particular, cloud–radiative effects in this regime appear insufficient to substantially alter the large-scale dynamical state or to induce new dynamical features. At the same time, our single-band, gray radiative transfer with Rosseland-mean opacities does not account for spectral windows that could probe deep, cloud-rich layers and may therefore underestimate cloud radiative feedback in those bands. In our simulations, however, the condensate mixing ratios keep the clouds optically thin near the photosphere, so their radiative influence remains small. Our conclusions should thus be interpreted as applying to atmospheres with relatively low cloud optical depths and do not preclude stronger feedback in cases with a larger deep vapor reservoir.

These results suggest that the dominant mechanisms behind Y dwarf variability remain uncertain. Although our predicted inhomogeneities are slight, the sensitivity of JWST means that contrasts of a few Kelvin can affect the spectra and should be detectable. A more sophisticated modeling approach, with higher resolution dynamics, improved cloud microphysics and radiative coupling, and disequilibrium chemistry, will help clarify the picture, but progress depends on future JWST observations analogous to recent variability studies (e.g., Biller et al. 2024; Chen et al. 2025) of L–T dwarfs.

References

Appendix A Adaptation of the tracer equations

Our 2-species tracer-based relaxation model can be generally written as $$\frac{Dq}{Dt}=\widehat{\mathcal{S}},$$(A.1)

where $\frac{D}{Dt}=\frac{\partial }{\partial t}+v·\nabla$ is the material derivative, q is the MMR of the vapor or the condensate, v the wind field and $\widehat{\mathcal{S}}$ are the related source terms. We would like to rewrite these equations as a standard continuity equation with a source term: $$\frac{\partial \rho q}{\partial t}+\nabla ·(\rho vq)=\mathcal{S}.$$(A.2)

In this form, we can pass the tracer MMRs directly to the dynamical core (as explained in Section 2.4.4), since the entropy equation solved by the dynamical core has the same form and is readily discretized on the icosahedral grid. The source terms can then be integrated in the physics modules and passed to the dynamical solver as correction terms, for it to advect tracer concentrations accordingly. The tracer equations for the vapor and the condensate species are given by: $$\frac{D{q}_v}{dt}=(1-s)\frac{\min \left(q_s-q_v,q_c\right)}{{\tau }_c}-s\frac{\left(q_v-q_s\right)}{{\tau }_c}-\frac{q_v-q_{\text{deep\hspace{0.17em}}}}{{\tau }_{\text{deep\hspace{0.17em}}}},$$(A.3) $$\frac{D{q}_c}{dt}=s\frac{\left(q_v-q_s\right)}{{\tau }_c}-(1-s)\frac{\min \left(q_s-q_v,q_c\right)}{{\tau }_c},$$(A.4)

Expanding the material derivative terms, we can write, $$\frac{D}{Dt}\left(\rho q_c\right)=\frac{\partial \left(\rho q_c\right)}{\partial t}+v·\nabla \left(\rho q_c\right),$$(A.5) $$=\frac{\partial \left(\rho q_c\right)}{\partial t}+\nabla ·\left(\rho q_{c}v\right)-\rho q_c(\nabla ·v),$$(A.6) $$\Rightarrow \frac{D}{Dt}\left(\rho q_c\right)=\frac{\partial \left(\rho q_c\right)}{\partial t}+\nabla ·\left(\rho q_{c}v\right)+q_c\frac{D\rho }{Dt},$$(A.7) $$\Rightarrow \rho \frac{D{q}_c}{Dt}=\frac{\partial \left(\rho q_c\right)}{\partial t}+\nabla ·\left(\rho q_{c}v\right),$$(A.8)

whereby we have used the vector calculus identity $$\nabla ·(\varphi v)=\varphi (\nabla ·v)+v·\nabla \varphi ,$$(A.9)

and the mass continuity equation $$\frac{\partial \rho }{\partial t}+\nabla ·(\rho v)=0.$$(A.10)

We can rewrite the expanded material derivative in the same form as the above-mentioned tracer equations A.3 & A.4, yielding, $$\begin{array}{l}\frac{\partial \left(\rho q_v\right)}{\partial t}+\nabla ·\left(\rho v{q}_v\right)=\rho [(1-s)\frac{\min \left(q_s-q_v,q_c\right)}{{\tau }_c}-s\frac{q_v-q_s}{{\tau }_c}\\ -\frac{q_v-q_{\text{deep\hspace{0.17em}}}}{{\tau }_{\text{deep\hspace{0.17em}}}}],\\ \frac{\partial \left(\rho q_c\right)}{\partial t}+\nabla ·\left(\rho v{q}_c\right)=\rho \left[s\frac{q_v-q_s}{{\tau }_c}-(1-s)\frac{\min \left(q_s-q_v,q_c\right)}{{\tau }_c}\right]\end{array}$$(A.11)

Fig. B.1 Vertical profiles of diabatic heating rates for the run case Teff = 600 K; Prot = 20 h ; log g = 3.5 at trun = 500 d. Overlaid are convective zones and cloud condensate layers. The smaller panels on the right-hand side shows the extent of each overplotted layer.

Appendix B Higher variability case results

Here we present a quick summary of the “higher variability” case from the Discussion Section. The only change in the setup of the model is a lower gravity of log g = 3.5. As discussed in the main text, this increases the radiative timescale and allows dynamical structures to form. In this case an equatorial jet emerges with u ∼ 30 m s⁻¹, and the temperature field shows a broad, nearly uniform, and generally warmer (about 0.5 K) equatorial region with similarly broad, uniform bands at higher latitudes. This supports the trend in Section 3 from small high-latitude vortices to larger zonal band structure focused near the equator.

The convective mixing and the extent of the cloud layers are narrower, as shown in Figure B.1. Together with the more zonally uniform temperature field, this suggests that the increased variability is driven by atmospheric wave activity rather than cloud feedback, which would produce more complex, small scale variations.

Appendix C Cloud distributions

In the parameter space we explore in this study, cloud radiative feedback remains ineffective in driving the atmosphere. For completeness, we report the horizontal cloud distributions at the level where the cloud condensates settle. Figure C.1 displays the VMRs of condensable cloud vapors and condensates at a pressure of p ≈ 1 bar for our run case T_eff = 600 K ; P_rot = 20 h. The top row, showing the vapor VMRs, closely follows the underlying temperature field, with amplitudes differing by orders of magnitude due to differences in condensation curves. The cloud condensates also track the temperature structure but show slight variations. Both KCl and Na₂S clouds display granular structures with higher overall VMRs than MnS, which shows a more continuous structure. This may reflect the higher VMRs and thus number densities of KCl and Na₂S, which strengthen local evaporation, cooling, and recondensation cycles. VMR variations remain below 1%, consistent with the overall ineffectiveness of cloud thermal feedback in these simulations. While we present a single case here, the described behavior captures the overall trends in our simulated sample: cloud structures largely mirror the underlying temperature field, with the degree of granularity set by the condensate VMRs at the cloud bases. These results largely indicate a lack of strong feedback mechanisms. The cloud patchiness is dictated by temperature variations due to the long-lived vortices present in the atmosphere.

Fig. C.1 Example cloud distribution plot for Teff = 600 K; Prot = 20 h from a snapshot at a runtime of trun = 500 d at p ≈ 1 bar.

Fig. C.2 Globally averaged vertical profiles of vapor (red) and condensate (blue) VMRs, together with the saturation vapor VMR (green), for KCl, Na2S, and MnS (left to right) in the Teff = 600 K and Prot = 20 h case at trun = 1000 d averaged over the final 10 days. Solid curves show the fiducial model including convective mixing; dashed curves show an advection-only variant with convective diffusion omitted.

Lastly, we use the same run case as an illustrative example to highlight the differences between pure advection, in which the convective diffusion terms in Eqs 23 & 24 are omitted and tracers are transported only by advection in the dynamical core and our fiducial cases presented in the main body of the paper including convective mixing of the tracers. Figure C.2 shows a comparison between these two cases and, other than the omission of the overlaid T–p profile, is identical to Figure 5 As stated in Section 4, the condensate layers peak at similar pressures, but their peak VMR values decrease by orders of magnitude when convective mixing is removed, demonstrating that convective mixing dominates the replenishment of condensates for the parameter space explored here.

All Tables

All Figures

Fig. 1 Summary of our calculations, presented in Section 2.1. Highlighted are the cases with gravity log g = 4.5, a wave-driving efficiency of η = 10−3 and a vertical scale of $\Delta \tilde{z}=2$ at the approximate photosphere level pphotosphere = 0.1 bar. Other possible combinations are overplotted in fainter colors to illustrate the complete space spanned by the range of values considered in this study. Dotted lines mark the edges of minimum and maximum values. The red star corresponds to (Teff = 400 K; Prot = 10 h; log g = 4.46) and is from a recent GCM study, included as a reference point.

In the text

Fig. 2 Various T–p profiles overplotted to illustrate the initial conditions from which the simulations were spun up. The solid lines correspond to the κlw = 2.5 × 10−3 m2 kg−1 curves with the lighter shades of each color illustration curves corresponding to the same internal temperature Tint with varying long-wave opacities of κlw = [2.5 × 10−5, 2.5 × 10−4, 2.5 × 10−2] m2 kg−1. Also overplotted are condensation curves for the cloud condensate species KCl, MnS, Na2S. The shaded horizontal region (p = 0.1 – 1 bar) corresponds to the approximate photosphere.

In the text

	Fig. 3 Effective radiating temperature maps of individual runs at trun = 1000 d. The shown results are time-averaged over the last 10 days of the given run time. Run cases vary in effective temperature from Teff = 400 K to Teff = 600 K from top to bottom and correspond to a rotation period of Prot = 2.5–20 h from left to right.
In the text

	Fig. 4 Vertical mixing profiles of all simulations. The T–p profile is overplotted with the convective mixing Kzz and the advective mixing contribution $H(z)\sqrt{{\langle w^2(z)\rangle }_h}$. The pressure extents of convective regions are shown as gray horizontal bands in each panel.
In the text

Fig. 5 Globally averaged vertical profiles of cloud condensate (blue) and vapor (red) VMRs, together with the saturation vapor VMR (green), for three cloud species (columns, left to right: KCl, Na2S, MnS) and three effective temperatures (rows, top to bottom: 400–600 K). Each panel shows results at a simulation time of trun = 1000 d, averaged over the final 10 days, for a fixed rotation period of Prot = 20 h. The horizontal gray shaded area marks the photosphere for the single-band thermal emission.

In the text

Fig. 6 Vertical profiles of diabatic heating rates for the run cases with Teff = [400, 500, 600] K and Prot = 20 h at trun = 1000 d. Overlaid are convective zones and cloud condensate layers, intended to give an overview of the physical processes that contribute to the radiative balance of the atmospheres under study. The smaller panels on the right-hand side shows the extent of each overplotted layer.

In the text

Fig. 7 Thermal Rossby number RoT plotted against the non-dimensional radiative timescale ${\widehat{\tau }}_{\text{rad\hspace{0.17em}}}$ for the complete sample of planets discussed in this study. The plotted quantities are calculated according to the definitions given in Sect. 2.1 and at the approximate photosphere corresponding to p = 0.7 bar. Chosen examples are highlighted together with their OLR profiles overplotted in the figure. The dashed lines correspond to a boundary between circulation regimes identified by Hammond et al. (2023). The red star is included as a reference point, the same as in Figure 1, corresponding to a run case of (Teff = 400 K; Prot = 10 h; log g = 4.46) from a recent GCM study.

In the text

	Fig. B.1 Vertical profiles of diabatic heating rates for the run case Teff = 600 K; Prot = 20 h ; log g = 3.5 at trun = 500 d. Overlaid are convective zones and cloud condensate layers. The smaller panels on the right-hand side shows the extent of each overplotted layer.
In the text

	Fig. C.1 Example cloud distribution plot for Teff = 600 K; Prot = 20 h from a snapshot at a runtime of trun = 500 d at p ≈ 1 bar.
In the text

Fig. C.2 Globally averaged vertical profiles of vapor (red) and condensate (blue) VMRs, together with the saturation vapor VMR (green), for KCl, Na2S, and MnS (left to right) in the Teff = 600 K and Prot = 20 h case at trun = 1000 d averaged over the final 10 days. Solid curves show the fiducial model including convective mixing; dashed curves show an advection-only variant with convective diffusion omitted.

In the text