Atmospheric circulation of brown dwarfs and directly imaged exoplanets driven by cloud radiative feedback: effects of rotation

Observations of brown dwarfs (BDs), free-floating planetary-mass objects, and directly imaged extrasolar giant planets (EGPs) exhibit rich evidence of large-scale weather. Understanding the mechanisms driving the vigorous atmospheric circulation of BDs and directly imaged EGPs and its effects on their observed lightcurve variability and spectral properties is a pressing need. Our previous work has suggested a strong role of radiative cloud feedback on inducing a spontaneous time evolution in a simple one-dimensional framework. Yet the radiative cloud feedback in a three-dimensional (3D) dynamical framework remains unexplored for conditions relevant to these objects. Here we present a series of atmospheric circulation models that self-consistently couple dynamics with idealized cloud formation and its radiative effects. We demonstrate that vigorous atmospheric circulation can be triggered and self-maintained by cloud radiative feedback. The circulation is dominated by cloud-forming and clear-sky vortices that evolve over timescales from several to tens of hours. The typical horizontal lengthscales of dominant vortices are strongly constrained by the rotation, showing a linear dependence on the inverse of rotation rate with stronger rotation leading to thinner clouds. Domain-mean outgoing radiative flux exhibits variability over timescales of tens of hours due to the statistical evolution of storms. The circulation driven by cloud radiative feedback represents a natural mechanism generating significant surface inhomogeneity as well as irregular flux variability. Our results also have important implications for near-IR colors of dusty BDs and EGPs, including the scatter in the near-IR color-magnitude diagram and the viewing-geometry dependent near-IR colors.


INTRODUCTION
Active weather is likely common among brown dwarfs (BDs) as has been indicated by several lines of evidence. First, atmospheric circulation drives temperature anomalies and inhomogeneous cloud coverage on a global scale, which is responsible for the observed lightcurve variability of many L and T dwarfs (e.g., Artigau et al. 2009;Radigan et al. 2012;Buenzli et al. 2012;Apai et al. 2013;Buenzli et al. 2014;Wilson et al. 2014;Metchev et al. 2015;Yang et al. 2016;Leggett et al. 2016b;Miles-Páez et al. 2017;Apai et al. 2017;Manjavacas et al. 2017;Zhou et al. 2018;Vos et al. 2019;Eriksson et al. 2019;Lew et al. 2020;Zhou et al. 2020;Bowler et al. 2020;Hitchcock et al. 2020, see also recent reviews by Biller 2017 andArtigau 2018). Second, spectra and near-IR colors of many L E-mail: xianyu.tan@physics.ox.ac.uk dwarfs suggest the presence of thick clouds in their photospheres (e.g., Chabrier et al. 2000;Allard et al. 2001;Tsuji 2002;Burrows et al. 2006;Saumon & Marley 2008;Charnay et al. 2018). The existence of thick clouds indicates the presence of atmospheric dynamics against gravitational settling. Third, the abrupt transition from L to T dwarfs is likely caused by either a sudden change of cloud patchiness or thickness (Ackerman & Marley 2001;Burgasser et al. 2002;Knapp et al. 2004;Marley et al. 2010), for which atmospheric dynamics likely plays a role. Fourth, large-scale circulation provides a source of vertical mixing that helps to explain the inferred chemical disequilibrium in a wide range of BDs (e.g., Saumon et al. 2006;Stephens et al. 2009;Leggett et al. 2016aLeggett et al. , 2019Miles et al. 2020), especially in stratified atmospheres where convection does not play a direct role in mixing. Other techniques detecting or constraining the presence of global circulation of BDs include Doppler imaging (Crossfield et al. 2014), simultaneous tracking of near-IR and radio variability (Allers et al. 2020) and precise near-IR polarization measurements (Millar-Blanchaer et al. 2020). The latter three techniques may be extended to a larger sample of BDs in the future.
Extrasolar giant planets (EGPs) detected by the direct imaging technique so far are mostly young, hot and relatively distant from their host stars, and therefore they receive negligible stellar bolometric flux compared to their internal heat flux. Their atmospheric structure and dynamics are likely determined mostly by internal luminosity, and thus they fall into the same category of field BDs in terms of atmospheric characteristics. Similar to BDs, the near-IR spectrum and color of most directly imaged EGPs suggest the presence of thick clouds and possibly significant chemical disequilibrium in their photospheres (e.g., Currie et al. 2011;Barman et al. 2011a,b;Marley et al. 2012;Oppenheimer et al. 2013;Ingraham et al. 2014;Rajan et al. 2015;Chauvin et al. 2017;Stolker et al. 2020). Although current instruments place a high threshold for the detection of lightcurve variability of directly imaged EGPs orbiting bright stars, their planetary-mass, free-floating counterparts commonly exhibit lightcurve variability (Biller et al. 2015;Zhou et al. 2016;Vos et al. 2018;Biller et al. 2018;Schneider et al. 2018;Manjavacas et al. 2019;Miles-Páez et al. 2019).
Understanding the atmospheric circulation of BDs and directly imaged EGPs is a pressing need given the large body of observational constraints. Unlike planets in close-to-modest distances from their host stars whose atmospheric circulation is primarily driven by the strong horizontal differential stellar insolation, atmospheres of relatively isolated BDs and directly imaged EGPs lack horizontal differential heating from external sources, and their atmospheric circulation is driven almost solely by internal heat flux. So far, two major categories of sources have been proposed to drive the global circulation in the stratified, observable atmospheric layers.
The first is a mechanically forced scenario: convection interacts with the overlying stratified atmospheres and generates a wealth of waves and turbulence. These atmospheric eddies propagate upward and interact with the large-scale flows, driving large-scale circulation. There have been a few studies in this direction. Local hydrodynamic simulations by Freytag et al. (2010) show that gravity waves generated by interactions between the convective interior and the stratified layer can cause mixing and lead to small-scale cloud patchiness. Showman & Kaspi (2013) presented global convection models and analytically estimate typical wind speed and horizontal temperature differences driven by the absorption and breaking of atmospheric waves in the stably stratified atmosphere. By injecting random material forcing to a shallow-water system, Zhang & Showman (2014) showed that weak radiative dissipation and strong forcing favor large-scale zonal jets for BDs, whereas strong dissipation and weak forcing favor transient eddies and quasi-isotropic turbulence. Using a general circulation model coupled with parameterized thermal perturbations resulting from interactions between convective interior and the stratified atmosphere, Showman et al. (2019) showed that under conditions of relatively strong forcing and weak damping, robust zonal jets and the associated meridional circulation and temperature structure are common outcomes of the dynamics. They also demonstrated that long-term (multi months to years) quasi-periodic oscillations on the equatorial zonal jets, similar to the Quasi-biennial oscillation observed in Earth's stratosphere, can be driven by the thermal perturbations.
The second scenario is a thermally-driven mechanism linked to cloud radiative feedback . Clouds are critical in shaping the thermal structure, near-IR color and spectral properties of substellar atmospheres via large opacity loading (see recent reviews by Helling & Casewell 2014, Marley & Robinson 2015and Helling 2019. Cloud radiative feedback is similarly important in driving a vigorous global circulation and atmospheric variability. Imagine an atmosphere consisted of patchy clouds. In the cloudy region, less thermal radiation escapes to space from to top of the cloud where it is relatively cold, whereas, in the cloudless region, more radiation to space occurs from much deeper and hotter regions. The two regions will, therefore, experience extremely different vertical profiles of radiative heating and cooling, which will lead to horizontal temperature differences on isobars. These horizontal temperature contrasts will drive an overturning circulation that, in turn, can advect cloud condensate vertically, and in principle might be capable of maintaining the cloud patchiness. There are two distinct cloud radiative feedbacks. The first one involves interactions between cloud formation, cloud radiative feedback and small-scale convective transport of tracers. It can result in spontaneous variability of both the cloud and thermal structures locally in an atmospheric column that occupies a small horizontal area. This has been extensively demonstrated in  using a simple one-dimensional (1D), time-dependent model that couples radiative transfer, cloud formation, and small-scale convective mixing. In a large-scale sense, this generation of spontaneous variability is intrinsically 1D, without requirement of explicit 3D large-scale flows. The second feedback with the large-scale circulation requires intrinsically multi-dimensional flows. When clouds are advected by large-scale motions, and the large-scale motions are driven by the radiative heating or cooling associated with cloud properties, such a system could be linearly unstable, providing an energy source to drive the circulation (Gierasch et al. 1973).
The latter multi-dimensional cloud radiative feedback has never been investigated for giant planets in the fully nonlinear regime. There is a strong motivation to examine this cloud feedback in a full numerical model and explore its dynamical properties. In reality, this form of feedback might be dominant in certain cases, for example, for some L dwarfs or late T dwarfs in which clouds condense in the upper stratified atmospheres where convection does not directly provide mixing (e.g., Tsuji 2002;Burrows et al. 2006;Morley et al. 2012). In this situation, the intrinsic 1D variability driven by cloud feedback would not occur. The evolution of clouds and temperature anomalies would rely on the explicit large-scale flows.
In this study, we numerically investigate effects of the cloud radiative feedback on driving a large-scale circulation in the context of BDs and directly imaged EGPs using general circulation models that couple cloud formation and their radiative feedbacks. We tune the model parameters such that clouds form in the stratified layers and the 1D intrinsic variability would not occur, leaving us a clean environment to understand the multi-dimensional cloud radiative feedback.
Our models assume a Cartesian geometry, periodic horizontal boundary conditions and a constant Coriolis parameter f across the whole model domain-the so-called f −plane approximation where f = 2Ω sin φ measures the local vertical rotation rate of the thin atmosphere, φ is latitude and Ω is the rotation rate. BDs are likely rapid rotators. Doppler broadening of spectral lines (Reiners & Basri 2008) and rotation periods inferred from lightcurve variability (e.g., Artigau et al. 2009;Radigan et al. 2012;Metchev et al. 2015;Apai et al. 2017;Allers et al. 2020) indicate typical rotation period of less than two hours to slightly more than 10 hours for field BDs. Planetary-mass, free-floating giant planets and Directly imaged EGPs likely rotate rapidly as well (Snellen et al. 2014;Zhou et al. 2016;Bryan et al. 2018). From equator to pole, atmospheric dynamics of BDs and directly imaged EGPs is likely affected by a wide range of local f , and the key dynamical lengthscales that are closely linked to the local planetary rotation may vary substantially. Therefore, f will be systematically varied while keeping other parameters the same in our models to thoroughly investigate effects of rotation on the circulation driven by cloud radiative feedback.
The reason we start our investigations with a f −plane approximation instead of full global models is the following. Turbulence is horizontally homogeneous under the constant− f assumption, which is considerably simpler than that in the full spherical geometry in which the latitudinal-dependent f induces horizontal anisotropy in the turbulence (Rhines 1975;Vallis & Maltrud 1993). This strategy provides a clearer context to explore effects of varying rotation on turbulence and cloud formation. Indeed, in a long history of investigations of geophysical turbulence (e.g., De Verdiere 1980;McWILLIAMS et al. 1999;Arbic & Flierl 2003;Arbic et al. 2007) and Earth's tropical cyclones (e.g., Held & Zhao 2008;Zhou et al. 2014), theories and models using the f −plane assumption have yielded significant insights on the properties of turbulence before considering a spherical geometry. Dynamics in global models will be investigated in a further study.
We have three basic conclusions from this study: 1) vigorous atmospheric circulation can be triggered and self-sustained by the cloud radiative feedback, providing a mechanism for surface inhomogeneous and its short-term evolution for BDs and directly imaged EGPs; 2) typical horizontal lengthscales of storms are closely related to the Rossby deformation radius, which is linearly proportional to 1/ f when other parameters are the same; 3) the vertical extent of clouds decreases with increasing rotation rate. This paper is organized as follows. Section 2 introduces the numerical model. Section 3 describes results from models with varying Coriolis parameter f . In section 4 we discuss our results and implications, then finally draw conclusions in section 5.

