Atmospheric circulation of brown dwarfs and directly imaged exoplanets driven by cloud radiative feedback: global and equatorial dynamics

Brown dwarfs and directly imaged exoplanets exhibit observational evidence for active atmospheric circulation, raising critical questions about mechanisms driving the circulation, its fundamental nature, and time variability. Our previous work demonstrated the crucial role of cloud radiative feedback on driving a vigorous atmospheric circulation using local models that assume a Cartesian geometry and constant Coriolis parameters. In this study, we explore the properties of the global dynamics. We show that, under relatively strong dissipation in the bottom layers of the model, horizontally isotropic vortices are prevalent at mid-to-high latitudes while large-scale zonally propagating waves are dominant at low latitudes near the observable layers. The equatorial waves have both eastward and westward phase speeds, and the eastward components with typical speeds of a few hundred m/s usually dominate the equatorial time variability. Lightcurves of the global simulations show variability with amplitudes from 0.5 percent to a few percent depending on the rotation period and viewing angle. The time evolution of simulated lightcurves is critically affected by the equatorial waves, showing wave beating effects and differences in the lightcurve periodicity to the intrinsic rotation period. The vertical extent of clouds is the largest at the equator and decreases poleward due to the increasing influence of rotation with increasing latitude. Under weaker bottom dissipation, strong and broad zonal jets develop and modify wave propagation and lightcurve variability. Our modeling results help to explain the puzzling time evolution of observed lightcurves, a slightly shorter period of variability in IR than in radio wavelengths, and the viewing angle dependence of variability amplitude and IR colors.


INTRODUCTION
Brown dwarfs (BDs) are objects with masses intermediate between stars and giant planets, but are analogous to giant planets in temperature, composition and size (Burrows et al. 2001). Isolated BDs are easier to observe than exoplanets orbiting bright stars, making them ideal targets to investigate physical, chemical and dynamical processes in the context of planetary atmospheres (Apai et al. 2017). Growing observations of BDs in the past decade have revealed that active atmospheric circulation is common among BDs. Here we summarize several key types of observations that set direct constraints on the large-scale circulation: • Broadband lightcurve time variability at infrared (IR) wavelengths is commonly observed for field BDs (e.g., Gelino et al. 2002;Clarke et al. 2008;Artigau et al. 2009;Gillon et al. 2013;Buenzli et al. 2014;Wilson et al. 2014;Metchev et al. 2015;Leggett et al. 2016;Miles-Páez et al. 2017;Apai et al. 2017; ★ E-mail: xianyu.tan@physics.ox.ac.uk 2017; Vos et al. 2019;Eriksson et al. 2019;Bowler et al. 2020;Hitchcock et al. 2020;Vos et al. 2020, see also recent reviews by Biller 2017 andArtigau 2018). Most of them are likely caused by the rotation of large-scale surface inhomogeneities related to different cloud opacity and/or temperature that move in and out of view as the objects rotate. The shapes of some lightcurves are irregular and change over timescales comparable to or slightly longer than rotational timescales, indicating rapid evolution of the large-scale spatial patchiness in some cases. Apai et al. (2017) summarized a few types of irregularities of the lightcurve variability and suggested a possible solution by invoking differential zonal propagation of waves.
• Multi-wavelength near-IR lightcurve variability provides additional information on the vertical structure of the surface patchiness because different wavelengths probe different atmospheric pressures (Buenzli et al. 2012;Radigan et al. 2012;Biller et al. 2013;Apai et al. 2013;Yang et al. 2015Yang et al. , 2016Lew et al. 2016;Zhou et al. 2018;Zhou et al. 2020;Lew et al. 2020). The amplitude of the variability changes with wavelength, showing smaller amplitude near the water absorption bands and larger amplitude near spectral windows. These are typically consistent with the scenario that the surface patchiness results from a combination of different cloud thicknesses rather than a combination of a single cloud type and a completely clear-sky region (Buenzli et al. 2012;Apai et al. 2013;Yang et al. 2015;Lew et al. 2016Lew et al. , 2020. Alternative scenarios invoking variation of temperature in clear-sky atmospheres have also been proposed (Robinson & Marley 2014;Morley et al. 2014;Tremblin et al. 2020). The phases of lightcurve rotational variations can be different at different wavelengths, and the extent of this difference varies with spectral types (Yang et al. 2016).
• Viewing-angle dependence of variability and near-IR color of field BDs has been suggested by Vos et al. (2017), Vos et al. (2018) and Vos et al. (2020), who showed that BDs viewed more equatoron tend to have higher variability amplitudes and redder near-IR colors than those viewed more pole-on. These results indicate the possible presence of systematic equator-to-pole differences of cloud properties, temperature and/or chemical composition.
• Doppler imaging could provide perhaps the most direct constraints on the presence of large-scale surface patchiness and their morphology. This has been done for Luhman 16B (the spectral type T1 component in the BD binary system Lhuman 16AB) and yielded surface flux differences between the patchiness of a few hundred K in terms of effective temperature (Crossfield et al. 2014).
• Simultaneous tracking of near-IR and radio variability, in which the former measures the period at which atmospheric features rotate in and out of the view and the latter likely reflects the rotation period of the interior, could track possible differential rotation between the atmosphere and the interior. Recently, Allers et al. (2020) applied this technique to a nearby T6 brown dwarf and showed that the period of the near-IR variability is slightly shorter than that of the radio emission, suggesting that the dominant atmospheric features travel eastward relative to the interior with a zonal speed of a few hundred m s −1 .
• Net near-IR polarization could be due to scattering by cloud particles together with either the presence of oblateness due to fast rotation or cloud inhomogeneity (e.g., Sengupta & Marley 2009;Marley & Sengupta 2011). Recently, Millar-Blanchaer et al. (2020) unambiguously detected net near-IR polarization as well as its time variability on Luhman 16B. The time-mean polarization indicates a latitudinally dependent cloud structure, and its time variation indicates longitudinal cloud patchiness on top of the zonal cloud structure.
Directly imaged extrasolar giant planets (EGPs) are mostly young, hot and relatively distant from their host stars such that they receive negligible external stellar irradiation compared to their interior heat flux. Their spectrum and near-IR colors show similarities to the dusty field BDs (e.g., Currie et al. 2011;Barman et al. 2011;Marley et al. 2012;Rajan et al. 2015;Chauvin et al. 2017;Stolker et al. 2020). Near-IR lightcurve variability has been observed on directly imaged EGPs and planetary-mass, free-floating objects 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). From a meteorological point of view, directly imaged EGPs resemble lowgravity versions of BDs and fall into the same category as field BDs in terms of atmospheric dynamical regime.
These observations motivate the investigation of global atmospheric dynamics of BDs and directly imaged EGPs, their fundamental properties and effects on cloud formation and chemistry. There have been several studies investigating atmospheric dynamics appropriate for these objects (Freytag et al. 2010;Zhang & Showman 2014;Tan & Showman 2017;Showman et al. 2019;, and see recent reviews by Showman et al. 2020 andZhang 2020). Cloud radiative feedback has been proposed as a robust and novel mechanism to drive spontaneous atmospheric variability and dynamics in these atmospheres . This mechanism works as follows: imagine two adjacent atmospheric columns, one of which has a thick cloud layer and the other has a relatively thin cloud layer. The two columns experience different radiative heating and cooling rates due to the different cloud opacity, which then generate isobaric temperature differences. This drives an overturning flow that maintains the patchy clouds against settling, potentially sustaining the circulation. In a previous work, we confirmed the viability of this mechanism and explored dynamical properties of the circulation using local three-dimensional models that assume a Cartesian geometry and a constant Coriolis parameter across the domain . We demonstrated the importance of rotation in regulating the large-scale atmospheric dynamics and cloud formation, showing that the typical horizontal length scales of vortices and cloud patterns are inversely proportional to the Coriolis parameter if other parameters are held fixed, and that the mean cloud vertical extent decreases with decreasing rotation period.
The latitudinal variation of the Coriolis parameter in global geometry, the so-called effect where = / and is distance increasing northward, introduces additional dynamical behaviors compared to those in the -plane models . For example, the effect can introduce horizontal anisotropy in the turbulence (Rhines 1975;Vallis & Maltrud 1993) as well as large-scale atmospheric waves (Holton & Hakim 2012), both of which play crucial roles in the global circulation of giant planets (see reviews by Ingersoll et al. 2004;Vasavada & Showman 2005;Del Genio et al. 2009;Showman et al. 2018). In addition, the regional models in  occupy only a limited fraction of the surface area of BDs and EGPs, and therefore qualitative comparisons between model outputs and observed lightcurve variability were lacking.
In this study, we extend the modeling framework of  to a global geometry to investigate the global atmospheric circulation driven by the cloud radiative feedback and the resulting lightcurve variability. We show that when the bottom frictional dissipation is relatively strong, horizontally isotropic storms and turbulence are prevalent at mid-to-high latitudes while zonally propagating waves are present at low latitudes near the observable layer. The differential propagation of equatorially trapped waves induces short-term evolution of simulated lightcurves, analogous to the wave beating effects described in Apai et al. (2017). Eastward propagating equatorial Kelvin waves are sometimes dominant, causing a slightly shorter rotation period of the atmospheric features relative to the underlying planetary rotation period, which agrees well with observational results by Allers et al. (2020). When the bottom dissipation is weak, strong and broad zonal jets develop and modify wave propagation and lightcurve variability. We find systematic equator-to-pole differences of clouds and temperature structures due to the latitudinal variation of the Coriolis parameter , supporting recent observations by Vos et al. (2017) and Vos et al. (2020). Largescale equatorial disturbances may help to explain the nature of longitudinal variation in Doppler mapping (Crossfield et al. 2014) and time-varying polarization (Millar-Blanchaer et al. 2020) of Luhman 16B.
The paper is organized as follows. Section 2 briefly introduces the global numerical model. Section 3 describes results of models with different rotation periods and drag timescales as well as their resulting lightcurve variability. In Section 4 we discuss implications of the results and draw conclusions.