GENERAL CIRCULATION MODEL
We solve the standard 3D hydrostatic primitive equations in pressure coordinates (see, e.g., Holton & Hakim 2012). These are standard dynamical equations used in dynamical meteorology for a stratified atmosphere with horizontal lengthscales that greatly exceed vertical lengthscales (for reviews, see Vallis 2006), as appropriate to the global-scale atmospheric flow in photospheres of BDs and giant planets. Two tracer equations representing condensible vapor and clouds are simultaneously integrated. The horizontal momentum, hydrostatic equilibrium, continuity, thermodynamic, and tracer equations governing condensible vapor and clouds are as follows, respectively, where v is the horizontal velocity vector on isobars, ω = dp/dt is the vertical velocity in pressure coordinates, f is the Coriolis parameter, Φ is the geopotential,k is the local unit vector in the vertical direction, ρ is the density, ∇ p is the horizontal gradient in pressure coordinate, d/dt = ∂/∂t + v · ∇ p + ω∂/∂p is the material derivative, θ = T(p 0 /p) R/c p is the potential temperature, p 0 = 1 bar is a reference pressure, R is the specific gas constant and c p is the specific heat at constant pressure. The ideal gas law p = ρRT is used for the equation of state for the atmosphere. R v represents a Rayleigh frictional drag applied to horizontal winds in the deep atmosphere to crudely represent the effects of momentum mixing between the weather layer and the quiescent interior where flows are likely to experience significant magnetohydrodynamic drag. The drag linearly decreases with decreasing pressure and takes the same form as Liu & Showman (2013): is a pressure-dependent drag coefficient, which decreases from 1/τ drag (where τ drag is a characteristic drag timescale) at the bottom pressure boundary p bot to zero at certain pressure p drag,top : In all simulations, we fix p drag,top to 5 bars, which is much deeper than cloud forming regions and does not directly affect cloud formation. Kinetic energy dissipated by the frictional drag is converted to thermal heating by the term R θ . Due to the unknown nature of interactions between the interior and the weather layer in BDs and EGPs, the characteristic drag timescale τ drag is treated as a free parameter to explore possible circulation patterns. For major sets of simulations presented below, we adopt τ drag = 10 5 s. We assume that the temperature at the model bottom boundary remains fixed during evolution, mimicking an atmosphere attached to a convective interior with a specific entropy that does not evolve over short timescales. As in , in calculating the net thermal radiative flux F we assume a grey atmosphere (with a single broad thermal band) in a plane-parallel, two-stream, multiple-scattering approximation. The radiative transfer equations in an absorbing, emitting and multiple scattering media with the δ-function adjustment for scattering are solved using the efficient numeral package TWOSTR (Kylling et al. 1995, their equations (7) and (8) are solved in our model.). The background model atmosphere uses a frequency-averaged gas opacity, the Rosseland-mean opacity κ R,g , from Freedman et al. (2014) with solar composition. The Rosselandmean opacity gives a good estimation of radiative fluxes in the optically thick limit. In this study, clouds can extend to the upper atmosphere where it is optically thin by the gas opacity alone. In these regions there is no good choice a priori for a single opacity in the grey approximation. Therefore, we impose a constant minimal opacity κ min in the whole atmosphere, such that the gaseous opacity is κ gas = max(κ R,g , κ min ). In most of the cases we chose κ min = 10 −3 m 2 kg −1 . 1 This value has been used previously for the thermal opacity of typical hot Jupiter's atmospheres (Guillot 1 We have tested additional models with κ min = 0, κ min = 10 −4 , and κ min = 10 −2 m 2 kg −1 . The resulting typical circulation pattern, mean temperature structure and outgoing thermal flux are qualitatively very similar among cases with κ min = 0, κ min = 10 −4 and κ min = 10 −3 m 2 kg −1 . This is expected because when κ min 10 −3 m 2 kg −1 , the overall ther-2010) whose atmospheric temperatures are close to those of BDs. The total atmospheric opacity κ is simply the sum of the gas and cloud opacities κ = κ gas + κ c , the latter being determined by instantaneous cloud mixing ratio as a function of time and location as will be described below.
q v is the mass mixing ratio of condensable vapor relative to the dry background air and q c is the mixing ratio of cloud particles. The terms δ(q v − q s )/τ c and (1 − δ) min(q s − q v , q c )/τ c are sources/sinks due to condensation and evaporation, respectively. q s is the local saturation vapor mixing ratio and τ c is a conversion timescale between vapor and condensates. We set δ = 1 when vapor is supersaturated and δ = 0 otherwise. The conversion timescale τ c is very short compared to dynamical timescales in conditions relevant to BDs and EGPs (Helling & Casewell 2014), and here we set τ c = 10 2 s which is slightly longer than a dynamical time step for numerical stability. Cloud and dust formation in substellar atmospheres is highly complex (see reviews by, for example, Helling & Fomins 2013;Helling 2019), and many models adopt the idealized, chemical equilibrium framework (see summaries in Ackerman & Marley 2001;Helling et al. 2008;Marley & Robinson 2015). Here we adopt an even more simplified approach-the saturation vapor mixing ratio q s is assumed to depend on pressure alone. In this way, dynamical cloud properties (horizontal lengthscale and vertical extent of clouds) depend on dynamics alone (more specifically, rotation). We specify a uniform pressure p cond , lower than which q s decreases with decreasing pressure and higher than which q s is arbitrarily large such that no condensation would occur: where q deep is the deep vapor mixing ratio. In this study p cond is assumed to be 0.5 bar, and q deep is assumed to be equal to 2 × 10 −4 kg/kg. Both choices are turned to satisfy conditions that convection is absent in the cloud forming regions while maintaining a vigorous circulation. The power law exponent of 3 in Equation (8) is chosen non-rigorously: as long as q s sharply decrease with decreasing pressure, the resulting cloud formation is representative enough for the purpose of investigating dynamics. We have tested additional models with power law exponents of 5 and 7, and they show qualitatively similar results. The term Q deep = −(q v − q deep )/τ deep represents replenishment of condensable vapor by deep convection over a characteristic timescale τ deep , which is applied only at pressure deeper than 5 bars. τ deep is generally set to 10 3 s, broadly consistent with mixing timescales over a pressure scale height near the upper convective zone (Showman & Kaspi 2013). We assume a constant cloud particle number per dry air mass N c (in unit of kg −1 ) throughout the atmosphere, then use N c to determine local cloud particle size with given time-and locationdependent amount of condensate q c . Cloud particles are assumed to have a single size locally in each grid point, and the particle size mal structure is determined by the combination of the cloud opacity and the Rosseland-mean gas opacity. Cloud opacity usually dominates over the gas opacity, which is why clouds can drive circulation. The case with κ min = 10 −2 m 2 kg −1 exhibits a slightly different mean thermal structure because κ min is large enough to affect the overall temperature-pressure profile even without the cloud opacity. We concluded that as long as κ min 10 −2 m 2 kg −1 our results are not sensitive to the choice of κ min . r c is then determined via where ρ c is the density of condensate. Radiation interacts with cloud particles by absorption and scattering, which we parameterize by the extinction coefficient Q ext , scattering coefficient Q scat and asymmetry parameterg. As in , the total cloud extinction opacity given a particle size κ c (r c ) is obtained by averaging over all wavelengths λ using the Rosseland-mean definition: where B λ is the Plank function and κ ext (λ) = N c πr 2 c Q ext (r c , λ) is the cloud opacity at λ. Assuming spherical particles, the Rosselandmean Q ext , Q scat andg as a function of temperature and particle size are calculated with Mie theory using the numerical package written by Schäfer et al. (2012) 2 . We pre-calculate tables containing these parameters as functions of temperature and pressure, and read them into the GCM using linear interpolation during the simulations. In this study, we use enstatite (MgSiO 3 ) to represent properties of the cloud species, including a density ρ c = 3190 kg m 3 and the refractive index of enstatite obtained from Jäger et al. (2003). The results are not sensitive to the choice of cloud species because the essence of clouds is generating radiative heating/cooling that drives the dynamics. In all simulations we assume N c = 10 11 kg −1 , which means with the specified deep vapor abundance the condensed particle size is around 0.5 µm, consistent with the expected sub-micron particles in L dwarfs (Burningham et al. 2017). With a given particle size, the cloud settling velocity V s as a function of pressure and temperature is calculated using Eqs. (3)-(7) of Parmentier et al. (2013).
The deep layers of our models reach the convectively unstable region. We parameterize effects of rapid convective mixing using a simple convective adjustment scheme as in the NCAR Community Atmosphere Model (Collins et al. 2004, see their Section 4.6). If any adjacent two layers within a vertical atmospheric column are unstable, they are instantaneously adjusted to a convectively neutral state while conserving total sensible heat ∆pT, where ∆p is the layer thickness in pressure. In a single dynamical step, the whole column is repeatedly scanned until convective instability is eliminated everywhere. Tracers are also well homogenized within the adjusted domain during adjustment. There is no adjustment in the horizontal direction.
The dynamical equations (Equations 1-6) are solved using an atmospheric general circulation model, the MITgcm (Adcroft et al. 2004, see also mitgcm.org). A standard fourth-order Shapiro filter is applied in the horizontal velocity and temperature fields to maintain numerical stability (Shapiro 1970), which smooths out grid-scale variations but has minimal effect on the large-scale structure. The model is in a Cartesian geometry with constant f across the domain. We assume periodic horizontal boundary conditions. For most models presented below, we assume a horizontal domain size of 3 × 10 4 km. The resulting area is relatively small compared to the planetary surface assuming a Jupiter sized object such that the constant f approximation holds, while remaining large enough to ensure proper statistical analysis in some cases. The pressure domain is between 10 bars and 10 −3 bar, which is discretized evenly into 55 layers in log pressure. Temperature at the bottom boundary (10 bars) is fixed at 2600 K, resulting in atmospheric temperatures comparable to those of L dwarfs. For the major models shown below, horizontal resolution is 150 km or 100 km per grid cell depending on f , which is sufficient to fully resolve dynamics within a Rossby deformation radius. We adopt physical parameters relevant for BDs, including the specific heat c p = 1.3 × 10 4 Jkg −1 K −1 , specific gas constant R = 3714 Jkg −1 K −1 , and surface gravity g = 1000 ms −2 .
In section A of the Appendix, we perform a sensitivity test of horizontal resolution on a 3D model with f = 4 × 10 −4 s −1 . Statistical results of the test show good convergence, and that our current horizontal resolution is adequate to capture the essential dynamics.