MODEL
The general circulation model (GCM) used in this study is fully described in , and here we apply the GCM to a global domain. Key elements of the model are summarized below. We solve the standard 3D hydrostatic primitive equations using an atmospheric GCM, the MITgcm (Adcroft et al. 2004). Two tracer equations which represent the mass mixing ratio of condensable vapor and clouds are integrated simultaneously. We assume the ideal gas law for the equation of state. The deep layers in our models have reached the convective zone which is typically at pressures larger than a few bars for typical L and T dwarfs (Burrows et al. 2001). Effects of rapid convective mixing on both entropy and tracers are parameterized using a simple convective adjustment scheme that instantaneously homogenizes entropy and tracers in the convectively unstable regions.
A Rayleigh drag is applied to the horizontal winds at pressures higher than 5 bars to crudely represent interactions between the modeled atmosphere and the quiescent interior where large-scale flows are likely retarded due to significant magnetohydrodynamic dissipation. Rapid convective mixing of specific entropy and fast rotation may lead to the Taylor-Proudman effect that could be efficient in slowing down large-scale winds in the shallow outer layer (see the detailed argument in Schneider & Liu 2009). Because of the much hotter interior of L and T dwarfs and likely strong magnetic fields of BDs (∼kG, Kao et al. 2016Kao et al. , 2018, the "quiescent" region inside BDs is expected to extend to a much shallower depth than that of Jupiter. Therefore the drag is applied uniformly in the horizontal direction. The form and strength of the drag are highly uncertain because of the unknown nature of interactions between the interior and the shallow outer layer. Yet no theoretical study systematically investigate how such interactions should be parameterized in GCMs of gaseous planets. Nevertheless, this simple drag serves to dissipate flows in our bottom model layers, and similar drag schemes have been used in previous studies of Jupiter (Schneider & Liu 2009) and hot Jupiters (e.g., Liu & Showman 2013;Tan & Komacek 2019;Carone et al. 2020).
The drag is the strongest at the bottom boundary and is characterized by a drag timescale drag , then it linearly decreases with decreasing pressure until it reaches zero at 5 bars. There is no drag at pressures lower than 5 bars. As in , we set a relatively strong drag timescale drag = 10 5 s for the major suit of models. We have tested models with stronger drags ( drag = 10 4 and 10 3 s), and they are quantitatively similar to those with drag = 10 5 s. Therefore, drag = 10 5 s is chosen to represent dynamics in the "strong-drag" regime. For drag sufficiently larger than 10 5 s, dynamics is in the weak-drag regime. The drag is said to be "strong" in a practical sense that the resulting zonal-mean zonal flows near the bottom layer are rather weak (with velocities much smaller than 100 m s −1 ) based on our modeling results discussed in Section 3.2, whereas the drag is said to be "weak" when the resulting zonal flows near the bottom layer are strong (with velocities much greater than 100 m s −1 ). Using scaling argument in the convective layers, Showman & Kaspi (2013) argue that the characteristic speeds of large-scale zonal flows near the top of convective zone range from only a few to tens of m s −1 depending on the number of jets. We therefore expect that the strong-drag regime might be appropriate for BDs and hot directly imaged EGPs. Nevertheless, in section 3.2, we show results with weaker drags of drag = 10 6 and 10 7 s to explore possible circulation patterns for two reasons. First, it is theoretically motivated to understand dynamics in different regimes in spite of its applicability in realistic situations. Second, recent gravity mea-surements of Jupiter and Saturn suggest that significant meridional density gradients (and therefore the vertical wind shears associated with the zonal jets) exist deep in the convective zone of Jupiter and Saturn Guillot et al. 2018;Iess et al. 2019). This indicates that strong organized zonal flows near the top of convective zone remain possible for BDs.
Atmospheric motions are driven by the horizontal pressure gradient that is rooted from horizontal differential radiative heating and cooling. In calculating the radiative flux, we assume a gray atmosphere with a single broad thermal band for simplicity and computational efficiency. The gas opacity in our model atmosphere is gas = max( R,g , min ), where R,g is the Rosseland-mean opacity from Freedman et al. (2014) assuming solar composition. The Rosseland-mean opacity gives a good estimation of radiative flux in the optically thick limit. The atmosphere above about 1 bar is optically thin by the gas opacity alone, and there is no good choice a priori for the opacity in the grey approximation. A minimal opacity min is then imposed in the gas opacity. We assume min = 10 −3 m 2 kg −1 , which is consistent with that used in semi-grey models in Guillot (2010) for hot giant planets. The total atmospheric opacity is simply the sum of the gas and cloud opacities = gas + , the latter of which is determined by the cloud mixing ratio as a function of time and location.
The cloud formation scheme is the same as that in . Cloud forms when the mixing ratio of condensable vapor exceeds the prescribed saturation mixing ratio . Otherwise, evaporation occurs when the vapor mixing ratio is less than . The saturation mixing ratio is assumed to depend on pressure alone: is deep at a condensation pressure cond and then rapidly decreases when pressure is less than cond . At pressures larger than cond , is assumed to be arbitrarily large such that no condensation would occur. The condensation pressure cond is assumed to be 0.5 bar, roughly consistent with conditions of typical mid-L dwarfs in which clouds form at and above ∼ 1 bar (e.g., Tsuji 2002;Burrows et al. 2006). The deep mixing ratio deep is assumed to be 2 × 10 −4 kg/kg. At pressures higher than 5 bars, the vapor field is relaxed towards deep over a timescale of 10 3 s. The assumed deep in this study is less than that of silicate vapor in atmospheres with solar abundance. This value is tuned to generate cloud opacity that greatly exceeds the gas opacity, but not to trigger convection within the cloud-forming region. Convection within the cloud forming region complicates the dynamics as it would introduce instability that does not totally rely on resolved large-scale flow . We aim at exploring large-scale cloud-driven dynamics in a cleanest possible context, therefore we do not chose a realistic deep .
We assume a constant cloud particle number per dry air mass N (in unit of kg −1 ) throughout the atmosphere. Cloud particles are assumed to be spherical and have a single size locally in each grid cell, and the particle size is then determined via a given N and timeand location-dependent amount of condensate (using the scheme of Tan & Showman 2020). We assume N = 10 11 kg −1 , which results in typical cloud particle size around 0.5 m given the assumed deep , consistent with that expected for L dwarfs (Saumon & Marley 2008;Burningham et al. 2017). Optical properties of cloud particles, including the extinction coefficient, scattering coefficient and asymmetry parameter are averaged using the Rosseland-mean definition. Tables containing these parameters as functions of temperature and pressure are pre-calculated for the GCM to read in during the integration. We use enstatite (MgSiO 3 ) to represent properties of the cloud species, including a density = 3190 kg m −3 and the refractive index. Enstatite is one of the most representative cloud species in atmospheres of L and early T dwarfs, and is appropriate for our atmospheric conditions (Ackerman & Marley 2001;Saumon & Marley 2008;Marley & Robinson 2015). However, the choice of cloud species is not critical in this study as long as the cloud opacity exceeds the gas opacity.
The global models in this study have different horizontal resolution depending on the rotation period-the shorter the rotation period, the smaller the deformation radius and thus higher resolution is needed. Our horizontal resolution of the cubed-sphere grid is C96 (equivalent to 384×192 grid points in longitude and latitude), C128 (512×256), C192 (768×384) and C256 (1024×512) for rotation period of 20, 10, 5 and 2.5 hours, respectively. These horizontal resolutions are among the highest published so far for exoplanet and brown dwarf global models 1 . The radius of the object is assumed to be 7 × 10 7 m, similar to the Jovian radius . Although radius of real BDs varies from more than 2 to slightly less than depending on their masses and ages, our choice of radius should be representative for field BDs. A standard fourth-order Shapiro filter is applied in the horizontal velocity and temperature fields to maintain numerical stability (Shapiro 1970). The pressure domain is between 10 bars and 10 −3 bar, which is discretized evenly into 55 layers in log pressure. The model pressure domain is deep enough to reach the convective zone, where we argue that rapid convective mixing leads to small horizontal temperature gradient. The upper boundary is high enough for models to safely capture dynamics associated with clouds. The temperature at the bottom boundary (10 bars) is fixed at 2600 K, resulting in atmospheric temperatures comparable to those of L dwarfs. We adopt physical parameters relevant for BDs, including the specific heat = 1.3 × 10 4 Jkg −1 K −1 , the specific gas constant = 3714 Jkg −1 K −1 , and a surface gravity = 1000 ms −2 . Models with drag = 10 5 s equilibrated after ∼ 100 simulation days. After equilibration, we integrated them for additional 300 days, and the time-mean statistical results were obtained with outputs of the last 200 days. Models with drag = 10 6 and 10 7 s equilibrated after ∼ 200 and ∼ 1000 simulation days, respectively, and their statistical results were obtained similarly to those with drag = 10 5 s.