Cloud radiative instability and initial evolution
The concept of cloud radiative instability which was first introduced by Gierasch et al. (1973) is particularly relevant in the clouddriven circulation of BDs and directly imaged EGPs. Large-scale dynamical instability may occur when the radiative heating depends on cloud properties, and when cloud properties depend on the large-scale vertical motion driven by the radiative heating rate. The essence of this instability can be illustrated using a simple linearized thermodynamic system as described in Gierasch et al. (1973). Let's assume that clouds are optically thick and that the atmospheric column was originally at rest and radiative equilibrium. The change of outgoing thermal flux could be due to brightness temperature deviations that are caused by either actual temperature variation or the cloud-top altitude variation. A slight perturbation of cloud-top altitude results in a change in the outgoing thermal flux. Let's suppose that the perturbation moves the cloud top upward to a colder altitude. The outgoing thermal flux decreases, and the atmospheric column is then no longer in equilibrium and radiative heating occurs. Large-scale vertical upwelling occurs in response to the heating, advecting the cloud top further upward (here the cloud settling speed is assumed to be smaller than the flow vertical velocity). This causes even less thermal flux emitted to space, which induces stronger heating and ascending, providing a positive feedback to the system. Gierasch's theory predicts an initial growth rate ∼ Γ c /(Γτ) for the vertical velocity, where τ = Gierasch et al. (1973) further showed that unstable linear modes are possible when the cloud radiative instability is coupled to the linearized dynamical systems. In the absence of rotation, 2D hydrostatic gravity waves have a set of pure unstable growing modes (no propagation) and sets of attenuating, eastward and westward propagating modes.
The linear model of Gierasch et al. (1973) serves as a valuable starting point. To better appreciate the transition from an initially linear to a nonlinear state, we start off our numerical investigation with an initial value problem. The simplest starting point is a twodimensional (2D) model in length-pressure domain with a periodic horizontal boundary condition and without planetary rotation. The system is initially at radiative-convective equilibrium in a cloudfree condition. Without clouds it will remain motionless (notice that effect of convection is parameterized as instantaneous adjustment to a convectively neutral state and thereby has no effect on driving large-scale dynamics in our models). We initialize a patch of clouds from 0.5 to 0.2 bar centered at the middle of the domain with an exponential decay in the form of q deep e −(x /1000 km) 2 where x is the distance from the domain center. Vapor concentration is initially set to q deep and is homogeneous at pressures greater than 0.5 bar but is zero at pressures less than 0.5 bar. The patch of cloud exerts substantial radiative heating at the cloud base and cooling at the cloud top, which subsequently initiate the circulation. Figure 1 exhibits the time and spatial evolution of cloud mixing ratio on the left column, radiative heating and cooling rates on the middle column and temperature perturbations relative to the initial temperature profile in the right column. The time sequence starts from time zero in the top row to 37.5 hours of simulated time in the bottom row as indicated in the left panels.
Interestingly, the dynamical evolution quickly differs from that predicted by the linear theory of Gierasch et al. (1973) wherein the cloud patch was expected to initially undergo exponential growth in thickness without wave-like propagation. Instead, radiation immediately drives a meridional circulation that splits the cloud patch into two parts that propagate in opposite directions. After the splitting, the cloud layer then starts growing linearly with time while propagating in both flanks. After about 10 hours, a secondary, thinner cloud layer develops at the center of the pattern, which is due to the convergent flow towards the center that uplifts condensable vapor. Eventually when the propagating clouds hit the model boundary after around 20 hours, the model symmetry breaks and the evolution becomes chaotic. The predicted linear growth rate Γ c /(Γτ) based on our modeled atmospheric conditions and a reasonable lapse rate for clouds is on the order of 10 −4 to 10 −3 s −1 , greater than the growth rate of cloud thickness while they propagate horizontally (from 2.1 to 18.8 hours in Figure 1). Although the linear theory predicts a continuous spectrum of unstable modes with similar growth rates, the growing patterns in our simulation are sparse and tend to be large-scale.
The horizontal propagation of the cloud pattern is significantly slower than adiabatic free gravity waves. A crude estimation of phase speed of the adiabatic wave is N H where N is the Brunt-Vaisala frequency and H is scale height, yielding ∼ 1000 m s −1 applying our model condition. This is much faster than the propagation of the cloud pattern, which is only roughly 250 m s −1 . A close inspection in the second row of Figure 1, especially in the heating rate and temperature perturbation panels, shows that a fast component already propagated more than 10 × 10 3 km horizontally away from the center in both directions at 2.1 hours. This fast component is roughly consistent with the dry waves that are triggered by the cloud radiative heating. Strikingly, in later times, vertical wave-like patterns that alternate with decreasing pressure are present above the cloud layers (see panels of heating rate and temperature perturbations at 6.2, 11.8 and 18.8 hours). These patterns do not propagate away like the fast component at 2.1 hours, but rather have a similar horizontal phase speed as the slowly moving cloud patterns and relatively stationary vertical phases-analogous to the forced, quasi-stationary waves.
The differences between our numerical results and the linear theory of Gierasch et al. (1973) in the non-rotating 2D model perhaps partly stem from the intrinsic nonlinearity of the diabatic system-the strong cloud radiative heating depends on all modes that advect the tracer, and is non-separable in some sense. This in-Figure 1. Time evolution of cloud distributions (left column), radiative heating and cooling rates (middle column) and temperature perturbations (right column). Each row represents the atmospheric state at different model time indicated in the left panels. These are results of the initial evolution of a non-rotating 2D simulation. The model is initially at rest, with an uniform initial cloud-free radiative-convective equilibrium. A patch of cloud is initially placed from 0.5 to 0.2 bar centered at the middle of the domain with an exponential decay in the form of q deep e −(x /1000 km) 2 where x is the distance from the domain center (as depicted in the upper left panel). The vapor concentration is initially q deep and homogeneous below 0.5 bar but is zero above 0.5 bar. The initial cloud patch generates radiative heating and cooling and drives subsequent evolution of the system. dicates that different linear modes can be "blended" together by the diabatic heating. This perhaps is why the growing and propagating cloud patterns in Figure 1 are analogous to the triggering of free gravity waves that propagate along opposite directions-the damped but propagating modes and the growing but non-propagating modes may be blended together, giving rise to the cloud patterns that simultaneously propagate in the horizontal directions and grow in the vertical direction. In addition, the nonlinearity also comes in when eddies grow sufficiently large. As seen in the middle panels in Figure 1, the strongest heating/cooling usually occurs at the edges of cloud patterns. Strong and complex local overturning circulations associated with the heating seem important in affecting the phase speed and growth rate of the pattern. Understanding these detailed complexities is outside of the scope of this study.

At statistical equilibrium
The circulation does not decay away but eventually evolves towards a chaotic state (see the bottom row of Figure 1). It is self-sustained and is in a statistical equilibrium in which gravitational settling of cloud particles is balanced by net upward transport of clouds and condensable vapor. The mean vertical extent of clouds does not change when cloud particles reach low pressures where gravitational settling can no longer be sustained by upwelling.
Available potential energy (APE) is generated from the spatially inhomogeneous radiative heating/cooling associated with partial cloud coverage, which is converted to kinetic energy (KE) associated with the vigorous circulation. The circulation in turn maintains patchy clouds that are responsible for generation of APE. Kinetic energy is removed mostly by the deep frictional drag, returning to the system via dissipated heat in the deep layers. The fixed bottom boundary temperature (in reality the hot interior of BDs or EGPs) ensures continuous energy supply that offsets radi-ation to space. The fact that clouds must eventually settle out at low pressures is essential to maintain the circulation-otherwise without settling, clouds will be eventually homogenized above the condensation level and therefore no patchiness would be available to generate APE. The concept of heat engines may be appropriate for the energetics of this system-it has been applied to similar topcooling-bottom-fueling Earth moist convection (Rennó & Ingersoll 1996) and tropical cyclone systems (Emanuel 1986).
The typical isobaric temperature perturbations in our model can exceed 100 K, and local horizontal wind speeds can exceed 1000 m s −1 . Clouds can be easily transported over 3 scale heights above the condensation level. Winds, temperature perturbations and clouds exhibit a wide range of spatial patterns, with horizontal lengthscales ranging about 1000 km to that comparable to the domain length, 3 × 10 4 km. Structures with lengthscales over 2 × 10 4 km are dominant in the vertical transport of clouds and kinetic energy, which evolve over a characteristic timescale of more than 10 hours. The seemingly dominant structure is easier seen in 2D simulations with extended horizontal domains up to 48 × 10 4 km shown in Appendix A. The dominant lengthscale for cloud patches is between 2×10 4 to 4×10 4 km, superposed with numerous smaller scale structure.
One could make use of a quasi-balance point of view to qualitatively understand why the dominant structure emerges at a horizontal lengthscale of a few 10 4 km. In statistical equilibrium, largescale subsidence through a hydrostatic, stratified atmosphere must be accompanied by diabatic cooling (an obvious example would be the Hadley circulation in Earth's troposphere, e.g., Pierrehumbert 2010). In a presumed quasi-stationary overturning circulation, air in the cloudy regions is heated and rises, and it must subsequently cool somewhere and descend due to the requirement of continuity. The rate of radiative cooling determines the velocity of subsiding air, which can then be used to constrain the typical horizontal lengthscale of the overturning circulation via continuity. Now we carry out a simple diagnostic scaling exercise based on the governing dynamical equations (Eq. [1] to [4]) to quantify the above picture. More detailed explanation of the scaling based on a similar set of dynamical equations can be found in, e.g., Komacek & Showman (2016) for application on hot Jupiter circulation. It is convenient to cast the primitive equations in log-pressure coordinates (e.g., Andrews et al. 1987;Holton & Hakim 2012). We assume statistical balance (thus the time-dependent terms disappear), and a symmetry between the ascending and the descending branch of the circulation, that they have similar speed, magnitude of heating/cooling rate and thus the fractional area. First, the continuity equation can be expressed as an order-of-magnitude estimation: where U is the characteristic horizontal velocity, L is a characteristic lengthscale of the circulation, W is a characteristic vertical velocity with a unit m s −1 , and H is scale height. In the angular momentum equation, because of the absence of rotation, the major balance is expected between the pressure gradient and advective forces, which in order-of-magnitude is U 2 /L ∼ ∆Φ/L (e.g., Tan & Showman 2017). Hydrostatic balance relates the isobaric temperature gradient to the isobaric geopotential gradient. The vertical difference of geopotential δΦ in a single column is δΦ = −RT δ ln p, where T is an appropriate vertically-averaged temperature and δ ln p is a characteristic thickness of the atmospheric column that is affected by the cloud radiative feedback in log pressure. Thus, the characteristic isobaric geopotential differ-ence between two columns that have a characteristic isobaric temperature difference ∆T is ∆Φ ∼ R∆T δ ln p. Combining with the angular momentum equation, we have The balance in the thermodynamic equation is expected to be between radiative heating/cooling and advection, and is the same for both the ascending and descending branch: δF where δF is the radiative flux difference inand-out of the atmospheric column with mass δM. The two terms on the right hand side are horizontal and vertical heat transport, respectively. Based on our simulated results, the term N 2 H 2 /R is usually larger than ∆T, so together with the continuity Equation (11) we expect that the balance in the thermodynamics is between vertical advection and radiative heating/cooling: Using all the balances in equations (11), (12) and (13), we are in the position to estimate the typical lengthscale of the overturning circulation We apply the following numbers based on diagnostic result of our simulations: (N H) 2 ∼ 10 6 m 2 s −2 , δF ∼ 10 5 Wm −2 , δM = δp/g ∼ 5 × 10 5 /10 3 kg (the thermal response of an atmospheric column to clouds extends far deeper than the assumed condensate level of 0.5 bar, and here we set δp ∼ 5 bars based on our numerical results), ∆T ∼ 100 K and δ ln p ∼ 4. The resulting characteristic lengthscale is about L ∼ 1.9 × 10 4 km, consistent with numerical results. In a short summary, in the absence of planetary rotation and under our assumed atmospheric conditions, a typical horizontal lengthscale on the order of ∼ 2 × 10 4 km is required for the hot air to cool and return to the altitude where it was heated up, which is a requirement for there to be a closed thermodynamic loop of the air. Failing to provide a sufficiently large simulated domain in a nonrotating 2D system may result in suppression of the self-maintained dynamics. This is because there is not sufficient room for the air to cool and descend, and thus no flow can be lifted up to form new clouds due to continuity and the initial clouds would gradually settle down. As a result, the response of the simulated atmosphere to the initial perturbation is simply to cool down as a whole towards the cloud-free equilibrium. Indeed, we tested a model with domain size 1.5 × 10 4 km (half of the default value), and the model eventually became quiescent no matter how strong an initial perturbation was applied.

Initial evolution
Planetary rotation plays a fundamental role in shaping the largescale dynamics of rapidly rotating BDs and directly imaged EGPs (Showman & Kaspi 2013). One of the profound consequences of rapid rotation on the stratified, thin atmosphere is the emergence of a major balance between the rotation and horizontal pressure gradient, which strongly affects the scale of flows, transport of tracers and typical horizontal temperature variations. Additionally, rotation leads to a natural dynamical lengthscale -the Rossby deformation radius R d = c/ f where c is the phase speed of gravity wave, a lengthscale over which many types of interesting dynamical phenomenons occurs (e.g., Vallis 2006).
In the inertia gravity wave system of Gierasch et al. (1973), when f is large, unstable modes associated with cloud radiative instability are only possible for small horizontal lengthscales. For the quasi-balanced flow under rapid rotation (the so-called quasigeostrophic flow), unstable modes are possible for axisymmetric modes, i.e., flows associated with coherent vortices. Growth rates of the axisymmetric modes are greater for larger horizontal wavenumber (smaller horizontal lengthscale). The instability ceases when the lengthscale becomes much greater than the deformation radius, because over large lengthscales the vertical motions tend to be inhibited by rotation.
Like the exercise done in Section 3.1.1, we investigate the initial-value problem to see how rotation shapes the dynamical evolution, but now in a 3D model with a double-periodic horizontal boundary condition and a constant Coriolis parameter f = 6 × 10 −4 s −1 . Similar to the 2D model, the 3D model is initialized with a patch of clouds from 0.5 to 0.2 bar centered at the middle of the domain with an exponential decay in the form of q deep e −(r/1000 km) 2 where r is the horizontal distance to the domain center. Vapor concentration is initially q deep and homogeneous at pressures larger than 0.5 bar but is zero at pressures less than 0.5 bar. The spatial and time evolution of cloud mixing ratio at 0.23 bar is shown in Figure 2, in which each panel shows the instantaneous cloud mixing ratio in color contours and the wind field in arrows. The time of each frame is indicated above each panel.
The initial evolution is essentially a geostrophic adjustment, a process by which an initially unbalanced perturbation is adjusted towards the geostrophic balance -a balance between the pressure gradient force and the Coriolis force (e.g., Holton & Hakim 2012;Gill & Donn 2016). Radiative heating associated with the initial cloud patch generates a strong positive temperature anomaly over a short time after initialization, which drives a subsequent outflow from the cloudy region. The outward winds are then deflected towards their right by rotation, and a strong anticyclone forms around the cloud pattern. 3 In this configuration, the Coriolis force tends to balance the outward pressure gradient force. Residual outward flow still persists despite the major geostrophic balance that is quickly established. This residual flow is responsible for the continuous horizontal growth of the cloud pattern until it saturates at about 4.5 × 10 3 km in radius. In the classic geostrophic adjustment problem wherein an initially unbalanced field is freely evolved without forcing and dissipation, the final equilibrium height (or pressure) field is predicted to be characterized by an exponential decay with a characteristic lengthscale of a deformation radius (e.g., Kuo & Polvani 1997;Gill & Donn 2016). The deformation radius in conditions relevant to our simulations may be estimated using a phase speed of long vertical wave c ∼ 2N H, yielding R d ∼ 3.3 × 10 3 km, which is roughly consistent with the prediction by the much simpler classical geostrophic adjustment theory despite the more complex dynamical process at work.
The circulation and cloud pattern before about 40 hours are dominated by a basic state that is axisymmetric around the center of the initial cloud patch. This is consistent with predictions by the quasi-geostrophic cloud radiative instability (Gierasch et al. 1973). The basic state is superimposed by a seemingly wavenumber-4 nonaxisymmetric component starting from 16.7 hour, and this particular component does not amplify as the evolution goes on. Eventually after about 50 hours, the basic axisymmetric circulation breaks, and multiple vortices emerge from a wavenumber-3 non-axisymmetric component at the end of the sequence shown in Figure 2. The initial non-axisymmetric component may be a result of instability of the anticyclone. One possibility is that a vortex initially embedded in an environment at rest may be unstable due to instability analogous to the shear instability (e.g., see a recent work by Reinaud & Dritschel 2019 for quasi-geostrophic vortices). This type of vortex instability is usually difficult to quantify in our primitive-equation models. The other possibility is the inertial instability, which may occur when the absolute angular momentum M = rV + 1 2 f r 2 of the vortex decreases with increasing radial distance ( f ∂M ∂r < 0, e.g., Holton & Hakim 2012), where r is the radial distance and V is the azimuthal velocity of the vortex. We have confirmed that the criteria f ∂M ∂r < 0 is indeed satisfied at the pressure level shown in Figure 2 very soon after the initialization. The unstable region is mostly near the outer edge of the anticyclone where the speed of the clockwise azimuthal wind rapidly increases with r.
Compared to the non-rotating results in Figure 1, we find that rotation significantly confines the horizontal extent of the circulation and cloud pattern. For example, at about 16.7 hours after the initialization, the edge of the cloud pattern only extends out to roughly 3000 km away from the center in the case with f = 6 × 10 −4 s −1 , whereas that in the non-rotating 2D system already reaches more than 12000 km away from the center. Our finding is qualitatively consistent with the conceptual understanding of the quasi-geostrophic, cloud radiative instability theory of Gierasch et al. (1973), showing that the typical horizontal lengthscale of the cloud-radiative-driven circulation can be limited to a scale close to the deformation radius.

Dynamics at statistical equilibrium with varying f
Similar to the 2D system, the 3D system shown in Figure 2 eventually evolves to a highly nonlinear and chaotic state in statistical equilibrium . that is self-sustained by cloud-radiative feedback. To systematically investigate the effect of rotation on the equilibrated dynamics, we first performed a set of 3D simulations with varying Coriolis parameter from f = 0 to 1 × 10 −3 s −1 over a fixed domain size (30000 km × 30000 km). Selected results of the equilibrated simulations with f = 0, 3 × 10 −4 , 6 × 10 −4 and 1 × 10 −3 s −1 are shown in Figure 3, in which the left column shows instantaneous cloud mixing ratio as 0.23 bar with horizontal wind vectors represented by arrows, the middle column shows the corresponding instantaneous temperature field with wind vectors at 0.23 bar, and finally the corresponding top-of-atmosphere thermal flux is shown the right column. Panels at different rows have different Coriolis parameter f as indicated above each panel.
The typical sizes of storms and vortices decrease with increasing rotation, as can be seen in Figure 3. When there is no background rotation ( f = 0), dynamics in the 3D simulation are characterized by domain-scale convergent and divergent flows which are qualitatively similar to those in the non-rotating 2D simulation. When rotation is included, the simulated domain at 0.23 bar is populated with anticyclones that are warmer and cloudy, and cyclones that are associated with relatively cloudless and cooler regions. Significant isobaric temperature variations over 100 K are present in all models. There exists significant variation of the isobaric cloud mixing ratio, ranging from almost zero to more than the deep mixing ratio q deep . Cloudy regions are usually accompanied with vigorous upwelling and cloud-free regions are associated with strong down-welling. The outgoing thermal flux exhibits strong variation across the domain, with more than 6 × 10 5 Wm −2 at cloudless regions and only slightly more than 3 × 10 5 Wm −2 at cloudy regions. This is consistent with our expectation that cloud opacity determines the level from which thermal flux escapes to space. Interestingly, anticyclones with cloud formation almost always have larger sizes than the cloudless cyclones.
Vertical thermal profiles are shaped by both the cloud-radiative effect and the dynamics. Figure 4 shows several instantaneous temperature-pressure (T-P) profiles sampled in the cloudy (red lines) and cloudless (black lines) regions of the model with f = 6 × 10 −4 s −1 in the top panels and those sampled in the model with f = 0 in the bottom panel. In the rapidly rotating case, the cloudy regions exhibit a characteristic cold top and hot bottom structure, whereas the cloudless regions are more isothermal above 1 bar which is much closer to the thermal structure that would occur in cloud-free radiative-convective equilibrium. All profiles merge to the same deep adiabatic profile that is specified by our bottom boundary condition. In the cloudy regions, clouds form above the condensation level (which is assumed universally to occur at 0.5 bar), and their greenhouse effect warms up the air below the condensation level. The greenhouse effect extends down to about 5 bars, below which the temperature profiles merge to the same deep adiabatic profile. Above the condensation level, the atmospheric lapse rate d ln T/d ln p is larger than that in the cloudless regions because the cloud opacity is large and results in an optically thick atmospheric column above the condensation level. Clouds usually extend more than one scale height above the condensation level, which, together with the large lapse rate, result in a characteristic crossing point between two sets of T-P profiles above the convective zone. Thermal flux escapes from cloudless regions at pressures of several bars where is hot, and thermal flux of cloudy regions is emitted from the much colder cloud-top.
Dynamics is important to determine the thermal structure seen in Figure 4 both explicitly and implicitly. Explicitly, temperature fluctuations with smaller vertical wavelength are caused by the inertia gravity waves that propagate upward. They are present in both cloudy and cloudless regions. Implicitly, although cloud opacity is the direct cause of the typical thermal structures, both the cloud and temperature anomaly structures need to resist fast-travelling gravity waves which efficiently remove anomalies. Rotation plays an essential role in two key ways. First, quasi-geostrophic balance due to rotation helps to sustain a large isobaric temperature difference (Charney 1963). Second, the resulting quasi-balanced motions are likely important for vertical transport of clouds, and they typically evolve much slower than gravity waves. Therefore, the vertical cloud structure can be sustained long enough for the radiative effects to be efficient. Indeed, T-P profiles in the non-rotating simulations (bottom panel of Figure 4) have much smaller differences between cloudy and relatively cloudless regions, and seem to be affected more by gravity waves.
The time evolution of the dynamical system is intriguing. Vortices are vulnerable and undergo straining, merging and dissipation over time. An individual anticyclone or cyclone usually evolves over a timescale of several to tens of hours. Inertia gravity waves are characterized by smaller length scales than the dominant vortices, evolving with a faster frequency and phase speed than the major vortices. The thickest cloud decks are usually formed in the mature anticyclones where vigorous upwelling helps to sustain them against gravitational settling. A growing anticyclone usually originates from a small perturbation that triggers a small patch of cloud. Then the geostrophic adjustment process discussed in Section 3.2.1 Figure 3. Instant horizontal cloud mixing ratio at 0.23 bar on the left column, and the corresponding temperature at 0.23 bar on the middle, and the corresponding outgoing top-of-atmosphere thermal flux on the right. Arrows represent instant horizontal wind vectors. These are results from models with different Coriolis parameter f = 0 (top row), 3 × 10 −4 (second row), 6 × 10 −4 (third row) and 1 × 10 −3 s −1 (bottom row). The model domain size is fixed at 30000 km × 30000 km. Other parameters are the same among these models. is driven by the cloud radiative feedback, injecting kinetic energy to and promoting the growth of the vortex. However, not all "seed" clouds are able to grow. Most of the time, small cloud patches are sheared apart or strained away by the turbulent flow, and only the lucky ones survive and are able to grow.
The size of the dominant vortices approximately linearly decreases with increasing f (Figure 3), and now we further quantify the horizontal size distribution of various quantities. Individual vortices are chaotic and unpredictable, and their statistical properties can be uncovered using power spectral analysis in wavenumber space. Assuming that the flow is isotropic in the horizontal direction, we can express the power spectra in the total wavenumber space |k| = |k x | 2 + |k y | 2 where k x and k y are wavenumber vector in x and y direction, respectively. Figure 5 shows time-averaged power spectra for K E = (u 2 + v 2 )/2 in the upper left, variance of temperature perturbations T 2 in the upper right, and variance of total tracer mixing ratio perturbations q 2 at the lower left, all of which are analyzed at pressure 0.23 bar. On the lower right we show the power spectra for the variance of outgoing thermal flux perturbations F 2 . These are from models with Coriolis parameter f = 1 × 10 −4 , 2 × 10 −4 , 3 × 10 −4 , 4 × 10 −4 , 6 × 10 −4 , 8 × 10 −4 and 1 × 10 −3 s −1 . The horizontal axis is in unit of 2π/wavelength. Note that these simulations are carried out with the same simulated domain size as shown in Figure 3. The energy containing wavenumber, which is defined as , where E(k) is the power at wavenumber k (e.g., Schneider & Liu 2009), are plotted as vertical dashed lines for each power spectra for K E, q 2 and flux variance.
The most prominent feature of the power spectra for KE and q 2 at 0.23 bar and F 2 is that the peaks of the spectra systematically shift from the smallest wavenumber (the longest wavelength) to larger wavenumbers (shorter wavelength) as f increases. The energy containing wavenumber k e also systematically increases with increasing f , although the k e sometimes differ slightly from the peak wavenumber. With a fixed f , peaks of power spectra for the above three quantities are almost the same. For simulations with f = 1 × 10 −4 and f = 2 × 10 −4 s −1 , their dominant horizontal structure is comparable to the domain size and therefore their power spectra peak at the smallest wavenumber. The power spectra of T 2 at 0.23 bar differ from those for the other three quantities, especially for those with f ≥ 3 × 10 −4 . Not only the spectral peaks of T 2 generally differ with those of the other three quantities given a fixed f , but also the T 2 power spectra show double local peaks when f is relatively large. The similar spectral shape between F 2 and q 2 but not between F 2 and T 2 quantitatively demonstrates that patchy clouds (and therefore the cloud-top temperature variations) are the dominant factors shaping the outgoing thermal flux variability, and that isobaric temperature variation is a side effect. Indeed, analysis of observed near-IR spectral time variability for several BDs (e.g., Buenzli et al. 2012;Apai et al. 2013;Lew et al. 2016Lew et al. , 2020, in which detailed static 1D atmospheric models were fit to the spectral varability, have shown that changes of cloud vertical structures are essential, and that temperature anomalies as well as the corresponding gas chemical variation help to improve the fit. KE spectra also help to infer the typical behavior of the turbulence. Pure 2D turbulence tends to transfer kinetic energy from small to large scale starting from where energy is injected-the so-called upscale KE transfer of 2D turbulence. This results in a characteristic KE power spectral slope of −5/3: E(k) ∝ k −5/3 in the so-called inertial range with lengthscales larger than the energy injection scale. In the meantime, transport of enstrophy (the square of potential vorticity integrated over the domain) from large to small scales leads to a KE power spectral slope of −3: E(k) ∝ k −3 in the inertial range where the lengthscale is smaller than the injection scale. Largescale turbulence in the rapidly rotating, stratified atmospheres has similar properties to the 2D turbulence (Charney 1971). For a comprehensive tutorial of incompressible and geostrophic turbulence related to atmospheric applications, see Chapter 8 and 9 of Vallis (2006). In the KE power spectra at 0.23 bar shown in Figure 5, for f 3 × 10 −4 s −1 , the spectral peaks roughly locate at the scale of the internal deformation radius L d = c g / f is, where c g is the phase speed of gravity waves. The KE power then decreases with increasing wavenumber in a characteristic slope close to −3 (a blue dotted line indicating a −3 slope is shown in the top left of Figure  5). However, at some point the KE spectra flattens out, then the slow increases again at very large wavenumber. For relatively small f , the spectral range with a −3 slope is smaller than that of larger f , and the transition to flatter spectra occurs at a wavenumber that is closer to the peak for cases with smaller f . For f 4 × 10 −4 s −1 , the KE spectral slope can be flatter than a −5/3 power law (a red dotted line indicating a −5/3 slope is plotted) whereas those with f 6 × 10 −4 s −1 is closer to −5/3. Interestingly, the q 2 spectra similarly exhibit a slope transition between −3 to −5/3 only when f is quite large.
The spectral features suggest a few interesting dynamical processes at 0.23 bar. First, KE at 0.23 bar is likely injected directly around the lengthscale close to the internal deformation radius via geostrophic adjustment. Then enstrophy of the turbulence cascades from the deformation radius to smaller scales, showing a characteristic KE spectral slope of −3. Turbulence in the slope= −3 range has a quasi-geostrophic nature. Second, in spectral space where the slope is flatter than −3, non-2D turbulence (possibly inertia gravity waves) becomes energetically important in the KE power spectra. It has been known that the Earth's upper troposphere exhibits a transition between a −3 KE spectral slope at large scale to a −5/3 slope at mesoscale (Nastrom et al. 1984). Non-2D turbulence likely contributes to this transition (e.g., Dewan 1979;Lindborg 1999Lindborg , 2007. In our simulations, inertia gravity waves may play this role. Indeed, flows with scales close to the deformation radius are close to geostrophic balance, indicating a quasi-2D nature. Flows with much smaller scales have large components of imbalanced inertia gravity waves. This was confirmed by investigating the relative fractions of the rotational and divergent parts of the horizontal velocity at different lengthscales. For models with smaller f , the degree of geostrophic balance is weaker, and inertia gravity waves are likely more important in the energetics than in models with larger f . Interestingly, the transition between slope −3 to −5/3 depends on the forcing amplitude as well. In Section B in the Appendix, Figure  B1, we show that when the model with f = 4 × 10 −4 s −1 is forced progressively weaker, the KE spectra at 0.23 bar exhibits a slope −3 all the way to large wavenumber before numerical dissipation takes over. When the wavenumber continues increasing toward the largest value, numerical dissipation near the grid scale becomes important and the spectral slope deepens. Third, at wavenumbers smaller than that of the deformation radius (larger lengthscales), KE is transferred upscale. However, at large scales, KE is strongly dissipated by the bottom frictional drag, so that the KE power spectra decreases as the wavenumber decreases. We shall discuss this point more in Section 3.4.
Fianlly we summarize dynamics with varying rotation in Figure 6 with various statistical quantities for simulations as a function of f (see the caption of Figure 6 for a description). There exist interesting trends with f . Mean temperature at 0.23 bar monotonically decreases with increasing f , while the mean outgoing thermal flux and mean cloud top pressure monotonically increase with increasing f . These three trends are tightly linked. As rotation increases, the mean thickness of clouds decreases, resulting in an increase of the cloud-top pressure and thus an increase of outgoing thermal flux due to the higher cloud-top temperature. The radiative greenhouse effect due to clouds is weakened by a reduction of clouds, and thus the temperature at 0.23 bar decreases with increasing f .
The wind speed at 0.23 bar shows a steep increase from f = 0 to f ∼ 1×10 −4 s −1 , and then only slight increase till f ∼ 3×10 −4 s −1 , after which it mildly declines with increasing rotation rate. This trend probably indicates the degree of geostrophic balance established in systems with a fixed domain size. At f = 0, dynamics are dominated by gravity waves. Although significant isobaric temperature differences must exist due to the cloud radiative effect, gravity Errorbars in some panels are the variance of the horizontal mean quantities with respect to time, and therefore they represent the time fluctuation of the instantaneous horizontal-mean field. All data in blue are from models with a fixed model domain size of 30000 km×30000 km. Panel (a) is for the horizontal temperature at 0.23 bar; panel (b) is for the RMS of horizontal wind speed at 0.23 bar; panel (c) is for the mean cloud-top pressure, which is arbitrarily defined as where the horizontal-mean cloud mixing ratio is equal to 2 × 10 −6 kg kg −1 ; panel (d) is for the outgoing topof-atmosphere thermal flux; and finally panel (e) is for the dominant vortex size defined as π/k max , where k max is the wavenumber with the maximum KE power. Note that in panel (e), black circles for f = 1 × 10 −4 , 2 × 10 −4 and 3 × 10 −4 s −1 are from models with domain size properly extended, which is to properly capture the statistical vortex behaviors. In panel (e), the dashed curve is a fit to the vortex size using the deformation radius L d = c g / f with a fitted gravity wave phase speed c g , yielding c g = 2170 ms −1 .
waves are still efficient in removing large horizontal temperature gradients across most of the domain (see Figure 4). With a small f , a finite deformation radius emerges and vortices form which helps sustain larger isobaric temperature variation than the non-rotating case. In this case the deformation radius is larger or comparable to the domain size, such that slight increases of f greatly help sustain a larger temperature gradient. This may help to explain the rapid trend when f is small. When f is sufficiently large and so the deformation radius is smaller than the domain size, more vortex pairs are populated in the domain, and the trend of increasing temperature gradient flattens. When f is further increased, the deformation radius becomes a small fraction of the model domain. In this regime, the horizontal wind speed scales as U ∼ (R∆T δ ln p)/(L f ) where L is a characteristic length scale of the vortices. L is loosely proportional to the deformation radius and thus L f remains a rough constant. Therefore wind speed scales with the characteristic horizontal temperature variation. As has been shown, cloud radiative forcing declines with increasing f , and as a result, the decrease of ∆T then slows the wind speed with increasing f . As we expected, when the deformation radius is much smaller than the domain size, the horizontal scale of dominant vortices linearly depends on 1/ f (see panel [e] of Figure 6). The fit using a deformation radius L d = c g / f fits well to the vortex lengthscale with a gravity phase speed c g = 2170 ms −1 . This is consistent with our previous estimate using a long vertical wave phase speed 2N H ∼ 2000 ms −1 .

Vertical tracer transport
One of the profound impacts of rotation on the dynamical system is that the thickness of cloud layers decreases with increasing rotation. We emphasize that the cloud condensation level is intentionally fixed at 0.5 bar at all horizontal locations, such that the change of mean cloud thickness is solely due to the change of rotation. The left panel in Figure 7 shows time-and horizontal-mean total tracer (gas and clouds) mixing ratio as a function of pressure for simulations with f = 0 to 1 × 10 −3 s −1 . The right panel contains the timemean RMS of the isobaric total tracer mixing ratio perturbation q as a function of pressure. The overall mean tracer mixing ratio smoothly decreases with increasing rotation rate, and is the same with the RMS q . Intuitively, one might imagine that the stronger the rotation, the greater the suppression of the vertical velocity due to the higher tendency towards geostrophic flow, making the flow less efficient to vertically transport tracers against gravitational settling. This results in smaller horizontal temperature anomalies via the weakened cloud radiative feedback, which then give raise to a positive feedback to the reduced vertical velocity. In this section we quantify the net vertical transport of tracers with varying f and at different pressure. Net upward transport of tracers across an isobar relies on the positive correlation of tracer abundance and vertical velocity, i.e., having upwelling at regions where tracers are more abundant and downwelling at regions with fewer tracers (e.g., Holton 1986;Parmentier et al. 2013;Zhang & Showman 2018;Komacek et al. 2019). In a statistically balanced state, the horizontal mean total tracer abundance is set by a balance between the net vertical transport by large-scale motions and gravitational settling of clouds: where q = q c + q v is the total tracer and A denotes a horizontal average of quantity A over the domain. Integrating Equation (15) from very low pressure where tracers are negligible to an arbitrary level where clouds are abundant, one obtains ω q = −q c V s at that level. This states that the net upward transport of total tracers across that isobar balances the total settling flux above that isobar. Regions with abundant clouds likely have strong radiative heating near the cloud base and cooling near the cloud top, whereas relatively cloudless regions usually radiatively cool. This typically results in upwelling at cloudy regions and downwelling at cloud-free regions, which naturally represents a mechanism for net upward tracer transport against gravitational settling. The question is, what type of motions are responsible for the major vertical transport and how do they depend on rotation? As seen from Figure 3 and their time evolutions, as well as the power spectral properties of KE and q shown in Figure 5, the size of cloud patches are comparable to the dominant vortex scales and are highly correlated to the flow pattern. Thus we would expect that the dominant tracer transport near the condensation level is by the mean flow associated with cloud-forming vortices. To quantify this, as a standard exercise in meteorology, we calculate the cospectral power density for the quantity ω q , which is simply 2R(q k ω * k ) where q k and ω k are the coefficients at wavenumber k space for tracer and vertical velocity, and ω * m is the conjugate of ω m (e.g., Randel & Held 1991). Figure 8 shows the cospectral power density at a pressure of 0.43 bar (right above the condensation level), for simulations with f = 0 to f = 1 × 10 −3 s −1 but all with a fixed domain size in the upper panel.
When f = 0 or the vortex scale is comparable to the simulated domain size, the divergent and convergent flow occurs on the domain scale, and therefore the dominant transport power is at the lowest wavenumbers. As f increases, similar to the KE spectra, the peak of the ω q cospectral power density is at lengthscales close to the deformation radius, and as a result it systematically shifts to larger wavenumber (smaller lengthscale). The peak value of the cospectral power density also monotonically decreases with increasing f , consistent with the weakened transport of tracers shown in Figure  7.
Dashed lines in the lower panel of Figure 8 show three additional models with extended domain sizes with for f = 1 × 10 −4 , 2 × 10 −4 and 3 × 10 −4 s −1 . Cases with f ≥ 4 × 10 −4 s −1 and with a fixed domain size are also plotted in the same panel for comparison. This confirms that the dominant transport mode is correlated with rotation even at small f but is not limited by the domain size. As expected, ω q spectra of the additional models also peak at lengthscales close to the deformation radius, and the peak value also monotonically decreases with increasing f .
At lower pressures, however, the vertical transport by smallerscale motions start to be comparable or even dominate over that by the vortex-scale motions. Figure 9 shows the ω q cospectral power density at different pressure levels for the simulation with f = 0 in the top panel, with f = 4 × 10 −4 s −1 in the middle panel and with f = 1 × 10 −3 s −1 in the bottom panel. In the lower two panels, one can see that at lower pressures, the cospectra starts to split into two groups, one being the vortex-scale motions near the deformation radius and the other being the much smaller motions. The total transport power by the small-scale groups is comparable to or larger than that by the vortex-scale motions in lower pressure.
There is likely a dynamical reason for this transition. As discussed in Section 3.1.2, tracer transport by a direct thermally-driven circulation requires that rising air is heated and subsiding air cools. This works well in the non-rotating cases, which is why tracer transport is always dominated by the lowest wavenumber at all pressure in the case with f = 0 (top panel of Figure 9). For large f , the direct thermally-driven mechanism works well only when the air in the ascending cloudy regions is not cooler than air in the descending cloudless regions. For example, in the case with f = 4 × 10 −4 s −1 , the temperature in cloudy regions is generally lower than that in cloudless regions when the pressure is less than 0.17 bar. In the middle panel of Figure 9 with f = 4 × 10 −4 s −1 , the total transport power by smaller scales starts to be comparable to that by the vortex-scale motions at 0.17 bar. At an even lower pressure, because the direct thermal-driven circulation is likely limited due to thermodynamic constraints, eddy-driven (here the eddy refers to flows with lengthscale smaller than the dominant vortex scale) circulation becomes increasingly important in the vertical tracer transport mechanism.

Dynamics with varying bottom drag
Kinetic energy of the flow at lengthscales larger than the deformation radius is strongly dissipated by the frictional drag that extends from bottom up to 5 bars. Why do flows driven by the cloud formation well above 5 bars experience bottom drag? In quasi-geostrophic turbulence, if KE is injected from the baroclinic flow (which refers to structures with vertical variation), the upscale transfer of KE also occurs over the vertical direction in the horizontal scale close to the deformation radius (e.g., Rhines 1977;Salmon 1978Salmon , 1980Smith & Vallis 2002;Chemke & Kaspi 2015). KE in the baroclinic flow will be transferred towards the flow with greater vertical wavelength, and eventually to the barotropic flow (which refers to flows independent of pressure or height). KE in the barotropic flow then continues to transfer to larger scales. In our case, KE is injected by the baroclinic cloud-radiative-driven dynamics. When it is transferred towards the barotropic flow, the KE is deposited mostly in the deep layer simply because there is more mass there. With a strong bottom drag that directly removes KE of the deep layers, the rate of KE generation from the upper cloud level is not sufficient to maintain a strong barotropic flow. Therefore, in Figure 5, the strength of the domainscale flow associated with the barotropic flows at 0.23 bar is rather weak.
To investigate the turbulence properties when the damping on the barotropic flow is different, we performed additional experiments assuming different bottom drag timescales of τ drag = 10 4 , 10 6 and 10 7 s. This is also motivated by the fact that the bottom drag only rather crudely represents the effect of mixing with the deep interior. There is no justification for how strong the drag should be, and here we test the sensitivity of our results to the varying drag timescale. We assume f = 1×10 −3 s −1 in order to maximizes scale separation between the deformation radius and domain size. Results are shown in Figure 10, in which temperature maps at 0.23 bar are shown on the left column, maps of cloud mixing ratio are shown at 0.23 bar on the middle column, and the outgoing radiative flux is shown on the right column. The drag timescales from the top to bottom row are 10 4 , 10 5 , 10 6 and 10 7 s, respectively. Results with τ drag = 10 4 and 10 5 s appear to be very similar in terms of the typical storm sizes, amplitudes of temperature perturbations and cloud mixing ratio. With strong drags, winds correlate well with the cloud forming vortices. When τ drag = 10 7 s, individual cloud-forming storms still have sizes comparable to the deformation radius. However, the wind field exhibits a strong component that is comparable to the domain size, much larger than the deformation radius (see wind vectors). The domain-size patterns can also be seen in the cloud mixing ratio and the outgoing thermal flux map. The domain-scale flow is largely pressure independent. Results of the model with τ drag = 10 6 s is roughly between those with strong drag and τ drag = 10 7 s.
We perform spectral analysis to quantify the scale separation of the turbulence. Figure 11 shows the same quantities in spectral space as in Figure 5 but for the four different τ drag . In the KE spectra, both case with τ drag = 10 4 and 10 5 show a single peak at 1.2 × 10 −6 rad m −1 , and then the KE decreases at both longer and shorter wavelengths. At larger wavelengths (smaller lengthscales), the case with τ drag = 10 4 s shows a lower energy density than the case with τ drag = 10 5 s. With longer τ drag , the KE spectra show two peaks, one at 1.2 × 10 −6 rad m −1 and a stronger one at the smallest wavenumber (domain size). The KE spectra at wavenumbers larger than 3 × 10 −6 rad m −1 are similar among the four cases. The local peak of KE spectral density near 1.2 × 10 −6 rad m −1 in cases with larger τ drag is smaller than that with smaller τ drag . Similar to the KE spectra, spectra of q 2 , T 2 and F 2 are all affected by increasing τ drag .
As the bottom frictional drag becomes weaker, the rate at whici it dissipates KE decreases, and the KE associated with the barotropic flow accumulates. Dynamics of the nearly barotropic flows are more akin to the 2D flow, and its KE in the f −plane can be transferred to the domain scale if the dissipation is weak. Indeed, the simulation with τ drag = 10 7 s shows a pair of cyclone and anticyclone comparable to the domain size. They appear as a pair because of angular momentum conservation. After the vortices grow to the domain size, they cease to grow in size but their KE keeps increasing until an equilibrium between KE upscale transfer and removal by drag is reached. In this situation, the KE power spectrum steepens towards the largest scale.
The domain-scale vortices exert feedbacks to the formation of storms on scales close to the deformation radius. The first prominent feature in Figure 10 is the existence of domain-scale isobaric Figure 10. Instant horizontal cloud mixing ratio at 0.23 bar in the left column, the corresponding temperature at 0.23 bar in the middle column, and the corresponding outgoing top-of-atmosphere thermal flux in the right column. Arrows represent instant horizontal wind vectors. These are results from models with different bottom frictional drag τ drag = 10 4 s (top row), 10 5 s (second row), 10 6 s (third row) and 10 7 s (bottom row). The Coriolis parameter is f = 1 × 10 −3 s −1 for all models. All other parameters are the same among these models. temperature difference, cloud patterns and outgoing thermal flux. The domain-scale cyclone is slightly cooler, less cloudy and emits more flux to the space, whereas the anticyclone is the opposite. Not like cloud formation with sizes near the deformation radius, these domain-size structure is not driven by the geostrophic adjustment. Effects of the bottom frictional drag-analogous to Ekman pumping and suction-acts against this configuration. For the nearly barotropic flow, the pressure gradient is mainly balanced by the Coriolis force associated with the wind in the absence of drag.
In the presence of frictional forces, the three-way force balance induces a drift in the direction from high pressure to low pressurethis means a convergent bottom flow that pushes air upward in the cyclone and divergent flow that sucks air downward in the anticyclone. This effect should in principle promote domain-scale cloud formation in the cyclone and suppress cloud formation in the anticyclone, which is the opposite to our results. The responsible mechanism is likely the migration of vortex under inhomogeneous background vorticity. Vortices tend to migrate to regions where the background vorticity matches the vortex's absolute vorticity. An important example is the migration of vortices under the planetary vorticity gradient-cyclones migrate poleward and anticyclones migrate equatorward (e.g., Adem 1956;LeBeau Jr & Dowling 1998;Scott 2011). In our model with τ drag = 10 7 s, the magnitude of relative vorticity anomalies associated with the domain-size vortices is comparable to f , and therefore the background vorticity is no longer homogeneous. Cloudlesss cyclones have relatively high absolute vorticity and cloudy anticyclones have relatively low absolute vorticity. When small cyclones are generated inside the domain-size anticyclone which has low absolute vorticity, they tend to migrate towards regions with high absolute vorticity-the domain-size cyclone. Whereas cloud-forming anticyclones tend to move to the domain-size anticyclone. Overtime, clouds tend to accumulate in the domain-size anticyclonic region and the domain-size cyclonic region is filled with cloudless air.
Secondly, the strong shear flows associated with the strong domain-size vortices decrease the efficiency of cloud formation. Newly formed storms near the strong shear regions are sometimes disrupted and sheared apart before they become mature. This likely causes smaller peaks around the internal deformation radius in the KE power spectra for models with τ drag = 10 6 and 10 7 s than in those with τ drag = 10 4 and 10 5 s (see Figure 11).
Clouds induce radiative heating to the domain-scale anticyclonic region, whilst relatively cloud-free air cools the domain-scale cyclonic region. This drives upwelling in the anticyclonic region and downwelling in the cyclonic region, providing extra tracer transport in addition to the vortex migration. As a result of cloud radiative feedback, the domain-scale vortices have temperature variations on isobars. This effect is prominent in the case with τ drag = 10 7 s. This is dynamically interesting because, by definition, the barotropic flow is independent of pressure (or height). Given that the bottom temperature is almost horizontally uniform in our model, a "true" barotropic flow should not have isobaric temperature variation at any pressure. In our case, the domain-size flow is only quasi-barotropic due to the cloud radiative feedback.

Observational implications
Most observed lightcurve variability of BDs and free-floating EGPs are thought to be caused by rotational modulation of surface inho-mogeneity (Biller 2017;Artigau 2018). The variability often appears to be periodic, indicating the presence of quasi-stationary surface features. However, a fraction of the variability exhibits changes over short timescales comparable to the rotation period, suggesting that the surface features evolve over timescales comparable or slightly longer than the rotation period (Artigau et al. 2009;Wilson et al. 2014;Metchev et al. 2015;Apai et al. 2017). A small fraction of them even show long-period (over 20 hours), non-sinusoidal variability (Metchev et al. 2015), and it is unclear whether rotational modulation is the typical cause. The stochastic evolution of storms driven by cloud radiative feedback helps to explain these irregularity in the lightcurve variabilities. Figure 12 shows the domain-averaged outgoing thermal flux as a function of time for three simulations with f = 2 × 10 −4 , 4 × 10 −4 and 8 × 10 −4 s −1 and with a fixed domain sizes of 30000km × 30000 km. The variability in these simulations are solely caused by the statistical fluctuation of the storm system. The case with f = 2 × 10 −4 s −1 exhibits the largest fractional flux variability up to about 7% peak-to-peak variation, while that with 8 × 10 −4 s −1 is the smallest among these three cases (about 2% peak-to-peak variation). The typical evolution timescale of the flux is tens of hours, longer than the overturning timescale of an individual storm (which can be roughly estimated as on the order of ∼ 10 4 s from the dominant lengthscales and wind speed of the storms). Figure 12 suggests that the more storms in the domain, the less effect that the fluctuation of a single storm would have in the total flux variability.
In the planetary surface, the latitudinal-dependence of the Coriolis parameter f indicates latitudinal-dependent sizes of local storms. Near the equator where f is small, storms might be able to grow to a lengthscale that is a nontrivial fraction of the planetary radius. For example, a BD with a 5-hour rotation has a Coriolis parameter at f = 1.2 × 10 −4 s −1 at 10 • latitude and f = 7 × 10 −4 s −1 at the pole. Given our model conditions, this results in a typical storm diameter of about 18 × 10 3 km at 10 • latitude and 3 × 10 3 km on the poles. Assuming a Jupiter radius, a single storm near the pole covers less than ∼ 0.1% of the disk whereas that at 10 • latitude covers ∼ 1.7% of the disk. If the BD or directly imaged EGP is observed from pole on, statistical fluctuation of storms is expected to cause negligible flux variability. However, if the view from equator on, not only rotational modulation is maximized, but the evolution of the larger storms may induce additional variability over a timescale on the order of tens of hours. This provides a reasonable explanation for the irregular variability of many observed lightcurves of BDs (e.g., Metchev et al. 2015).
When vertical transport of clouds is dominated by large-scale flows via cloud radiative feedback, one would expect that clouds are thicker at low latitudes and thinner at high latitudes for a BD with certain rotation rate. If the BD rotates sufficiently slow, the equator-to-pole cloud thickness variation may be small. Under similar atmospheric conditions -including temperature, gravity and metallicity, BDs that rotate faster are expected to have overall thinner clouds than those rotate slower. These consequences have implications for observed near-IR colors for BDs and directly imaged EGPs as the thickness of clouds affects the spectral properties. For example, different overall cloud thickness caused by different rotation rate may contribute to the near-IR color scattering of mid-to-late L dwarfs (e.g., Faherty et al. 2016). Recently, Vos et al. (2017) suggested that BDs that are viewed from more equator-on tend to exhibit redder near-IR colors and larger flux variability than those viewed more pole-on, for which our numerical results naturally support.

Unresolved issues and outlook
We have shown in Figure 7 that cloud thickness decreases with increasing rotation rate. One would intuitively expect that faster rotation can lead to a more stringent balanced state of the flow given roughly the same amount of cloud radiative feedback, and therefore weaker vertical motions to transport tracers against particle settling. However, this picture itself would not automatically imply a mechanism. This is a central unresolved issue of this study. Revealing the detailed mechanism likely requires an understanding of the eddies and their effect on the mean flows (here, eddies refer to motions with scales smaller than the vortex scale, and the mean flows refer to the vortex motions). If the vortex is roughly circular, the flow can be decomposed into axisymmetric and non-axisymmetric components. The axisymmetric flow will be largely in balanced state, with a residual unbalanced component providing divergence and convergence that are essential to drive cloud formation. The balance of this residual overturning circulation is primarily between angular momentum transport by eddies and the Coriolis force associated with the mean residual overturning flow if the vortex is close to quasi-geostrophic balance (e.g., Showman & Kaspi 2013). Therefore, the cloud-forming overturning circulation is essentially both thermally driven and eddy regulated. Furthermore, near the cloud top, as we have shown in Figure 9, eddies contribute significantly to the net upward transport of clouds. Cloud vertical transport near the cloud top is essential because the cloud-top pressure determines the outgoing thermal flux and therefore the strength of the column atmospheric heating. In summary, understanding eddies is essential to pin down the mechanism by which stronger rotation leads to vertically thinner clouds. Future models adopting various level of complexity are needed to investigate this mechanism.
Our idealized models are dedicated to a clear understanding of dynamical mechanisms. An obvious next step within this framework is to extend the model domain to a global geometry which includes effects of the latitudinal variation of f -the so-called β effect where β = df /dy. The β effect plays a central role in driving zonal banding and jets in Jovian and Saturnian atmospheres (see reviews by, for example, Vasavada & Showman 2005;Showman et al. 2018). This would help to clarify whether BDs and directly EGPs exhibit zonal jets like Jupiter and Saturn as indicated by long-term lightcurve monitoring (Apai et al. 2017) and simultaneous tracking of radio and near-IR flux variability (Allers et al. 2020). Future developments beyond the current idealized framework include adopting a realistic radiative transfer scheme coupled with chemistry and cloud formation-this allows direct comparisons between global models and the rich observed spectrum, near-IR colors and spectral time variability of BDs and directly EGPs.

CONCLUSIONS
Cloud radiative feedback likely plays an essential role in driving vigorous atmospheric circulation in self-luminous, substellar objects including brown dwarfs (BDs) and directly imaged extrasolar giant planets (EGPs). In this study, we have numerically investigated the atmospheric circulation in conditions relevant to these objects using a general circulation model that is self-consistently coupled with idealized cloud formation and its radiative feedback. As a first step in this line of study and to better understand the effects of rotation on the turbulent cloud formation, our models adopt a constant Coriolis parameter f for the whole domain in a Cartesian geometry with a double periodic horizontal boundary condition. We have reached the following key conclusions: (i) Vigorous atmospheric circulation can be triggered and selfsustained by interactions with cloud radiative feedback in conditions appropriate for BDs and directly imaged EGPs. In a local, constantf approximation, the circulation is dominated by turbulent vortices, with thick clouds forming in the anticyclones and thin clouds or clear sky in the cyclones. Local wind speeds can reach 1000 m s −1 and isobaric temperature differences can be a few hundreds of Kelvin with plausible physical parameters. Fractional outgoing thermal flux differences between cloudy and cloudless regions can be comparable to 1. This is a natural mechanism to generate significant surface inhomogeneity that is responsible for the observed lightcurve variability of BDs and EGPs.
(ii) In the presence of strong rotation and strong deep frictional drag (which crudely represents interactions between the weather layer and the deep quiesent interior), the characteristic horizontal lengthscales of dominant vortices are close to the deformation radius L d = c g / f where c g is the phase speed of gravity waves. In the absence of rotation or when rotation is weak, the circulation is characterized by convergent and divergent flows with horizontal lengthscales regulated by the radiative timescale of the atmosphere.
(iii) Stronger rotation leads to vertically thinner clouds, and the cloud thickness is greatest in the absence of rotation. Both the vortex-scale flow and smaller-scale inertia gravity waves contribute to the vertical transport of vapor and clouds. The mechanism by which rotation regulates the cloud thickness is likely related to interactions between vortices and smaller-scale flows.
(iv) When the deep frictional drag is weak, turbulence can grow to ever larger horizontal lengthscales that are comparable to the simulated domain size, whereas the sizes of storm-forming vortices remain close to the deformation radius. Strong domain-scale turbulence could affect efficiency of smaller-scale cloud formation and induce horizontal drifts of smaller-scale cloud-forming vortices.
(v) Turbulent storms driven by cloud radiative feedback evolve over timescales of several to tens of hours, and the statistical fluctuation of the ensemble of storms could induce variability of the domain-mean flux. The flux variability caused by the dynamical evolution of storms helps to explain the irregular lightcurve variability observed in a fraction of BDs. The change of cloud thickness with different f indicates that the cloud thickness may be different at different latitudes of the BDs or directly imaged EGPs. The global-mean cloud thickness might also be different for objects with different rotation rate. These may contribute to the observed scattered near-IR colors of dusty L dwarfs and the viewing angle dependent near-IR colors of BDs.   Figure 1 shows a test of the non-rotating 2D models with different domain length. Dominant aspects of the circulation, including the mean kinetic energy, mean cloud thickness and wind speeds are not sensitive to the chosen domain size in Figure 1. All simulations in Figure 1 exhibit characteristic large-scale structures with a size of about 2 × 10 4 to 4 × 10 4 km. Figure 2 shows the cospectral power power density for ω q at 0.43 bar for the 2D simulations shown in Figure 1 with along an additional simulation with 48 × 10 4 km, demonstrating that flows with size ∼ 3 × 10 4 km dominate the formation and vertical transport of clouds in the non-rotating 2D systems. Then we perform a horizontal resolution test for the 3D model with a Coriolis parameter f = 4 × 10 −4 s −1 and fixed domain size 30000 km × 30000 km by conducting three simulations with the number of grid points in the horizontal of 100 × 100, 200 × 200, and 300 × 300 (corresponding to horizontal resolutions of 300, 150 and 100 km per grid cell). The resulting horizontal mean total tracer and RMS of tracer deviations as a function of pressure are shown in the left-most panel of Figure 3 for the three simulations. KE power spectra and outgoing thermal flux variance spectra are shown in the middle and right column of Figure 3. The vertical tracer structures are similar among three cases. The power spectra show convergence, especially among cases with 200 and 300 grid points per horizontal axis. This test suggests that the adopted numerical resolution in this study is adequate to resolve major dynamical processes.

APPENDIX B: DIFFERENT FORCING AMPLITUDES
Simulations shown in Figure 5 are somewhat strongly forced such that inertia gravity waves can be energetically important in the power spectra at 0.23 bar when the wavenumber is sufficiently larger than that of the deformation radius. We investigated effects of changing the forcing amplitude via adjusting the abundance of deep condensable vapor -the smaller the condensable vapor, the less forcing. In practice the settling velocity of cloud particles is artificially adjusted smaller; this is because as the circulation weakens due to the smaller forcing, it may not be able to keep cloud particles aloft to sustain the circulation. Figure B1 shows the KE power spectra for different forcing amplitudes. As the forcing becomes weaker, the transition from slope −3 to −5/3 gradually vanishes. In the two least forced cases, their KE power spectra exhibit a slope of −3 from the deformation radius all the way down to the scale at which KE is affected by the numerical dissipation, and no transition to −5/3 occurs. When the forcing is stronger, the KE power spectrum flattens out at relatively high wavenumber. Interestingly, the peak of the KE spectrum systematically moves to smaller wavenumber (larger size) as the forcing increases. This paper has been typeset from a T E X/L A T E X file prepared by the author.