Temperature and cloud patterns
Results of global models with rotation periods of 20, 10, 5 and 2.5 hours and with a drag timescale drag = 10 5 s (note that the frictional drag is applied only at pressures greater than 5 bars) are shown in Figure 1. Figure 2 shows the corresponding outgoing top-of-atmosphere thermal flux. These results are obtained after the models reach statistical equilibrium. Pressures of thermally radiating levels are sensitive to cloud top pressures-the radiating pressure could be lower than 0.05 bar in cloudy regions and could be about 3 bars in clear-sky regions. Dynamics at 0.23 bar is representative for layers near and slightly below the cloud forming region. Quantities at slightly higher or lower pressures are qualitatively very similar to those at 0.23 bar.
There is a lack of strong zonal jets, in the sense that the zonalmean zonal winds are much weaker than the eddy velocities, in all models shown in Figure 1 and 2 due to dissipation of kinetic energy by the strong bottom drag. The circulation in the cloud forming 1 If the horizontal resolution is insufficient to fully resolve dynamics within the deformation radius of local regions (especially near the polar regions), storms will cease to exist and there would not be dynamics in these regions. region is dominated by active vortices, turbulence and waves. At 0.23 bar, typical isobaric temperature differences are above 100 K, and local horizontal wind speeds can exceed 1000 m s −1 . Mid-tohigh latitudes are filled with cyclones 2 which are less cloudy than their surroundings, and anticyclones which are associated with thick clouds. In the cloudy areas, less thermal radiation escapes to space due to large cloud opacity, and this cloud greenhouse effect warms up the lower cloud forming region. In the cloud-free areas, more radiation can escape to space from the hotter, deeper region due to the lack of cloud opacity, efficiently cooling down the atmospheric column. This patchy-cloud configuration is self-maintained by the circulation driven by clouds themselves . Positive buoyancy is generated in the cloudy regions, and it maintains the cloudiness against cloud gravitational settling by transport of condensable vapor upward to the condensation level; whereas negative buoyancy in the cloud-free regions generates downwelling that advects cloud-free air from above, maintaining the IR cooling. In rapidly rotating conditions, the tendency of geostrophy implies that the warm, cloudy regions are anticyclonic and cool, cloud-free regions are cyclonic. The spatial patterns of the outgoing top-ofatmosphere thermal flux are highly correlated with that of cloud patchiness. Other than the vortices, turbulence and transient waves with smaller horizontal length scales and higher oscillation frequencies are also present. The basic findings from these global models agree well with those of box models with a fixed Coriolis parameter across the domain . At a given rotation period, typical sizes of vortices are generally larger at lower latitudes and smaller at higher latitudes. This is more prominent in rapidly rotating models with rotation periods of 5 and 2.5 hours. The overall sizes of vortices are larger with longer rotation period. This is because the typical sizes of vortices driven by cloud radiative feedback are close to the Rossby deformation radius = /| | , where is the phase speed of gravity waves, = 2Ω sin , Ω is planetary rotation rate and is latitude. The dependence of on Ω and naturally leads to the spatial and rotation-period dependence of vortex sizes. Individual vortices are short lived and experiences merging, disruption, shearing and straining over typical timescales of tens of hours. There is no systematic migration of individual vortices after they form. This is likely because the atmosphere is filled with vortices, and individual vortex is significantly disrupted by adjacent vortices before they can migrate over a long distance.
As shown by , greater | | leads to overall thinner clouds. Qualitatively, stronger rotation leads to a larger degree of geostrophic balance in the flows which exerts greater suppression to the vertical velocity, making the flows less efficient to vertically transport tracers against cloud gravitational settling. Therefore, the variation of the Coriolis parameter from the equator to the poles as well as the different rotation period have significant consequences for the vertical extent of clouds. In Figure 3 we show the time-and zonal-mean cloud mixing ratio as a function of pressure and latitude for models with four rotation periods. With a given rotation period, the vertical extent of clouds is largest at the equator and decreases poleward. The cloud thicknesses at the equator are almost the same among models with different rotation periods. At the poles, the zonal-mean cloud thickness is smaller when the rotation period is shorter. For rapidly rotating models, spans a larger range and is responsible for the greater change of the zonal-mean cloud thick-  ness. As a result of the equator-to-pole cloud thickness variation, the time-and zonal-mean temperature-pressure (T-P) structures exhibit a systematic equator-to-pole variation due to the cloud radiative feedback. Figure 4 shows the time-averaged zonal-mean T-P profiles at selected latitudes for the model with a rotation period of 2.5 hours. The isobaric temperature variation can reach more than 100 K from the equator to the poles. On top of this systematic equator-to-pole zonal-mean variation, instantaneous differences in the T-P structures associated with the cloudy and clear-sky regions are also present, similar to that shown in Tan & Showman (2020).

Equatorial wave properties
Dynamics near the equator shows qualitative differences to that at mid-to-high latitudes, which can be seen from two features. The first is the morphology of the eddies. Eddies near the equator are likely more elongated in the zonal direction, which is in contrast to the horizontally isotropic vortices at mid-to-high latitudes. This feature is more prominent in models with rotation periods of 5 and 2.5 hours. The second difference is the time evolution of the eddies. Equatorial eddies exhibit systematic eastward or westward propagation, while trajectories of mid-to-high-latitude vortices are more akin to random  walks and show no systematic migration along certain directions. Figure 5 shows the Hovm oller diagrams (longitude-time sequence) of outgoing thermal flux at the equator as a function of longitude and time for models with different rotation periods. All models exhibit obvious characteristic eastward propagating patterns (those moving towards the upper right of the panels) with zonal speeds on the order of a few hundred m s −1 . On the contrary, Figure A1 in Appendix A similarly shows the Hovm oller diagrams at 45 • latitude where stochastically evolving vortices dominate, and there is not evidence of eastward or westward propagation of the eddies. These propagating eddies at low latitudes have much larger zonal wavelengths than the mid-to-high-latitude vortices (comparing Figures 5 and A1), showing dominant features characterized by zonal wavenumbers of a few. Westward propagation of eddies is also present with a slower phase speed. The time evolution of the zonally propagating eddies is not always coherent, and often shows complications as they evolve. For example, some existing perturbations disappear when they propagate while some new perturbations are generated. The propagating eddies at low latitudes are likely equatorially trapped waves, and they exert significant effects on the short-term evolution of lightcurve variability as will be discussed in Section 3.3. Based on both the morphology and time evolution of the eddies, latitudinal boundaries that separate the two groups of atmospheric motions seem to emerge. The boundary is closer to the equator with shorter rotation period. This boundary is likely related to the equatorial deformation radius eq = √︁ / . Poleward of this distance from the equator, there is sufficient room for mature vortices to develop and the dynamics is dominated by vortices. Within this scale from the equator, dynamics are shaped by equatorially trapped waves that are coupled to the cloud radiative effect. Using ∼ 2000 m s −1 obtained from fitting theoretical off-equatorial deformation radius to the measured sizes of dominant vortices as a function of in the constant− models of , the equatorial deformation radius eq extends to about 24 • , 17 • , 12 • and 9 • away from the equator for models with rotation periods of 20, 10, 5 and 2.5 hours, respectively. This is roughly consistent with that shown in Figure 1 and 2 and their time evolution (not shown). Now we turn to characterize properties of equatorial waves that are present in our simulations. We perform a spectral analysis at the equator in the wavenumber-frequency domain, similar to that performed in Wheeler & Kiladis (1999) and Showman et al. (2019). The brief procedure is the following. We perform two-dimensional fast Fourier transforms on the outgoing thermal flux as a function of longitude and time near the equator to obtain the raw Fourier coefficients in the wavenumber-frequency space. These raw coefficients are then heavily smoothed in the wavenumber-frequency space to generate the background spectrum. Finally, the raw coefficients are divided by the background spectrum to obtain the relative power spectrum, in which significant signals will appear to have values greater than 1. The eddy fields are further decomposed into symmetric and antisymmetric components about the equator because it helps to clarify the wave properties (Wheeler & Kiladis 1999;Kiladis et al. 2009). The relative power spectrum of models with different rotation periods are shown in Figure 6 for the symmetric components and Figure  7 for the anti-symmetric components. Our space-time spectral analysis demonstrates the robust existence of groups of zonally propagating waves. In the symmetric components, a group of eastward waves with zonal wavenumbers between 2 to 10 are present in all models. Weaker signals of westward waves with zonal wavenumbers between −10 and −2 exist in models with rotation periods of 5 and 2.5 hours but not in those with rotation periods of 20 and 10 hours. In the anti-symmetric components, models with rotation periods of 20 and 10 hours exhibit evidence of both eastward and westward waves with zonal wavenumbers in between about −5 to 5, and those with 5 and 2.5 hours host westward waves with zonal wavenumbers between −10 to 1.
We compare the signals in Figure 6 and 7 to dispersion relations of adiabatic, freely propagating equatorial waves derived from the shallow-water system (Matsuno 1966). These waves are characterized by an equivalent depth ℎ , such that the gravity wave speed (which is the same as the Kelvin wave speed) = √︁ ℎ where is the surface gravity. 3 Although waves in our simulations are highly diabatic and tightly coupled to radiative effects of clouds, this serves as a baseline comparison. Solid lines in Figure 6 represent the dispersion relations of adiabatic free equatorial Rossby waves with a meridional mode number = 1 and with three equivalent depths of 300, 150 and 50 m (equivalent to gravity wave phase speeds of 548, 387 and 224 m s −1 ). Dashed lines in Figure 6 represent dispersion relations of Kelvin waves with the same equivalent depths. In the anti-symmetric components shown in Figure 7, = 0 westward mixed Rossby-gravity (MRG) waves and = 0 and eastward inertia gravity (EIG) waves with three equivalent depths of 400, 200 and 100 m (gravity wave phase speed of 632, 447 and 316 m s −1 ) are plotted as dotted lines.
In the symmetric components ( Figure 6), evidence of Kelvin waves are quite strong for all models as the spectral powers follow the Kelvin wave dispersion relations (dashed lines) reasonably well. Amplitudes of the Kelvin waves are stronger among low zonal wavenumbers between 2 − 6. The wave speeds of our simulated Kelvin waves are somewhat dispersive. Kelvin modes with lower zonal wavenumbers (longer zonal wavelengths) typically have slower phase speeds than those with higher zonal wavenumbers. A phase speed of ∼ 550 m s −1 brackets the upper limit of phase speeds in all models. The phase speeds of waves with a zonal wavenumber of 2 in cases with rotation period of 10 and 5 hours are much slower than 220 m s −1 . There is no evidence of equatorial Rossby waves in models with rotation periods of 20 and 10 hours, and the evidence is tentative in models with 5 and 2.5 hours. However, even in the case with 2.5 hours, the westward branch with relatively high frequencies > 0.5 cycle per day (CPD) obviously deviates from the theoretical Rossby-wave dispersion relation. There are no wave signals at frequencies higher than those shown in Figure 6 in the symmetric components, suggesting the absence of high-frequency inertia gravity waves.
In the anti-symmetric components (Figure 7), MRG and = 0 EIG waves are evident in cases with rotation periods of 20, 10 and 5 hours. They have either eastward or westward phase velocities with speeds larger than that of the Kelvin waves. In the case of 2.5 hours, although the wave signals lie in between the dispersion relations of theoretical MRG waves assuming different equilibrium depth, the signals show somewhat constant frequency between zonal wavenumber -7 to -3, making the interpretation less obvious. Nevertheless, these diagrams ( Figure 6 and 7) quantify the propagation of waves and help to recognize the wave types.
Because eastward Kelvin waves seem dominant in both the equatorial time evolution of out-going thermal flux ( Figure 5) and the spectral analysis (Figure 6), we further illustrate their properties in the physical domain by extracting from the total dataset through filtering in the wavenumber-frequency domain, similar to the procedure in Wheeler & Kiladis (1999). At a given latitude, we choose a signal with a zonal wavenumber and a frequency in the wavenumberfrequency domain, then project it back to the longitude-time domain. We repeat this process for all latitudes near the equator and obtain the propagating wave in the physical space. Panel (a) in Figure 8 shows the snapshot of an eastward symmetric wave with a zonal wavenumber of 4 and a frequency of about 0.15 CPD that shows a strong signal in the model with a rotation period of 10 hours (the upper right panel in Figure 6). An example of adiabatic free Kelvin wave is shown in panel (b) as a comparison to our simulated wave. The adiabatic free Kelvin wave has eddy zonal velocity in phase with the pressure anomaly, and the disturbances decay away from the equator. Our simulated wave shows certain similarities to the adiabatic free Kelvin wave -velocities have much larger zonal components than the meridional components in most places, and the geopotential anomalies decay away from the equator. Note that cloud abundances correlate well with the convergence/divergence of the wind field (see a comparison between vectors and the colors that represent the cloud abundance).
Meanwhile, in the simulated wave, the eastward zonal velocity pattern exhibits a moderate eastward phase shift compared to the geopotential anomalies in the zonal direction (whereas those in the adiabatic free Kelvin wave are in-phase). Moreover, the simulated wave exhibits a northwest-southeast phase tilt in the northern hemisphere and a southwest-northeast phase tilt in the southern hemisphere for the geopotential anomalies. This likely reflects a tendency that the circulation resembles wave patterns excited by a stationary heat source that is symmetric about the equator (Matsuno 1966;Gill 1980). An example of such a stationary wave solution subjected to damping on both velocity and layer thickness in the shallowwater system adopted from Showman & Polvani (2011) is shown in panel (c). This so-called Matsuno-Gill pattern is characterized by an eastward shift of the Kelvin component at the equator and a westward shift of Rossby components off the equator. Therefore it shows the northwest-southeast phase tilt in the northern hemisphere and a southwest-northeast phase tilt in the southern hemisphere, as well as that the maximum positive zonal velocity is east of the height maximum at the equator. In our simulations, when cloud forms by convergence at the equator, the accompanied heat source drives the flow towards the Matsuno-Gill shape. However, because clouds themselves are also advected by flows over a timescale that is comparable to that required to form the Matsuno-Gill pattern, the advection of the heat source disrupts the complete formation of the Matsuno-Gill pattern. This is why our simulated wave still shows quantitative discrepancies to that in panel (c).
It is likely that our simulated eastward, Kelvin-like waves resemble properties of both the adiabatic free Kelvin waves and the forced-damped Matsuno-Gill circulation. Recently, using a quasilinear moist shallow-water system, Vallis & Penn (2020) showed that the Matsuno-Gill-like circulation is excited by latent heating (which itself is coupled to the flow, somewhat analogous to the cloud radiative heating in our cases) at the equator, and a slow eastward migration of the whole pattern is driven by the interactions of the flow and moisture field. It is unclear how much the mechanism shown in Vallis & Penn (2020) participates in the eastward propagation of our simulated waves in addition to that related to the free Kelvin wave. Teasing out the detailed mechanisms in our simulations is beyond the scope of this paper as they are highly nonlinear. Here we still refer these zonally traveling disturbances as waves, partly because these disturbances follow wave dispersion relations reasonably well (as shown in Figures 6 and 7).
The origin of the equatorial waves in our simulations is likely related to self-excitation by cloud radiative feedbacks. However, even treating them as waves, mechanisms controlling the wave properties are unclear and little previous work exists. Linear stability analyses on the cloud radiatively coupled dynamics was carried out for inertia gravity waves in a -plane system and a quasi-geostrophic system by Gierasch et al. (1973), and they showed that unstable modes are possible as a result of cloud radiative instability. Here, we extend their theory to the equatorial -plane which is appropriate for equatorially trapped waves. The derivations are shown in Appendix B. In the linear theory, we find that unstable modes are possible for a set of Kelvin and = 0 MRG modes, suggesting a source of kinetic energy on the equatorial waves shown in our models. However, the unstable modes from the theory do not propagate. All other modes that resemble the classic adiabatic free equatorial waves discovered by Matsuno (1966), including propagating Kelvin, Rossby, MRG, eastward and westward inertia gravity modes, show damping on the eddy amplitudes due to thermal radiation. The linear theory is not a total failure because it still predicts the kinetic energy sources for most equatorial waves in our simulations, but it cannot explain the propagation of the simulated waves. This situation is similar to that in the non-rotating two-dimensional system shown in .
The wave speeds shown in Figure 6 and 7 are significantly slower than expected from adiabatic waves with long vertical wavelengths in conditions appropriate for our simulated atmospheres, the latter of which is characterized by a much larger dry equivalent depth ℎ ,dry = 2 ,dry / ∼ 4000 m, where ,dry ≈ 2 ∼ 2000 m s −1 , is the Brunt-Vaisala frequency and is the scale height. This is similar to equatorial waves in Earth's troposphere which are affected by latent heat released from moist convection (Wheeler & Kiladis 1999). To understand the reduced phase speed of tropical waves by diabatic effects, the following idealized framework has been proposed (see a review by Kiladis et al. 2009). Suppose that the large-scale diabatic heating and cooling is included in the linearized thermodynamics system as the following where is the geopotential (note that / is proportional to the temperature perturbation), = − log( / ) is the log-pressure coordinate, is pressure, is a reference pressure, and is vertical velocity in log-pressure coordinates. Furthermore, if the heating and cooling are proportional to the vertical velocity such that = 2 (where is an arbitrary constant), then Equation (1) becomes If is positive and less than 1, diabatic heating and cooling effectively reduce the stability of the atmosphere and therefore reduce the wave phase speeds. For Earth's tropical waves, the problem reduces to determine theoretical values based on parameterized moist convection for to explain the observed wave speeds (e.g., Neelin & Held 1987;Emanuel et al. 1994;Haertel & Kiladis 2004), but the problem has not been completely solved (Kiladis et al. 2009). For waves coupled purely with cloud radiative feedback, the relation = 2 likely holds as well. In regions with upwelling, vapor is advected above the condensation level and clouds form, which generates warming near the condensation level due to cloud radiative effect. In regions with downwelling, cloud-free air is advected downward and clears out the region, enhancing the IR flux to space and inducing cooling. To quantitatively confirm the positive correlation between upwelling and heating, or downwelling and cooling, we calculate the cospectral power density for the quantity in the cloud forming layer, which is simply 2R( * ) where and are the coefficients at wavenumber space for vertical velocity and heating rate, and * is the conjugate of . A detailed description of similar exercises in constant− models is referred to . The cospectral power density shows positive values for almost all zonal wavenumbers, proving that = 2 with a positive holds well in our models. Therefore, we expect the above idealized framework may be used to qualitatively understand the reduced wave speeds seen in our simulations (however it does not address questions as why and how certain wave modes are excited but others are not). As for the detailed determination of , we defer to a future study.

Results with weaker bottom drag
We examine properties of the circulation when the bottom frictional drag becomes weaker by performing two additional models with drag timescales of drag = 10 6 s and 10 7 s and a rotation period of 5 hours. Note that in all models, the frictional drag is applied only at pressures greater than 5 bars, well below the cloud condensation level of 0.5 bar. Figure 9 shows the time-averaged zonal-mean zonal wind as a function of pressure and latitude for models with drag timescales of drag = 10 5 , 10 6 and 10 7 s from the top to the bottom panel. When the drag is strong ( drag = 10 5 s), the zonal-mean zonal wind is much weaker than local horizontal wind speeds of vortices and  geopotential anomalies from 0.5 × 10 4 to 3.5 × 10 4 m 2 s −2 and thick dashed lines are negative geopotential anomalies from −3.5 × 10 4 to −0.5 × 10 4 m 2 s −2 (with H marking the highgeopotential center and L marking low-geopotential center). Arrows represent wind vectors, and colours represent cloud mixing ratio associated with this wave. Panel (b): an example of an adiabatic free Kelvin wave as a function of non-dimensional zonal and meridional distances. Similar to panel (a), arrows represent wind vectors and thick lines represent height perturbations in the shallow-water system. Panel (c): An analytic forced-damped stationary wave solution forced by a zonal-wavenumber-1 stationary forcing centered around (0,0). Light colours represents larger values in the height field and arrows are eddy wind vectors. This is adopted from Showman & Polvani (2011), which the details of the solution is referred to.
turbulence. The zonal-mean zonal flow at low latitudes exhibits an interesting vertical wind shear, with a westward mean flow centering at around 0.5 bar and an eastward mean flow around 0.01 bar. The equatorial eastward jet centering around 0.01 bar corresponds to a superrotation, which requires up-gradient angular momentum transport to the equator by eddies. Our analysis (not shown) suggests that horizontal eddy momentum transport by transient waves is responsible for this local superrotation. This might be somewhat analog to those proposed for superrotation in solar-system bodies, such as Venus, Titan, tropospheres of Jupiter and Saturn (see a recent review by Imamura et al. 2020). In the model with drag = 10 6 s, a broad equatorial westward jet and high-latitude eastward jets with speeds ∼ 400 m s −1 emerge. At low latitudes, the vertical wind shear is similar to that with drag = 10 5 s despite the overall equatorial jet velocity being westward. The jet speeds are comparable to the horizontal mean eddy velocity, but the vortices are still nearly isotropic at mid-to-high latitudes. Near the equator, the equatorial waves are Doppler shifted by the westward jet whose speed is comparable to that of wave propagation. In the model with drag = 10 7 s, the jet speeds become much stronger, reaching about −2200 m s −1 at the equator and more than 1500 m s −1 at about ±50 • latitude. Compared to the jet structure with drag = 10 6 s, the meridional width of the equatorial jet with drag = 10 7 s is very similar and only the jet speed increases. However, at high latitudes, the cores of eastward jets are closer to the equator than those with drag = 10 6 s. There is a strong barotropic (pressure-independent) component for zonal jets with drag = 10 7 s.
The RMS horizontal eddy velocity in the cloud forming regions is about 400 − 550 m s −1 , which is much larger than the jet speeds in the model with drag = 10 5 s. Even in the model with drag = 10 6 s, the jet speeds are only comparable to the eddy velocities. In the following analysis, we focus primarily on the model with drag = 10 7 s in which the jet speeds well exceed the eddy velocities and the dynamics exhibits obvious differences to those with drag = 10 5 s. Cloud formation is affected by the presence of strong zonal jets in the case with drag = 10 7 s. Wind shears at the flanks of the midlatitude jets create strong cyclonic regions poleward of the jets and strong anti-cyclonic regions equatorward of the jets. Anti-cyclones which are associated with cloud formation tend to be more vulnerable in the strong cyclonic zones, in the sense that they are sheared apart and destroyed more easily than in other regions. On the other hand, formation of cyclones which are associated with thin clouds are more prevalent in the cyclonic zones. In addition, vortices tend to migrate towards regions with background absolute vorticity closer to that of vortices (e.g., Scott 2011;O'Neill et al. 2015). There is a tendency that cyclones formed near the mid-latitude jets migrate to the poleward flanks of the jets where it is cyclonic, whereas anticyclones tend to migrate to the equatorward flanks of the jets. These behaviors are observed in time evolution of storms in the simulation with drag = 10 7 s (not shown). All these imply that, clouds associated with anti-cyclones are relatively depleted in the polar flanks of the mid-latitude jets and are more abundant in the equatorial flanks of the jets. The vortex behaviors influenced by the strong zonal jets described above may be one of the mechanisms responsible for the cloud mixing ratio at 0.23 bar in panel (a) of Figure 10. There are also zonal-mean upwelling equatorward of the mid-latitude jets where it is cloudy and zonal-mean downwelling poleward of the midlatitude jets where it is less cloudy (not shown). This mean meridional circulation may or may not be a mechanism driving the meridional cloud gradient. It could simply be radiatively driven given the existing meridional cloud gradient. Shapes of vortices are deformed due to the significant horizontal wind shear, especially near latitudes between 15 • −45 • and 50 • −70 • (see Figure 10). Interestingly, there is a strong polar cyclone at each pole with the center of the cyclone not aligned with the pole. Clouds are also strongly depleted in the cyclones, and the same reasons mentioned above for cloud depletion poleward of the mid-latitude jets might also be relevant. An instantaneous temperature map at 0.23 bar is shown in panel (b) and instantaneous absolute vorticity + is shown in panel (c) of Figure 10, where = k · ∇ × v is the relative vorticity, v is horizontal velocity vector, k is the local upward unit vector on the sphere and ∇ is the horizontal gradient in pressure coordinates. Due to the radiative effects of clouds, strong cyclonic regions where clouds are relatively depleted are systematically colder than anti-cyclonic regions.
The development of robust zonal jets with a weaker bottom frictional drag is consistent with that found in Showman et al. (2019), who studied atmospheric circulation of brown dwarfs using horizontally isotropic, randomly evolving thermal forcing. Storms and vortices excite Rossby waves, and their generation, propagation and wave breaking interact with the mean flow, driving zonal jets (Dritschel & McIntyre 2008). Without efficient removal of kinetic energy by the strong bottom frictional drag, zonal jets can be pumped up and maintained by wave-mean-flow interactions. Some features of the jets in our simulations are interesting. First, the jets are quite meridionally broad, with the equatorial westward jet extending to ±40 • latitude and subsequent eastward jets at about ±70 • latitude for the case with drag = 10 6 s and at about ±48 • latitude for the case with drag = 10 7 s. As a comparison, Jupiter has about 7 subtropical zonal jets in each hemisphere (e.g., Ingersoll et al. 2004). Classic turbulence-driven jet theory predicts that the meridional jet spacing is related to wind speed via the Rhines scale √︁ / where is a characteristic wind speed, indicating that the number of jet on a sphere is roughly given by jet ∼ √︁ 2Ω / where Ω is the rotation rate and is planetary radius (see reviews by, e.g., Vasavada & Showman 2005;Dritschel & McIntyre 2008;Showman et al. 2010). Given a typical wind speed of several hundreds of m s −1 in the case with drag = 10 6 s and ∼ 2000 m s −1 in the case with drag = 10 7 s, one would expect a number of zonal jets of ∼ 9 in the case with drag = 10 6 s and ∼ 5 in the case with drag = 10 7 s, which are obviously greater than that seen in our simulations. Second, the equatorial jets are westward in our simulations, whereas those in Jupiter and Saturn are eastward. The broader-than-expected jets are probably related to the efficient horizontal mixing of potential vorticity (PV) caused by strong vortices, which tends to break the PV staircase and smooth out the jet structure. One interpretation of the Rhines jet scaling is that the jets are a natural result of PV staircases. Eastward jets correspond to the PV discontinuity whereas the westward jets correspond to the PV homogenization, and jet speed is determined by the meridional width of the PV staircase (e.g., Dritschel & McIntyre 2008;Scott & Dritschel 2012). But if the magnitude of vorticity sources is much larger than the background PV discontinuity, mixing of PV between the staircases leads to destruction of the PV staircase, therefore the jet scaling does not apply well in this situation (Scott & Dritschel 2012). In strongly forced and damped cases of Showman et al. (2019), the PV structure is disrupted by eddy vorticity and the resulting zonal jets are less sharp and meridionally broader than those that are weakly forced and damped. In our case, the vorticity mixing is even more extreme. We find that in our simulations, the vorticity associated with mid-to-high-latitude vortices are typically much stronger than either the background planetary vorticity or that associated with the zonal jets. These vortices occurs somewhat randomly in space and timetheir occurrence is not particularly constrained by the presence of zonal jets-and they can easily disrupt the zonal-mean PV structure. In panel (c) of Figure 10, we can see that the eddy relative vorticity is far greater than the background planetary vorticity in the case with drag = 10 7 s. Figure 11 shows the vertically averaged zonal-mean zonal wind in the left panel and its corresponding absolute vorticity (together with the planetary vorticity) in the right panel. Because the barotropic component dominates the jet structure in the case with drag = 10 7 s (panel (c) in Figure 9), the absolute vorticity is a good representation of PV. The eastward jets correspond to the sharp gradient of absolute vorticity at ±48 • latitude. There is a lack of other absolute vorticity staircase in other regions. Meanwhile, absolute vorticity between ±40 • latitude tends to be homogenized (although still far from being completely well mixed). In some cases, the cross-equator homogenization of absolute vorticity can lead to a strong equatorial westward jet (Dunkerton & Scott 2008). Such a tendency may be responsible for the broad, strong equatorial westward jets in our simulations.
There is a local absolute vorticity maximum northward of the jet at 48 • latitude and a local minimum southward of the jet, indicating that the jet violates the barotropic stability criteria. This is associated with the accumulation of cyclones northward of the jet and anticyclones southward of the jet. A positive feedback likely takes place in this case: the jet structure triggers accumulation of vortices, and the accumulation of those vortices further strengthens the jet structure that already promotes the accumulation. Similar influences of the migration and accumulation of vortices on the jet structure has been proposed in previous modeling studies of Jovian atmospheric dynamics (Thomson & McIntyre 2016). A balance may be reached between the accumulation of vortices and the instability associated with the jet that tends to restore the jet structure towards stability.
In box simulations with constant Coriolis parameter performed in , we showed that when the bottom frictional drag is weak ( drag = 10 7 s), kinetic energy accumulates and a pair of a cyclone and a anticyclone forms with sizes comparable to the simulated domain size. In the global domain, the presence of -effect excites Rossby waves, and eddy-mean-flow interactions channel the kinetic energy to the zonal direction, forming zonal jets instead of ever larger vortices. Our simulations demonstrate the importance of the effect on the formation of zonal jets.
Finally, equatorially trapped waves also exist in the weak-drag models but their propagation is influenced by the equatorial jets. The influence is most prominent in the case with drag = 10 7 s, in which equatorial waves are Doppler shifted by the equatorial jet with a zonal-mean zonal wind ∼ −2200 m s −1 . In the wavenumberfrequency analysis for this case, we remove the Doppler shifts associated with the jet for which we assume a uniform zonal-mean zonal wind of −2130 m s −1 according to Figure 11. The results are shown in Figure 12. Part of the signals shown in these panels are likely associated with waves propagating relative to the zonal-mean zonal wind. We find evidence of eastward Kelvin waves and westward Rossby waves in the symmetric component around the equator. A family of MRG waves is also indicated in the analysis of the antisymmetric component although their signals slightly deviate from the theoretical ones. However, in reality, the equatorial jet has horizontal shears, whose influences on waves cannot be removed by simply assuming a uniform zonal-mean zonal wind in the analysis shown in Figure 12. Such influences may be responsible for the abnormally strong westward Rossby wave signal (which is only tentative in the strong-drag cases shown in Figure 6) in the left panel, as well as for signals around zonal wavenumber 3-5 and frequency 0.1-0.2 CPD in the right panel, which is absent in the strong drag cases (Figure 7). Indeed, when we slightly change the assumed uniform zonal-mean zonal wind in the wavenumber-frequency analysis, the position or strength of the above "abnormal" signals is sensitive to the assumed zonal-mean zonal wind. But the signals associated with the Kelvin waves and MRG waves are only moderately affected, indicating that Kelvin waves and MRG waves seem robust.

Simulated lightcurves
Due to the significant inhomogeneity in the outgoing thermal flux, lightcurve variability is expected from our simulations. Figure 13 shows the simulated lightcurves normalized to their time-mean values as a function of time that is normalized by the rotation periods of the models. Panels a to d are lightcurves from models with  Figure 6 for a rotation period of 5 hours. On the right panel, dotted lines are theoretical MRG waves, same as those shown in Figure 7 for a rotation period of 5 hours. drag = 10 5 s and rotation periods of 20, 10, 5 and 2.5 hours, respectively. Panels e and f are lightcurves from models with a rotation period of 5 hours and drag = 10 6 and 10 7 s, respectively. Black, red and blue lines in each panel represent viewing angles of 0 • , 45 • and 90 • relative to the equator, respectively-0 • means an equator-on view and 90 • means a northern-pole-on view. The former maximizes the rotational modulation of the lightcurve whereas in the latter there are no rotational effects. The atmospheres in our models are statistically symmetric between the northern and southern hemisphere, and therefore simulated lightcurves viewed from two hemispheres are qualitatively similar.
The simulated lightcurves exhibit several important characteristics. For all models, the amplitudes of lightcurve variability are maximized when viewed equator-on, whereas they are minimized when viewed pole-on. Part of the reason is because the horizontal length scales of storms are larger at low latitudes, which causes larger flux perturbations when the object is viewed equator-on. This is consistent with the viewing angle dependence of near-IR flux variability amplitude found by Vos et al. (2017). The peak-to-peak amplitude of the equator-on lightcurve is typically a few percent and decreases with decreasing rotation period-the amplitude is almost 4% for the case with rotation period of 20 and 10 hours, and slightly more than 2% for the case with 5-hour rotation, and finally slightly more than 1% for the case with 2.5-hour rotation. These characteristics are further summarized in Figure 14, which shows the peak-to-peak amplitudes of normalized lightcurves with different viewing angles.
Lightcurves with rotation periods of 10, 5 and 2.5 hours exhibit clear periodicity related to the rotation period at least during certain times in their evolution. However, at some points the rotational periodicity can be complicated by the time evolution of the surface patchiness. For instance, in the case with a rotation period of 5 hours, the 5-hour periodicity is clear before the 12th rotation period, but transition to double peaks within a rotation period after that. By eye, the shape of lightcurves for a rotation period of 20 hours is more complicated by short-term irregularities. This is perhaps because the rotation period is long and the typical sizes of vortices are relatively large, which together indicate that evolution of individual vortices (typically over a timescale of several tens of hours) may significantly impact the lightcurve evolution over rotational timescales.
Viewed from pole-on, lightcurves are much smoother and show much longer periods of variation. The amplitude of poleon lightcurves decreases with decreasing rotation period. This lightcurve variability is caused by the statistical fluctuations of the ensemble of storms as discussed in . The fewer storms projected in the disk, the larger the effects that the evolution of an individual storm can have on the the lightcurve. Sizes of storms near the poles inversely decrease with decreasing rotation period, therefore the above trend found in pole-on lightcurves emerges. Even though smaller than those viewed equator-on, the peak-to-peak amplitudes of pole-on lightcurves can still reach almost 2% for cases with 20 and 10 hours and 1% for the case with 5 hours, which are detectable given sensitivities of current instruments (e.g., Wilson et al. 2014;Metchev et al. 2015). This may contribute to some long-term variations in observed lightcurves that are not easily explained by rotation modulation (e.g., some samples can be seen in Metchev et al. 2015).
Weak-drag models with drag = 10 6 and 10 7 s similarly exhibit significantly time-varying waves and vortices, and some properties of their lightcurves are similar to those with drag = 10 5 s, including peak-to-peak normalized lightcurve amplitude of a few percent, certain irregular time variability, viewing-angle dependence of variability amplitude, and complications on the periodicity due to evolution of cloud patchiness. However, because weak-drag models develop meridionally broad, strong zonal jets, zonal advection of clouds by the jets can modify the periodicity of the lightcurves. As shown in Figure 9, the equatorial jets are both westward with speeds of a few hundred m s −1 and about 2000 m s −1 in models with drag = 10 6 and 10 7 s, respectively. Although Kelvin waves have eastward phase speeds relative to the mean flow, in the rotational frame they travel to the west due to the strong westward equatorial jet. Their lightcurves, especially viewed equator-on, are expected to show periodicity longer than the rotational period of the model. This is most prominent in the case with drag = 10 7 s in panel f of Figure 13, wherein the lightcurve shows only 18 periods over 20 rotational periods.
We divide contributions to the lightcurve variability by surface inhomogeneities into two groups of dynamical features: the zonally propagating equatorial waves and mid-to-high-latitude vortices (the latter do not migrate along the zonal direction). We diagnose the effects of the two groups in the lightcurve by the following process.
First, in the model outputs, we artificially hold the equatorial region to be static. A synthetic lightcurve is generated based on this configuration, and the time evolution of the shape of the variability is caused only by evolution of mid-to-high-latitude vortices. Then, we artificially hold the mid-to-high-latitude regions to be static but do not hold the equatorial region static. The resulting time evolution of the shape of the variability represents only effects of the propagating equatorial waves. Figure 15 shows these experiments for four models with drag = 10 5 s, in which the thick grey lines are the original full lightcurves; the red lines are lightcurves wherein the equatorial regions are held fixed from time zero; and the blue lines are lightcurves wherein mid-to-high latitudes are held fixed from time zero. All cases are viewed equator-on. The equatorial region is defined as in between ±25 • latitude for the case with rotation period of 20 hours, between ±20 • latitude with 10 hours, and between ±15 • latitude with 5 and 2.5 hours. These latitudinal bands are chosen to safely include the equatorial trapped waves. Slightly changes of these latitudes do not influence our results and conclusions. In general, both the equatorial waves and mid-to-high-latitude vortices contribute to the short-term evolution of the full lightcurve variability-removing the time evolution of either component results in significant changes in the lightcurves. However, equatorial waves impact the time evolution of lightcurve shapes in more critical ways. First, the wave beating effect that causes the change of variability amplitude with time is much weaker when the equatorial regions are held fixed. This is most obvious in the case with a rotation period of 5 hours, in which the red curve (for which equatorial region is held time invariant) has almost a constant amplitude, and the splitting to double sub-peaks shown in the original lightcurve does not occur. Second, there are significant phase differences in the lightcurve variability between the full lightcurve and the lightcurve with equatorial regions fixed. This is obvious in all cases and we take the case with 2.5-hour rotation as an example: the red curve starts to show phase differences relative to the thick grey curve at a time of about the 7th rotation period; towards the end the red and grey curves show an almost 180 • phase difference. The blue curve, for which mid-to-high-latitude regions are held time invariant, mostly only show differences in the local peaks and troughs of the variability, but not in the long-term evolution of the amplitude and the phase of the variability. Our diagnostic analysis suggests that equatorial waves have major impacts on the lightcurve variability and the time evolution of the shape of the variability, and mid-to-high-latitude vortices contribute to the local features of the lightcurves. This is in good agreement with the analysis of long-term observed lightcurves of a few BDs by Apai et al. (2017), in which they found that waves can explain the major evolution of lightcurves, and "spots" are needed to fit the remaining local inconsistency between data and the wave model. Both eastward and westward waves are present at low latitudes, and the eastward Kelvin waves are sometimes more dominant as visually shown in Figure 5 and in the wavenumber-frequency analysis ( Figure  6 and 7). This implies a faster rotation of the equatorial features than the intrinsic planetary rotation, which may lead to a slightly shorter period of the lightcurve variability than the underlying rotation period. This effect has been clearly shown in the case with 2.5-hour rotation in Figure 15, and now we quantify this using periodogram analyses of the simulated lightcurves shown in Figure 13. Figure 16 shows power densities (in arbitrary unit) by the periodogram analysis for normalized lightcurves viewed equator-on for our models.
Several features are interesting in the periodograms for models with drag = 10 5 s. First, some cases exhibit slight shifts of the major peaks to a period slightly shorter than the rotation period. In particular, cases with 20 and 2.5 hours show obvious shifts. This is because of the eastwardly propagating Kelvin waves with phase speeds of a few hundred m s −1 . Second, the power densities exhibit major peaks very close to (in some cases exactly the same as) the underlying rotation period of the models. However, these peaks are broadened, indicating that the dynamical features, either eastward or westward with various phase speeds, could contaminate the periodicity of the lightcurve variability. Third, there is usually a secondary peak of the power density at a period of approximately half of the rotation period. This may be caused by the equatorial features with a zonal wavenumber 2 that appears twice to the "observer" as the objects rotate once. Finally, our models have shown that different equatorial waves that travel in different zonal velocities may contribute to the flux patchiness differently at different times ( Figure  5). This could cause time variability on the properties of lightcurves when the lightcurves are sampled over a finite period (as the real observations do). This is illustrated for the case with 5-hour rotation and drag = 10 5 s in panel c of Figure 16. We show two additional periodograms of lightcurves that sample the model outputs at different times, and each of them also samples about 20 rotational periods as the original one. Their relative shapes differ slightly, and importantly, their major power density peaks can be at periods both longer and shorter than the underlying rotation period. This suggests that not only these waves can cause differences between the "observed" lightcurve periodicity and the underlying rotation period, but also the degree of this deviation may vary with time.
For models with weaker frictional drags, the development of strong equatorial westward jets induce shifts in the period of variability.
Panel e and f in Figure 16 show periodograms for the case with drag = 10 6 and 10 7 s, respectively. The first-order feature of them is that their major peaks shift to longer periods compared to the rotation period, and the shift is stronger in the case with drag = 10 7 s.

DISCUSSION AND CONCLUSIONS
Existence of large-scale zonally propagating waves at low latitudes in our simulations opens up an avenue to understand weather on isolated BDs and directly imaged EGPs, and its consequences on the observed lightcurve variability. Rapidly evolving isotropic storms and vortices are prevalent at mid-to-high latitudes, contributing to lightcurve variability especially when the objects are viewed relatively pole-on (in which rotational modulation of lightcurve variability is diminished). Both dynamical features are driven by the cloud radiative feedback, providing an essential physical mechanism to explain several types of time evolution of lightcurve variability (see a summary and analysis in Apai et al. 2017). Eastward propagating Kelvin waves with phase speeds of a few hundred m s −1 are sometimes dominant in our simulations, and the existence of these waves in atmospheres of BDs may explain the shorter rotation period of the atmosphere than that of the interior observed for a nearby BD (Allers et al. 2020). Our models predict that different equatorial waves that may travel in different zonal velocities could influence the lightcurves differently at different times. The interesting consequence is that these waves can cause differences between the "observed" lightcurve periodicity and the underlying rotation period, and the degree of this deviation

Normalized Flux
Rotation period 2.5 hours, equator view Figure 15. Normalized lightcurves with an equator-on viewing angle for models with four different rotation periods of 20, 10, 5 and 2.5 hours (from the top to bottom panels) and with a drag timescale drag = 10 5 s. In each panel, the thick grey line is the full-model lightcurve. The red line is from thermal flux outputs wherein the equatorial region is held fixed at time zero and stays invariant, while the mid-to-high-latitude regions are not modified. On the contrary, the blue line is from outputs wherein the mid-to-high-latitude regions are held fixed at time zero. The equatorial regions here are defined as between ±25 • for the model with a rotation period of 20 hours, and between ±20 • for 10 hours, and between ±15 • for both 5 and 2.5 hours. may vary with time. It will be interesting that the same observations performed by Allers et al. (2020) could be repeated for the same system in the future to examine its long-term variability.
Recently, Vos et al. (2017), Vos et al. (2018) and Vos et al. (2020) suggested a viewing angle dependence of the observed near-IR colors and variability amplitude for a large sample of field BDs. Our dynamical models provide support for their observational results. Larger variability amplitude when viewed more equator-on is a natural result of the equatorial maximum of rotational variability along with our finding that cloud patches reach maximum sizes at low latitudes (Figure 1, 2, 13 and 14). The vertical extent of zonal-mean cloud mixing ratio is higher at lower latitudes (Figure 3), which could be responsible for the redder near-IR colors when viewed more equatoron. The near-IR colors show a wide scatter in the color-magnitude diagram for mid-to-late L dwarfs (see color-magnitude diagrams for a large sample of BDs in, e.g., Faherty et al. 2016;Liu et al. 2016). Due to the latitudinal variation of the cloud thickness at a given rotation period and the dependence of global-average cloud thickness on varying rotation period, the different viewing angles and the variation of rotation periods of the field BDs might contribute to the scatter of observed near-IR colors. Millar-Blanchaer et al. (2020) showed that assuming two broad zonal bands with different cloud properties in each hemisphere, they can reproduce time-averaged polarization measured for the nearby BD Luhman 16B. Our models do not show clear zonally banded cloud structure like those in Jupiter and Saturn, but exhibit smooth equator-to-pole cloud thickness variations. It is worthwhile to explore how can polarization be produced by the smooth equator-to-pole cloud thickness variation with different viewing angles and compare it to the measured polarization. Amplitudes of our simulated lightcurves typically range from 0.5 percent to a few percent depending on the rotation period and viewing angle, consistent with those found in the majority of observed lightcurves (Buenzli et al. 2014;Radigan et al. 2014;Radigan 2014;Wilson et al. 2014;Metchev et al. 2015;Vos et al. 2018). We did not reproduce variability with amplitude much greater than a few percent. Yet, several field BDs and free-floating low-mass objects have shown large variability with peak-to-peak amplitudes greater than 10 percent (e.g., Buenzli et al. 2012;Apai et al. 2013;Buenzli et al. 2015;Lew et al. 2016;Apai et al. 2017;Zhou et al. 2020). Either cloud patches with sizes larger than those in our simulations or a greater horizontal contrast of outgoing thermal flux are required to explain those very large varability amplitudes. Our models only occupy a very limited parameter space. Further exploring parameter space, including broad range of surface gravity, atmospheric temperature and varying cloud particle size, will yield richer dynamical behaviors.
Our models are highly idealized in the sense that radiative transfer and cloud formation (including the cloud microphysics and the treatment of sub-grid-scale structure) are vastly simplified in order to emphasize the role of atmospheric dynamics. Our modeled cloud structures capture the first-order behavior of cloud formation that has been shown by previous cloud formation models of BDs-a sharp cloud base typically emerges around the condensation level and the cloud mixing ratio smoothly decreases with decreasing pressure due to mixing (see reviews by, e.g., Helling & Fomins 2013;Marley & Robinson 2015;Helling 2019). Our cloud scheme neglects temperature feedback on cloud formation, which could influence locations of the cloud base, and it certainly does not capture all the sophisticated microphysical processes. Future endeavours should include better representation of radiative transfer, cloud microphysics and parameterization of sub-grid cloud structure.
Finally, we summarize our major findings in this study as follows: • Vigorous atmospheric circulation can be triggered and maintained by cloud radiative feedback in conditions appropriate for BDs and directly imaged EGPs. When the bottom frictional drag is strong, zonal flows in the deep layers of our models are weak (with speeds much smaller than 100 m s −1 ). In the observable layer where clouds form, mid-to-high latitudes are dominated by isotropic vortices, with thick clouds forming in anticyclones and thin clouds or clear sky in cyclones. This is consistent with the results of previous -plane models . At low latitudes, large-scale zonally propagating waves with both eastward and westward phase speeds are the dominant dynamical feature. At a given rotation period, sizes of storms and vortices are typically larger at lower latitudes than those at higher latitudes. The overall sizes of storms are larger when the rotation period is longer.
• Lightcurves from our simulations have amplitudes from 0.5 percent to several percent, consistent with the majority of observed lightcurves. For a given rotation period, the lightcurve amplitude decreases with increasing viewing angle (0 • means equator-on and 90 • means pole-on), while it typically increases with increasing rotation period at a given viewing angle. When the bottom drag is strong, zonally propagating waves at low latitudes have typical phase speeds of a few hundred m s −1 . They can cause short-term evolution of lightcurves via wave beating effects, qualitatively similar to the observed lightcurves of some field BDs as summarized in Apai et al. (2017). The eastward Kelvin waves can cause the equatorial flux inhomogeneity rotating faster than the underlying planetary rotation, which may help to explain the observed shorter rotation period of atmospheric features than that of the interior (Allers et al. 2020). Isotropic storms and vortices at mid-to-high latitudes also contribute to the lightcurve variability.
• Clouds are generally mixed to higher altitudes near the equator than at high latitudes due to the stronger effect of rotation at high latitudes. This supports the observed redder IR colors for objects viewed more equator on (Vos et al. 2017;Vos et al. 2020).
• We expect that the strong-drag models might be appropriate for BDs and directly imaged EGPs because efficient convective mixing in the interior is expected to suppress strong zonal flows near the top of convective zone. But we still perform experiments with weaker drags to explore dynamics in the weak-drag regime. We find that robust zonal jets with speeds from several hundred to more than 2000 m s −1 can form in our weak-drag models, with typically a broad westward equatorial jet and a high-latitude eastward jet in each hemisphere. Similar to the weak-drag cases, vortices form at mid-to-high latitudes and equatorially trapped waves form at low latitudes. Both the zonal propagation of equatorial waves and the spacial distribution of vortices are affected by the presence of strong jets. Simulated lightcurves show longer periodicity than the underlying rotation period due to the strong westward equatorial jets. The meridionally broad jet structure may be related to the efficient potential vorticity mixing associated with intense eddies.
• The origin of the equatorially propagating waves in our simulations is likely related to the self-excitation by cloud radiative feedbacks. Physical properties of the equatorially symmetric eastward waves resembles properties of both adiabatic free Kelvin waves and forced-damped waves triggered by a stationary equatorial heat source. Linear stability theory of waves coupled with cloud radiative feedback may help to explain the origin of our simulated waves but fails to explain their propagation.

APPENDIX A: HOVM OLLER DIAGRAMS AT MID LATITUDES
In Figure A1, we show the Hovm oller diagrams of outgoing thermal flux at 45 • latitude as a function of longitude and time for models with different rotation period. These regions are dominated by stochastically evolving vortices, and there is no systematic eastward or westward propagation seen in these diagrams, which is in stark contrast to the equatorial disturbances shown in Figure 5.

APPENDIX B: CLOUD RADIATIVELY INDUCED EQUATORIAL WAVES
Equatorially trapped waves in our simulations play vital roles in the evolution of simulated lightcurves, and may be responsible for the puzzling time evolution of some observed lightcurves (Apai et al. 2017). Given their central importance, we explore their nature using a linear wave theory that is coupled to cloud radiative effect. We will show that the linear theory may be used to explain the energetic sources of our simulated waves, but cannot capture the propagating nature of the simulated waves.
The analysis of cloud radiatively induced instability was first carried out by Gierasch et al. (1973) in the quasi-geostrophic system and inertia gravity wave ( -plane) system. 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 mode. In a −plane approximation, inertia gravity waves with small horizontal wavelengths can host purely unstable modes but unstable modes cease at larger horizontal scales. In a quasi-geostrophic flow, unstable modes are possible for axially symmetric flows and flows with nonzero (Gierasch et al. 1973). We now extend this theory to the equatorial plane.

B1 Cloud radiative effect upon thermodynamics
We first present the thermodynamics related to cloud radiative effect following Gierasch et al. (1973). Assuming that the change of outgoing thermal flux is due to brightness temperature deviation that is caused by either actual temperature variation or the cloud-top altitude variation, we can write the change of emitted flux as where is the brightness temperature, is the Stefan-Boltzmann constant, is the change of cloud-top altitude, and Γ = | |. Because of the change of outgoing radiative flux, the atmospheric column is no longer in equilibrium and so net heating/cooling must occur. Denoting as the heating rate per unit mass, the gas density and the cloud thickness, we have For analytic simplicity, is assumed to obey the following linear relation at any level: where = ∫ 0 . Moreover, clouds are assumed to be perfectly advected by the flow, i.e., sedimentation is negligible compared to the vertical advection, so that the rate of change in cloud-top altitude Then we assume a basic state of the atmosphere at rest with small horizontal temperature perturbations. The linearized equation of thermodynamics in height coordinates expressed using temperature is written: Combining with Eq. (B4) we have the linearized thermodynamic equation where Γ = + , and 1 = 4 3 . Here represents a radiative timescale of the cloud-forming atmosphere.
By just the terms involving vertical velocity in Equation (B6), Γ = 1 Γ , an instability is obvious with a growth rate of = Γ Γ . Now we consider the full linearized dynamical equations below.

B2 Rossby, mixed Rossby-gravity and inertia-gravity modes
We start with the linearized dynamical equations with a basic state at rest in an equatorial plane, where the Coriolis parameter is written = . The zonal and meridional angular momentum , hydrostatic balance, and continuity equations are, respectively where 0 is a reference temperature, = − ln( / 0 ) is the vertical coordinate at log pressure, = 0 / is the scale height, and is the geopotential. Equations (B6) to (B10) forms a closed set.
We seek wave-like solutions: where˜( ),˜( ),˜( ),˜( ) and˜( ) are functions of only. To satisfy the form of equatorially trapped waves, the boundary condition is typically applied to ensure that the disturbances vanish when | | → ∞. Repeating the derivations for the classic equatorial waves (e.g., Matsuno 1966;Holton & Hakim 2012), we arrive at an equation for˜( ): where 2 = 2 / 2 , 2 = Γ/ 0 , = /L , and L is the equatorial deformation radius modified by the diabatic cloud radiative effect: Note that in the equatorial theory of Hayashi (1970) that considered latent heating effects, a complex equatorial deformation was also possible. In the adiabatic limit of → ∞, the deformation radius recovers the classic definition L = √︁ / and the dispersion recovers the classic dispersion relation of adiabatic, unforced equatorially trapped waves (Matsuno 1966). In the limit of → 0, there is no solution that satisfies the boundary condition. The following relation has to be met for the given boundary condition: Then, the solution has the form ( ) =˜0H ( ) exp(− 2 /2), where˜0 is a constant with units of velocity, and H ( ) designates the th Hermite polynomial. In addition, solutions satisfy the boundary condition if exp(− 2 /2) diminishes when | | → ∞, which requires that the real part of L 2 is positive.
We first seek solutions with real, positive . Possible solutions should be in the range between 0 and . Outside this range, there is no solution that satisfies a real , which is required to have wavelike zonal disturbances. Equating the imaginary part in the LHS of Equation (B14) to zero, one obtains 3 + 1 2 + 2 2 − 2 2 = 0.
Only when = 0 (the mixed Rossby gravity modes), both equations (B16) and (B17) can be simultaneously satisfied, and there is a set of purely unstable, non-propagating = 0 modes. In this case, L 2 is purely imaginary and fails the strict meridional boundary condition. Nevertheless, by solving equation (B14) numerically, we cannot find other modes that contain a positive real part of . All other modes have both negative real parts and imaginary parts, and they resemble the classic equatorial waves-the propagating Rossby, MRG and inertia-gravity modes but with damping in the wave amplitudes due to thermal radiation. The dispersion relations for all possible solutions are plotted as solid curves in Figure B1, in which the top panel shows the − (in the same format as that plotted in classic literature of equatorial waves, e.g., Holton & Hakim 2012) and the bottom panel shows the , where = + . Despite that the special non-propagating unstable mode in Equation (B16) does not strictly satisfy the lateral boundary condition, at least their disturbances do not amplify with | | → ∞. So we still plot the dispersion relation as the dotted lines in Figure B1.

B3 Kelvin modes
The Kelvin mode is a special mode in which the meridional velocity is zero. The governing equations for the zonal and meridional angular momentum, and continuity are simplified to We make use of equations (B18), (B9), (B20) and (B6) and assume wave-like solutions of equation (B11), then we obtain the following dispersion relation 3 + 1 2 + 2 2 − 2 2 = 0.
This is the same as Equation (B16) for = 0 modes, and growing modes with positive exist. Additional constraints from the lateral where˜0 is the amplitude of the perturbation zonal velocity at the equator, and is written as = + . If one restricts˜to vanish with | | → ∞, < 0 needs to be satisfied. With this regard, the purely unstable modes fail to satisfy the boundary condition, and only the eastward propagating but decaying modes are valid solutions. This reaches the same conclusion as the = 0 unstable modes. The dispersion relation of the unstable Kelvin modes are also represented as dotted lines in Figure B1.
Our key finding is that there is only one set of growing but nonpropagating modes corresponding to the Kelvin and = 0 modes that may be marginally relevant. This may provide a way to excite the Kelvin waves, = 0 MRG waves and = 0 eastward inertia gravity waves seen in our simulations. Other propagating modes have properties quantitatively similar to the adiabatic free modes but with damping of their amplitudes due to thermal radiation. This paper has been typeset from a T E X/L A T E X file prepared by the author.