Jet Streams and Tracer Mixing in the Atmospheres of Brown Dwarfs and Isolated Young Giant Planets

Observations of brown dwarfs and relatively isolated young extrasolar giant planets have provided unprecedented details to probe atmospheric dynamics in a new regime. Questions about mechanisms governing global circulation remain to be addressed. Previous studies have shown that small-scale, randomly varying thermal perturbations resulting from interactions between convection and the overlying stratified layers can drive zonal jet streams, waves, and turbulence. Here, we improve upon our previous general circulation model by using a two-stream grey radiative transfer scheme to represent more realistic heating and cooling rates. We examine the formation of zonal jets and their time evolution, and vertical mixing of passive tracers including clouds and chemical species. Under relatively weak radiative and frictional dissipation, robust zonal jets with speeds up to a few hundred $\rm m\;s^{-1}$ are typical outcomes. The off-equatorial jets tend to be pressure-independent while the equatorial jets exhibit significant vertical wind shear. Models with strong dissipation inhibit jet formation and have isotropic turbulence in off-equatorial regions. Quasi-periodic oscillations of the equatorial flow with periods ranging from tens of days to months are prevalent at relatively low atmospheric temperatures. Sub-micron cloud particles can be transported to several scale heights above the condensation level, while larger particles form thinner layers. Cloud decks are inhomogeneous near their cloud tops. Chemical tracers with chemical timescales $>10^5$ s can be driven out of equilibrium. The equivalent vertical diffusion coefficients, $K_{\mathrm{zz}}$, for the global-mean tracer, are diagnosed from our models and are typically on the order of $1\sim10^2\rm m^2\;s^{-1}$. Finally, we derive an analytic estimation of $K_{\mathrm{zz}}$ for different types of tracers under relevant conditions.


INTRODUCTION
Our understanding of the atmospheres of brown dwarfs (BDs) has grown sophisticated owing to the increasing range of observational constraints. Several types of observations have suggested that their atmospheres are dynamic, inhomogeneous and in disequilibrium (see a summary in Tan & Showman 2021b). In particular, several recent surveys that aim at monitoring the near-infrared (IR) brightness of BDs offer perhaps the strongest support for the presence of atmospheric circulation (e.g., Vos et al. 2022Vos et al. , 2020Manjavacas et al. 2019a;Vos et al. 2018;Apai et al. 2017;Yang et al. 2016;Rajan et al. 2015;Metchev et al. 2015;Wilson et al. 2014;Radigan et al. 2014;Buenzli et al. 2014;Heinze et al. 2013;Khandrika et al. 2013). They revealed that low-amplitude (∼ 1%) flux variability is likely common among L and T dwarfs, while a small fraction of BDs exhibit large-amplitude flux variations (e.g., Artigau et al. 2009;Radigan et al. 2012;Yang et al. 2016;Apai et al. 2017;Lew et al. 2020;Zhou et al. 2020). These motivate a greater need to better understand their global atmospheric circulation and the consequences for variability, cloud formation, and ★ E-mail: xianyu.tan@physics.ox.ac.uk mixing of chemical species (Showman et al. 2020;Zhang 2020). Extrasolar giant planets (EGPs) discovered by direct imaging so far are mostly young, far away from their host stars or isolated. They share similarities with BDs in terms of their spectral type, near-IR colors, inference of clouds and chemical disequilibrium (e.g., Chauvin et al. 2017;De Rosa et al. 2016;Skemer et al. 2016;Macintosh et al. 2015;Bonnefoy et al. 2014;Marley et al. 2012;Barman et al. 2011;Currie et al. 2011). Near-IR brightness variability has also been observed on these young isolated EGPs (Biller et al. 2015;Zhou et al. 2016;Biller et al. 2018;Miles-Páez et al. 2019;Manjavacas et al. 2019b;Lew et al. 2020;Zhou et al. 2020).
Advancements of observations for BDs and relatively isolated young EGPs provide ever unprecedented details to probe atmospheric dynamics in a new dynamical regime. This regime is typically characterized by fast rotation, vigorous internal thermal forcing and negligible external forcing. However, these objects receive negligible external stellar irradiation, and therefore lack the large-scale stellar heating contracts that are responsible for driving the global circulation like those on close-in exoplanets and solar system planets (Showman et al. 2020). The first-level question needed to be addressed is, 2 X. Tan what are the possible mechanisms driving the global circulation and what are the resulting properties?
Convection in the interior of BDs is vigorous and this convection is expected to perturb the overlying, stably stratified atmosphere, generating large-scale atmospheric circulation that consist of turbulence, waves, and zonal (east-west) jet streams ). There have been several studies along this direction. Hydrodynamic simulations in a local enclosed area by Freytag et al. (2010) and Allard et al. (2012) showed that gravity waves generated by interactions between the convective interior and the stratified layer can cause mixing that leads to small-scale cloud patchiness.  presented global convection models, and they analytically estimated the typical wind speeds and horizontal temperature differences driven by the absorption and breaking of atmospheric waves in the stably stratified atmosphere. By injecting random forcing to a shallow-water system that represents the effect of convection perturbing the stratified atmosphere, 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 mean meridional circulation are common outcomes of the dynamics. They also revealed long-term (multiple months to years) quasi-periodic oscillations on the equatorial zonal jets, and the mechanism driving the oscillations is similar to that of the Quasi-biennial oscillation (QBO) in Earth's stratosphere.
An alternative class of mechanisms invoking radiative feedback by clouds has been shown to be essential to induce instantaneous atmospheric variability and vigorous global circulation when the atmospheres are dusty Tan & Showman 2021a,b). These models generate large-scale cloud and temperature patchiness, vortices, waves and jet streams, as well as systematic equator-to-pole variation of cloud thickness. Latent heat released from condensation of silicate and iron vapor also helps to organize zonal jets and patchy storms . Similar radiative processes associated with clouds or chemistry may also induce small-scale instability (Tremblin et al. 2019(Tremblin et al. , 2021. There are two main goals in this study. First, we further investigate the general circulation of BDs and isolated young EGPs driven by thermal perturbations using an updated general circulation model. This study is along the line of Zhang & Showman (2014) and Showman et al. (2019), but we utilize a more realistic model to examine the findings revealed by previous studies. Our models in Showman et al. (2019) implemented an idealized relaxation scheme with a uniform radiative timescale throughout the atmosphere. While that scheme provides a good control of models over the parameter space which is essential to demonstrate dynamical mechanisms, it oversimplifies the radiative heating and cooling rates, and the relevance of these previous results to real objects needs to be validated against a more realistic model setup. Here, we adopt an idealized radiative transfer scheme that calculates the three-dimensional heating and cooling rates in a more realistic and self-consistent way which is especially important in the current context that we investigate models with different effective temperatures. This enables us to examine the characteristic behaviors of jet formation and their long-term evolution under more relevant conditions of BDs and isolated young EGPs. The second goal is to characterize the vertical mixing of passive tracers, including cloud particles and chemical species, by large-scale circulation in the stratified layers where convection is not directly responsible for mixing. This has not been investigated in the context of BDs and isolated young EGPs. Our study will demonstrate a dynamical mechanism for tracer mixing in the stratified layers, which is relevant in understanding the dusty L dwarfs and the prevalent non-equilibrium chemistry in atmospheres of BDs and isolated EGPs. This paper is organized as follows. We introduce our general circulation model in Section 2; we then present results of zonal jet formation in Section 3 and explore vertical mixing of passive tracers in Section 4; finally we discuss and summarize our results in Section 5.

General Circulation Model
We model a global three-dimensional (3D) thin atmosphere using a general circulation model (GCM), the MITgcm (Adcroft et al. 2004, see also mitgcm.org), which solves the hydrostatic primitive equations of dynamical meteorology in pressure coordinates. Our model is similar to that of Showman et al. (2019) but we couple a grey radiative transfer (RT) scheme to the dynamical core instead of the idealized Newtonian cooling scheme. In brief, thermal radiative fluxes within the atmosphere are obtained by solving the plane-parallel, two-stream approximation of the RT equations using an efficient and reliable numerical tool twostr (Kylling et al. 1995). Then the radiative heating rates are calculated by taking the vertical divergence of fluxes to drive dynamics. The detailed implementation and demonstrations of our RT scheme is described in Komacek et al. (2017) and Tan & Komacek (2019) on the application of hot Jupiters and in Tan & Showman (2021a,b) on BDs. The lower boundary condition is a prescribed, globally uniform temperature at the bottom pressure which is intended to mimic a uniform-entropy deep layer resulting from efficient convective mixing. We include gas opacity only and omit scattering. The Rosseland-mean gas opacity as a function of temperature and pressure from Freedman et al. (2014) is coupled to the RT calculation. Our radiative-convective equilibrium RT calculations generate temperature-pressure (T-P) profiles with reasonable radiative-convective boundaries (RCBs) compared to the more sophisticated one-dimensional (1D) atmospheric models (e.g., Allard et al. 2001;Tsuji 2002;Burrows et al. 2006;Marley et al. 2021).
The deep layers in our simulation domain reach the convectively unstable regions. Small-scale convective transport of entropy and tracers that cannot be resolved in the model is parameterized using a simple column convective adjustment scheme that instantaneously adjusts the unstable region to be convectively neutral and homogenize tracers within the adjusted region. If any adjacent two layers within a vertical atmospheric column are unstable, they are instantaneously adjusted to a convectively neutral state while conserving total enthalpy. The whole column is repeatedly scanned until convective instability is eliminated everywhere. Tracers are also well homogenized within the adjusted domain. Similar simple schemes have been adopted in some GCMs for Earth (e.g., Collins et al. 2004), exoplanets (e.g., Deitrick et al. 2020;Lee et al. 2021) and our previous brown dwarf models (Tan & Showman 2021a,b).
A Rayleigh drag is 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 magnetohydrodynamics (MHD) drag. The drag linearly decreases with decreasing pressure and takes the form as F drag = − ( )v, where v is the horizontal velocity vector and ( ) is a pressure-dependent drag coefficient. ( ) decreases from 1/ drag at the bottom boundary bot to zero at certain pressure drag,top , where drag is a characteristic drag timescale. No drag is included at pressures lower than drag,top . In all simulations, we fix drag,top to 50 bars. The bottom pressure is 100 bars. Kinetic energy dissipated by the Rayleigh drag is converted to thermal energy. As in Showman et al. (2019), the drag timescale drag in this study is treated as a free parameter to explore possible circulation patterns. As shown below, a weak drag is essential for jet formation.
We parameterize the global thermal perturbations resulting from interactions between convection and the overlying stratified layers using the scheme in Showman et al. (2019). The essence of this scheme is the stochastic and isotropic nature of the forcing, such that the resulting jets and their time evolution are purely results of the internal dynamics rather than being related to the nature of the forcing. Similar types of schemes have long been used to understand turbulence and jet formation on solar-system giant planets (see reviews of Vasavada & Showman 2005;Galperin et al. 2008;Showman et al. 2018 and many references wherein). The work in Showman et al. (2019) and here is an extension into the brown dwarf regime.
The procedure is summarized as follows. The temperature field is consecutively forced by a global 3D forcing pattern ( , , , ) as a function of longitude , latitude , pressure and time . At each time step, a forcing function is constructed by superposing a sequence of fully-normalized spherical harmonics with total wavenumbers = − to = and the zonal wavenumbers from 1 to . The phases of each zonal mode are randomized, and their amplitudes are statistically equal, such that their superposition ( , ) is spherically isotropic with a characteristic wavelength determined by and . Then the forcing function is normalized to a unit of K s −1 by multiplying a forcing amplitude . The horizontal forcing pattern ℎ ( , , ) is marched forward with time using a Markov process: where = / is a de-correlation factor, is a dynamical time step and is a storm timescale which characterizes the evolution timescale of the perturbation pattern. A similar horizontal forcing was also used in shallow-water models of Zhang & Showman (2014). In the vertical direction, the thermal perturbations are applied in regions within two scale heights away from the RCB of the equilibrium profile. The forcing amplitude decreases as the pressure deviates from the RCB assuming a normalized vertical forcing profile ( ): where rcb is the pressure at the RCB. Finally, the 3D forcing pattern ( , , , ) is simply ( ) ℎ ( , , ). The stochastic perturbation scheme is needed to drive the circulation given the rest of the model setup. The convective adjustment scheme is essentially a damping scheme that avoids dry superadiabatic profiles in our hydrostatic GCM. With the isotropic bottom temperature and the lack of external irradiation on brown dwarfs, the temperature-pressure profiles are relaxed towards globally isotropic, thereby damping out isobaric temperature variations and dynamics. Indeed, experiments in which the stochastic forcing is turned off show essentially no dynamics but only radiative-convective equilibrium profiles.
In this study we mainly explore atmospheric circulation of three effective temperatures, eff = 1500, 1000 and 500 K, whose equilibrium profiles are shown in the left panel of Figure 1. These represent objects with temperature of typical L, T and early Y type spectrum (Kirkpatrick 2005;Cushing et al. 2011). The vertical forcing regions are indicated by the thick lines in Figure 1. The right panel shows the corresponding profiles of the Brunt-Vaisala frequency defined as = √︂ 2 − ln ln where is gravity, is the specific gas constant and is the specific heat at constant pressure. is nearly zero in the deep convectively neutral regions, then becomes positive above the RCB with a transition to √ in the nearly isothermal layers. Note that for the case with eff = 500 K, a detached stratified zone appears at around 50 bars. This is not uncommon in models for slightly cooler objects (e.g., Marley et al. 2021). We consider storm timescale = 10 4 , 10 5 and 10 6 s, and drag timescale drag = 10 5 , 10 6 and 10 7 s. We assume a surface gravity of 1000 m s −2 , a radius of 6 × 10 7 meters and a rotation period of 5 hours in simulations presented in Section 3. We adopt = 30 and = 1 for most simulations in order to have sizes of perturbations much smaller than the radius of the objects in an affordable model resolution.
We solve the equations using the cubed-sphere grids (Adcroft et al. 2004). A fourth-order Shapiro filter is applied in the momentum and thermodynamics equations to maintain numerical stability. The nominal horizontal resolution of our simulations is C96, which is equivalent to 384 × 192 in longitude and latitude. We have tested models with higher horizontal resolution, and the large-scale circulation features are insensitive to the resolution. The simulated pressure domain is between 10 −3 to 100 bars, and is evenly discretized into 80 layers in log-pressure space. The ideal gas law of a hydrogen-helium mixture is applied to the equation of state, with the specific heat = 1.3 × 10 4 Jkg −1 K −1 and specific gas constant = 3714 Jkg −1 K −1 . With somewhat strong forcing in many of our models, those with a strong drag (a drag timescale of 10 5 s) equilibrate rather quickly after about 300 simulation days. Models with longer drag timescales spin up slower; for example, the ones with a drag timescale of 10 7 s require more than 1500 simulation days to equilibrate. After the models reached statistical equilibrium, we ran an additional 500 to 1000 simulations days for time averaging. Because our models are somewhat strongly forced-dissipated, their spinup time is considerably shorter than models in the regime of solar-system giant planets ).

Passive Tracers
To study the effects of stratospheric large-scale circulation on the mixing of cloud particles and chemical species, we implement two types of tracers cond and chem in our models. The tracers are assumed to be passively advected by the flow and do not have feedbacks on the atmospheric structure and dynamics. In reality, cloud radiative effect and latent heating would excert feedbacks to shaping the temperature structure and dynamics Tang et al. 2021;Tan & Showman 2021a,b;Tremblin et al. 2021). The sinks and sources of tracers in our models are idealized (as will be described below) despite that processes governing clouds and chemistry in BDs and EGPs are highly complex (e.g., Helling et al. 2008;Visscher et al. 2010;Helling & Casewell 2014;Marley & Robinson 2015;Gao et al. 2021). In this work, we step aside these complexities just to focus on the passive tracer transport problem. This is a necessary step because only after we understand the passive tracer problem we would appreciate tracer transport in models with more self-consistent coupling. These simplifications follow the merits of conceptual modeling for 3D tracer mixing in exoplanet atmospheres (Parmentier et al. 2013  The first type of tracer cond represents the mixing ratio of a condensable species relative to its abundance in the deep atmosphere. A critical condensation pressure cond is prescribed, and at pressures higher than cond , cond is in the form of vapor and subjected to replenishment by efficient convective mixing from the interior. At pressures lower than cond , cond is in the form of solid particles and is subjected to gravitational settling. This is equivalent to assuming that the scale height of the saturation vapor pressure is much smaller than a pressure scale height. The governing equation of cond is: where / = / + v · ∇ + / is the material derivative, ∇ is the isobaric gradient, is the vertical velocity in pressure coordinates, deep is the relative mixing ratio in the deep layers and is assumed to be 1, rep is the replenishment timescale which is set to be 10 3 s, and is the particle settling speed in pressure coordinates as a function of particle radius, pressure and temperature (with positive pointing downward to larger pressure). We adopt the formula summarized in Ackerman & Marley (2001), Spiegel et al. (2009) and Parmentier et al. (2013) to calculate the particle settling velocity in height coordinates, and then is approximated using the hydrostatic balance. The density of cloud particles is assumed to be 3190 kgm −3 as that of MgSiO 3 , which is one of the most abundant condensates in L and T dwarfs (Ackerman & Marley 2001). In models with eff = 1000 K, the condensation pressure cond is assumed at 10 bars which is near the RCB. This setup ensures that any vertical transport at pressures lower than cond is due to resolved large-scale dynamics rather than parameterized convective transport. The particle radius is treated as a free parameter and varies from 0.05 to 5 m.
The second type of tracer chem represents the mixing ratio of hypothetical chemical species relative to its deep abundance. A chemical equilibrium profile eq is assumed for the tracer. When the chemical species is driven out of equilibrium due to dynamics, it is relaxed back to eq over a characteristic chemical timescale chem . The governing equation for chem is simply: eq is assumed to vary with pressure alone which is reasonable in conditions relevant to BDs and isolated EGPs. We define a piecewisecontinuous function for such a transitional behavior: eq is assumed to be eq,top at pressures lower than top and eq,bot at pressures higher than bot ; then log eq varies linearly with log in between top and bot . In this work, we fix the following parameters as eq,bot = 1, eq,top = 10 −6 , bot = 10 bars and top = 10 −3 bar. Following Zhang & Showman (2018a), in a main suite of models, the chemical timescale chem is assumed independent of pressure and vary from 10 4 to 10 7 s. For many species relevant to BDs and isolated EGPs, the chemical timescale chem is expected to be short in the deep hot atmosphere and long in the cooler upper atmosphere. Therefore, in another suite of models, chem is assumed to be a function like that of eq with chem,bot = 10 3 s and chem,top varying from 10 4 to 10 7 s.
Small-scale gravity waves that cannot be resolved in global models can be potentially important in transport of tracers and angular momentum via small-scale turbulence and interactions with mean flows. These effects can only be parameterized, a procedure common in many GCMs, albeit fraught with uncertainties associated with poorly constrained free parameters related to generation and properties of the waves, assumptions about their dissipation and interactions between the mean flow. In this study we omit gravity-wave parameterizations for two reasons. The first is to ensure a clean environment to understand the mixing by resolved large-scale advection. Secondly, the sub-grid gravity wave parameterization is yet highly uncertain in conditions appropriate for BDs and EGPs. The mixing of tracers from our current models may thus be considered as a baseline in real conditions.

Zonal jet streams
Kinetic energy injected by thermal perturbations organize itself towards large-scale, while imposed isotropic radiation and bottom frictional drag dissipate the kinetic energy. The organization of the largescale flow is determined by the relative strengths of the forcing and damping processes. In this section, we examine the formation of zonal jets spanning a wide range of parameter space and have the following key findings. First, although the forcing and damping are horizontally isotropic, robust zonal jet streams can emerge from the interaction of the waves and turbulence with the planetary rotation under relatively weak radiative and frictional dissipation. At a particular atmospheric temperature with a fixed thermal perturbation rate expected for real objects, a weak bottom frictional drag favours jet formation at all latitudes, while models with strong bottom drag exhibit isotropic turbulence at mid-to-high latitudes and jets only near the equator at altitudes above the drag domain. Then, with a fixed thermal perturbation rate and a relatively weak frictional drag, low atmospheric temperature (weak radiative damping) promotes an overall stronger circulation and zonal jets, while high temperature (strong radiative damping) leads to a much weaker circulation. With a weak frictional drag and a low atmospheric temperature, the jet speed increases and the number of jets decreases with increasing thermal perturbation rates. Finally, in models that form robust zonal jets, the off-equatorial jet structures have significant pressure independent components but the equatorial jets usually exhibit significant vertical jet shears. Because the atmosphere is rotation dominated at mid-to-high latitudes, the mean meridional temperature variations associated with the off-equatorial jets are small compared to smallscale variations. Most of our findings are qualitatively consistent with those found in Showman et al. (2019). Below we describe and discuss these results in detail.
We first present a set of models with an effective temperature eff = 1000 K and a thermal forcing amplitude = 5 × 10 −4 Ks −1 , but with three bottom drag timescales of drag = 10 5 , 10 6 and 10 7 s. At statistical equilibrium, this results in typical isobaric temperature variation of 80 ∼ 100 K near the RCB, consistent with the estimated values of 50 to 200 K for brown dwarf conditions (see appendix in Showman et al. 2019). Figure 2 shows the instantaneous temperature at 8 bars (a level slightly above the RCB) on the left column and zonal velocity at 8 bars on the right column. Quantities at lower pressures are qualitatively similar to those shown in Figure 2. Figure 3 shows the time-averaged zonal-mean zonal winds as a function of latitude and pressure. These results are obtained after the models reaching statistical equilibrium.
The model with a strong bottom drag of drag = 10 5 s does not form jets at pressures higher than a few bars as well as throughout mid-to-high latitudes (top panel of Figure 3). At pressures lower than about 1 bar, alternating eastward and westward zonal flows with amplitude ∼ 150 m s −1 emerge at low latitudes. At pressures around 8 bars, the circulation is mainly consisted of quasi-isotropic turbulence. Temperature anomalies are nearly isotropic and their spacial distributions are not very different to the forcing pattern at mid-to-high latitudes. At low latitudes within about ±10 • , amplitude of the temperature anomalies appears to be much smaller than that at mid-to-high latitudes, and the shape is distorted from the forcing pattern. The zonal velocity at 8 bars exhibits a larger magnitude at low latitudes than that at high latitudes.
Dynamics at mid-to-high latitudes is rotationally dominated on rapidly rotating BDs and giant planets .
The importance of rotation is characterized by the Rossby number = / , where is a characteristic horizontal wind speed, = 2Ω sin is the Coriolis parameter, Ω is the rotation rate and is a characteristic length scale of the flow. In the regime with 1, the Coriolis force mainly balances the horizontal pressure gradient force, and large horizontal temperature difference can be supported (Charney 1963). Near the equator where advection is significant in the horizontal momentum balance, 1 and we expect a "tropical" regime in which gravity waves efficiently remove temperature anomalies, leading to a weak horizontal temperature gradient (Sobel et al. 2001). Adopting the RMS of horizontal wind speed ∼ 300 ms −1 , ∼ 10 7 m from our models, the Rossby number is ∼ 0.5, 0.25 and 0.12 at ±5 • , ±10 • and ±20 • latitudes, respectively. This dynamical regime transition corresponds to the qualitative difference of temperature field between low and mid-to-high latitudes.
With weaker bottom drags of drag = 10 6 and 10 7 s, zonal jets emerge, and the jet speed increases with increasing drag (see Figures  2 and 3). In the case with drag = 10 6 s, eddy wind speeds are still comparable to the jet speeds of about 250 m s −1 , and the jets show significant distortions as seen in panel (d) of Figure 2. Kinetic energy associated the eddies is comparable to that associated with the jets. In the model with drag = 10 7 s, zonal jets are more robust and stable with a peak speed reaching up to 500 m s −1 , exceeding the RMS eddy speed. Kinetic energy is mostly contained in the jets. Vigorous eddy activities are still clearly seen along with the jets, and these are likely waves triggered by the thermal perturbations and breaking near the flank of jets. No stable long-lived vortex is maintained. The equatorial jets exhibit significant vertical wind shears, similar to that in the strong-drag model. However, the off-equatorial jets develop significant pressure-independent components as shown in Figure 3.
The isobaric temperature anomalies in the weak-drag models are qualitatively similar to those in strong-drag models. There is still subtle difference in the case with drag = 10 7 s: the morphology of temperature anomalies is sheared and their amplitude is weakened by the strong zonal jets compared to those with a strong drag. The temperature fields at off-equatorial regions do not exhibit the same zonal banded structure as in the zonal velocity field. This is a property of the thermal wind balance which holds well in the regime of 1 (e.g., Holton 1973, Chapter 3): where and are zonal-mean zonal wind and temperature, and is northward distance. The largely pressure-independent off-equatorial jets are in line with the small meridional temperature gradient at off-equatorial regions in our models. Next, we examine the effect of different atmospheric temperature (therefore different radiative damping rate) on the circulation pattern. We perform a set of models with a thermal forcing amplitude = 1.5 × 10 −4 Ks −1 , a moderate drag timescale of drag = 10 6 s, and with three effective temperatures eff = 500, 1000 and 1500 K. The radiative damping timescales near their photospheres can be approximated as rad ∼ 4 3 (e.g., Showman & Guillot 2002) where is the pressure near the photosphere, is gravity and is the Stefan-Boltzmann constant. We adopt pressures near their photospheres and ∼ eff , yielding rad ∼ 10 5 , ∼ 1.5 × 10 4 and ∼ 6 × 10 3 s for models with eff = 500, 1000 and 1500 K, respectively. The equilibrium isobaric temperature variations reach about 60 K near the RCB in the model with eff = 500 K, but drastically decrease with increasing eff due to the rapidly increasing radiative damping. Figure 4 shows the instantaneous isobaric zonal velocity of these (a) drag = 10 5 s.

Figure 2.
Snapshots of temperature (left) and zonal winds (right) for models with three different bottom drag timescales of drag = 10 5 s (top), 10 6 s (middle) and 10 7 s (bottom). The quantities are shown at 8.2 bars. These models have an effective temperature eff = 1000 K, a perturbation amplitude = 5 × 10 −4 Ks −1 , a rotation period of 5 hours, surface gravity of 1000 ms −2 , and a storm timescale = 10 5 s. Results are taken after the models reaching statistical equilibrium. The instantaneous fields with drag = 10 5 s were obtained at about 800 simulation days; those with drag = 10 6 s were obtained at about 1000 simulation days; and those with drag = 10 7 s were obtained at about 2500 simulation days. models after the they fully equilibrate. The overall wind speed decreases significantly with increasing eff . Zonal jets show a speed up to ∼ 200 m s −1 in the model with eff = 500 K, and only up to several tens of m s −1 in the model with eff = 1000 K. In the model with eff = 1500 K, no zonal jet forms but only weak eddies with local maximum speeds up to a few tens of m s −1 are present at low latitudes. In the latter case, radiative damping sufficiently dissipates energy injected by the thermal perturbations before they can organize into zonal jets.
Lastly, we show a set of experiments with a fixed effective temperature, a weak drag, but varying amplitude of the thermal perturbation amplitudes. We choose eff = 500 K, drag = 10 7 s, and = 1.5 × 10 −4 , 6 × 10 −5 , 3 × 10 −5 and 1.5 × 10 −5 Ks −1 . All models exhibit robust zonal jets due to the weak drag and low tem-perature. These jets also show strong pressure-independent components at off-equatorial regions and vertical shears at low latitudes. Figure 5 shows the mass-weighted, time-averaged zonal-mean zonal wind profiles for these models, and these profiles are shifted with dashed lines corresponding to their zero values. The zonal jet speed decreases and the number of zonal jets increases with decreasing : there are two eastward jets with speed up to more than 350 m s −1 in each hemisphere in the model with = 1.5 × 10 −4 Ks −1 , but 4 to 5 eastward jets with speed up to 100 m s −1 in each hemisphere with = 1.5 × 10 −5 Ks −1 . These demonstrate that, with the same drag timescale and radiative damping rate, the equilibrium between energy injection and frictional dissipation is obtained at a lower jet speed when the energy injection rate is low. Classical two-dimensional turbulence theory predicts that the meridional jet  width is correlated to an anisotropic lengthscale, the Rhines scale, Rhines 1975, and see a review by Vasavada & Showman 2005). The number of jets on a sphere can be expected as jet ∼ / ∼ √︁ Ω / , where ∼ 2Ω/ and is the radius. With a rotation period of 5 hours and = 6 × 10 7 m, the characteristic number of jets is expected to be about 7 for a wind speed = 400 m s −1 and 14 for a wind speed = 100 m s −1 . Although the trend that higher wind speed corresponds to less number of jets is confirmed in our models, the number of jets in our models is obviously smaller than that predicted by the theory given the wind speeds. This is related to the vorticity structure of our simulated flows as will be explained below.
The vorticity structure provides additional information on jet formation. Absolute vorticity, defined as + , where = k · ∇ × v is the relative vorticity and k is the local upward vector on the sphere, is an approximation to the potential vorticity (PV), a quantity materi-(a) eff = 500 K. ally conserved under frictionless, adiabatic motions (Vallis 2006). In models with robust zonal jets, the absolute vorticity at a given level tends to be organized into zonal strips with large meridional gradient of absolute vorticity at the edges of adjacent strips but much smaller gradient within the strips. The relationship between the zonal-mean vorticity and zonal-mean zonal winds suggests that the boundaries between strips correspond to the eastward jets, whereas the interior of the strips correspond to the westward jets (e.g., McIntyre 1982;Dritschel & McIntyre 2008;Dunkerton & Scott 2008). If the absolute vorticity is well homogenized within the strip (the so-called absolute vorticity staircases), it naturally results in jets whose meridional spacing and jet speeds are related by the Rhines scale (e.g., Scott & Dritschel 2012). However, in our cases, the homogenization of absolute vorticity within the zonal strips is not perfect, and the corre-

Latitude [deg]
Barotropic zonal-mean zonal wind Figure 5. Time-averaged, mass-weighted zonal-mean zonal winds as a function of latitude for models with an effective temperature of 500 K, a bottom drag timescale drag = 10 7 s and a rotation period of 5 hours. The thermal perturbation rates in these models (from left to right) are 1.5 × 10 −4 , 6 × 10 −5 , 3 × 10 −5 and 1.5 × 10 −5 Ks −1 , respectively. The jet profiles are shifted with dashed lines corresponding to their zero values. These results are obtained after the models reaching statistically equilibrium.
sponding jet structure is less sharp and the meridional jet widths are broader than those predicted from flows with perfect absolute vorticity staircases. This explains that the number of jets ( Figure 5) are less than that expected from the Rhines scale argument. The imperfect PV staircases seem to be the common outcome of the strongly forced and strongly dissipated systems (Scott & Dritschel 2012). In such situations, magnitude of the eddy absolute vorticity can be comparable to the difference between two adjacent zonal absolute vorticity strips and their random motions are able to smooth out the boundaries between the strips, preventing perfect staircases. Indeed, such significant turbulent, filamentary absolute vorticity structures due to turbulence and Rossby wave breaking are present in our models. Note that such zonal-mean PV structures in our models naturally lead to a non-violation of the barotropic stability criteria which is in contrast to jets in the Jovian atmosphere (Ingersoll et al. 2004). Rossby wave breaking and its eddy mixing on PV play a crucial role in jet formation of our models that are under isotropic forcing and damping. Rossby waves break more easily in regions of weaker meridional PV gradient, and the resulting turbulent mixing tends to homogenize the PV in the wave breaking regions and therefore further decreases the meridional PV gradient. This is a positive feedback mechanism to induce zonal jets. Given an initial meridional PV structure that has variations to the planetary vorticity , Rossby-wave breaking preferentially occurs at regions with weaker PV gradient, and the associated PV mixing weakens the meridional PV gradient even more, promoting further wave breaking in these regions. In the equilibrium state, zonal jets form as a result of PV zonal striping, and the balance in the jets is achieved in between the dissipation (mainly the frictional drag in some cases) and the PV mixing; for reviews of this mechanism, see Dritschel & McIntyre (2008) and . The diagnosis of our models and interpretations are qualitatively very similar to those shown in Showman et al. (2019), and we shall not repeat them here.
The formation of pressure-independent (the so-called barotropic) off-equatorial jets is likely a natural outcome of the tendency that kinetic energy is stored to the gravest scale in rapidly rotating systems. Such a tendency has also been observed in idealized terrestrial GCMs (e.g., Chemke & Kaspi 2015) and gas giant planet GCMs (e.g., Lian & Showman 2008;Schneider & Liu 2009;Young et al. 2019;Spiga et al. 2020). Early theoretical work using two-layer quasi-geostrophic (QG) models (Rhines 1977;Salmon 1980) and multiple-layer QG models (Smith & Vallis 2002) have shown that energy can be converted efficiently from the smaller-scale pressuredependent (the so-called baroclinic) modes to the barotropic mode. In our weak-drag models, the emergence of barotropic off-equatorial jets is not surprising as the QG approximation holds well at mid-tohigh latitudes. Baroclinic potential energy is injected and is partially converted to baroclinic kinetic energy; at least part of it subsequently cascades to the barotropic mode, driving the barotropic jets.
Another mechanism contributing to the formation of barotropic jets is the development of vertically in-phase jets by Coriolis force associated with the mean meridional circulation. This is the so-called downward control in Earth's stratosphere (Haynes et al. 1991) and has been proposed for the Jovian atmosphere (Showman et al. 2006). The principle is the following: under rotation-dominated conditions, the zonal acceleration of zonal winds confined within a certain pressure range causes a responsive mean meridional circulation with a Coriolis force counteracting the jet forcing. The overturning meridional circulation penetrates to layers without zonal forcing, causing a reverse meridional velocity. The corresponding Coriolis force drives the zonal winds to the same direction as the jets in the forced layers. Indeed, diagnosis (not shown) in our models suggest that zonal jets are forced by eddy convergence and balanced by Coriolis force in the perturbed layers, while in deeper layers, the Coriolis force associated with the mean meridional circulation helps to drive the deep jets to the same direction.
The deep layers are nearly convectively neutral ( Figure 1). The weak stratification greatly helps with driving the strong barotropic component of the off-equatorial jets by either the downward control mechanism (Showman et al. 2006) or turbulent energy transport (Smith & Vallis 2002). In general, in many GCMs for giant planets in which the deep layers tend to be neutral, the development of strong barotropic components of the off-equatorial zonal jets is very common (e.g., Lian & Showman 2008;Schneider & Liu 2009;Showman et al. 2019;Young et al. 2019;Spiga et al. 2020).
There is still a subtlety regarding the mean meridional temperature structure. If the atmosphere is externally forced by equator-to-pole irradiation difference as in the Earth, there will be a pressure-dependent jet component at mid latitudes regardless how strong the pressureindependent component is. In our cases, the imposed homogeneous temperature at the bottom boundary acts against this configuration. The question is, will dynamics organize itself to maintain a systematic meridional temperature variation (and thus a significant vertical jet shear at mid latitudes) against radiative damping? Such a picture was proposed by : when alternating jets develop, some waves dissipate preferentially at different latitudes. This drives overturning meridional circulation cells that advect high-entropy air down and low-entropy air upward, creating isobaric meridional temperature gradients. Indeed, there are obvious mean meridional temperature variations in some very weakly radiatively damped models in Showman et al. (2019) (see their Figure 2). In our models, the off-equatorial mean meridional temperature variation is much weaker than that of the local variations, presumably because of the strong radiative damping in conditions appropriate for BDs and isolated young EGPs.
Finally, we turn our attention to the equatorial jets. In almost all models including those with a strong bottom drag, there exist local equatorial eastward jets at altitudes well above the thermal perturbation regions (see Figure 3 for examples). These correspond to local superrotations, and they emerge from interactions between the vertically propagating equatorial waves and the mean flows. Diagnoses (not shown) of the cases in Figure 3 suggest that vertical eddy angular momentum transport exerts eastward torques in the superrotating re-gion, while horizontal eddy momentum transport counteracts it, i.e., 1 cos 2 ( cos 2 ) ∼ − ( ) where is latitude, is radius, and , and are deviations of zonal, meridional and vertical velocities from the zonal average, respectively. Momentum transport by the mean flow also contributes to the balance. Jets in cases shown in Figure 3 are stable overtime. But in some other cases, especially those with a lower temperature, equatorial jets can evolve overtime in a quasi periodic manner, which will be described in next section.
However, all our models with relatively weak drags form equatorial westward jets at deep layers (see examples in Figure 3). Despite having local superrotations at low pressures, all these models exhibit equatorial sub-rotation when their jets are mass-weighted. Part of these models are shown in Figure 5. The common outcome with westward equatorial flow is similar to previous investigations using shallow-water models with small-scale random forcing (e.g., Scott & Polvani 2007;Showman 2007;Zhang & Showman 2014) and our previous study using the same forcing scheme but a different radiative damping scheme ). On the other hand, the strong equatorial eastward jets in the Jovian and Saturnian atmospheres penetrate to levels much deeper than their cloud decks.
As yet, there is no universal theory to predict whether westward or eastward equatorial flow should occur under conditions appropriate for gas giants. Multiple mechanisms have been proposed, but the resulting equatorial flows are sometimes sensitive to the detailed model setup (see a recent review by Showman et al. 2018). Regardless the details, in a shallow atmosphere framework (as opposed to the deep models that consider the fluid interior), the competition between the equatorially driven and off-equatorial driven forces largely determines the direction of the equatorial flow: (e.g., Schneider & Liu 2009;Liu & Schneider 2011;Showman et al. 2015). Schneider & Liu (2009) and Liu & Schneider (2011) argue that if the convective forcing is isotropic horizontally, one might expect equatorial superrotation. This is because convective heating fluctuations induce fluctuations in the large-scale horizontal divergence, and these represent a source of vorticity fluctuations and thus a source of Rossby waves. The Rossby wave source can be expected to be large in the equatorial region because the Rossby number there is order one or greater and large-scale horizontal flow fluctuations induced by convective heating fluctuations are divergent at leading order. On the other hand, the Rossby number at high latitudes is small, and largescale flows are nondivergent at leading order. Eddy fluxes due to waves propagating out of the equatorial region are expected to be more than those propagating into the equatorial region from mid latitudes. The former induces eastward acceleration while the latter exerts westward acceleration to the equatorial flow. Therefore, the net eddy angular momentum convergence tends to drive an equatorial superrotation.
Surprisingly, this picture does not describe our results despite that our forcing is globally isotropic. Even when our models do not impose the drag at the equatorial region as those in Schneider & Liu (2009), there is still no mass-weighted equatorial superrotation. One possibility may be that waves triggered at the equator remain a baroclinic structure. In the presence of a finite equatorial deformation radius for the baroclinic flow, these waves tend to be trapped in the equatorial region (Matsuno 1966) and few of them propagate out meridionally. Instead, they are transmitted in the vertical direction and drive vertically shearing equatorial jets. However, off-equatorial Rossby waves can freely propagate to the equatorial region, dissipating there and driving the equatorial westward flow. Indeed, our diagnoses suggest that the mass-weighted equatorial westward flow is maintained by both the horizontal eddy and mean-flow momentum transport while the frictional drag balances the westward forcing. It is interesting to note that even in the absence of meridional propagation of equatorial waves, based upon a shallow-water system, Warneford & Dellar (2017) show that dissipation of these waves can still lead to a net torque on the equatorial flow and the jet direction is sensitive to the form of the dissipation. It is yet unclear how these shallow-water theories can be transformed into understanding of the continuously stratified atmospheres.

Quasi-Periodic Oscillations of the Equatorial Jets
Quasi-periodic oscillations (QPOs) in the equatorial region have long been recorded in atmospheres of solar system planets. These include the quasi-biennial oscillation (QBO) in Earth's equatorial stratosphere (Baldwin et al. 2001), the Quasi Quadrennial Oscillation (QQO) with a period around 4 to 5 years in the Jovian stratosphere (Leovy et al. 1991;Friedson 1999;Cosentino et al. 2017;Antuñano et al. 2021), and the semiannual oscillation (SAO) with a period of ∼15 years in the stratosphere of Saturn (Fouchet et al. 2008;Guerlet et al. 2011). These QPOs are not related to their seasonal cycles but are driven by wave-mean-flow interactions involving the upward propagation of equatorial waves generated in the lower atmosphere and their damping and absorption in the middle and upper atmosphere (Baldwin et al. 2001;Guerlet et al. 2018). These interactions result in downward migration of the vertically stacked, alternating eastward and westward equatorial jets overtime. The QPOs viewed at certain levels are simply a manifest of the jet migration cycles. Our precursor work has shown that QPOs may be a common outcome in atmospheric conditions appropriate for some BDs and isolated EGPs ) . Here we further explore properties of the QPOs in a somewhat strong forcing and damping regime over a wider parameter space using the updated GCM.
We first conduct experiments over a grid with three effective temperatures of eff = 500, 1000 and 1500 K and three storm timescales = 10 4 , 10 5 and 10 6 s. Note that characterizes the typical evolution timescale of the injected thermal perturbation pattern. The former affects the radiative damping on dynamics and the latter affects the wave properties. A bottom drag with drag = 10 5 s is applied in all models. The strong drag suppresses formation of strong jets such that the wave properties are not influenced by the deep jets. The thermal forcing amplitudes are = 1.5 × 10 −4 Ks −1 for eff = 500 K, 5×10 −4 Ks −1 for eff = 1000 K and 2×10 −3 Ks −1 for eff = 1500 K. The typical equilibrium isobaric temperature variations near the RCBs are between 60∼100 K in these models. Figure 6 shows the equatorial zonal-mean zonal wind as a function of time and pressure of this grid. All the time sequences start after the models reaching statistical equilibrium. QPOs are present in some cases. We take the case with eff = 500 K and = 10 5 s as an example to illustrate its basic property. The downward migration of jets in this case has a mean period about 200 (Earth) days. The jet downward migration starts near the upper boundary and ends slightly above 1 bar. The migrating eastward and westward jets have similar magnitudes, and have strong vertical wind shear between them. The QPO is horizontally confined within a narrow latitudinal range (± ∼ 10 • ) around the equator.
The important overall picture is that the QPOs occur preferentially in the low left corner of the grid-with low atmospheric temperatures and short storm timescales. In models with high atmospheric temperatures eff = 1500 K and models with long storm timescale = 10 6 s, the QPO never occurs. In models with eff = 1000 K, only the one with = 10 4 s has a QPO. Within the same atmospheric s =10 6 s, T eff =500 K  Figure 6. Zonal-mean zonal jet at the equator as a function of pressure and time for models in a grid with three effective temperatures of 500, 1000 and 1500 K (from the left to right) and three storm timescales = 10 6 , 10 5 and 10 4 s from the top to the bottom. The bottom drag timescale is drag = 10 5 s for all models. The time zero starts where the models have reached statistical equilibrium. temperature, the shorter the storm timescale, the faster the jets migrate, a phenomenon clearly illustrated in the column with eff = 500 K. Given a , the cooler the atmosphere, the faster the jets migrate as shown in the row of = 10 4 s (even considering that the hotter models have higher forcing rates). The jet migration in models with eff = 500 K is coherent, but that in the model with eff = 1000 K and = 10 4 s has a more complicated pattern, showing independent migration of two eastward jets in a complete cycle.
We now briefly explore the dependence of the QPOs on other parameters. A similar set of models are performed with eff = 500 K, = 10 4 , 10 5 , 10 6 s, drag = 10 5 s, but with thermal perturbation amplitude ten times smaller than those in the left column of Figure  6. Reducing simply results in less wave flux. Only the model with = 10 4 s exhibits a QPO, and it has a slower equatorial zonal wind speed, a longer oscillation period and a lower pressure where the jet migration ends. Next we conduct a similar set of models as those in the left column of Figure 6 with eff = 500 K, = 10 4 , 10 5 , 10 6 s but with a longer drag timescale drag = 10 7 s that promotes formation of strong deep equatorial jets. At the jets develop, no model can sustain the QPO. However, during the spin-up phase wherein the deep equatorial jet is still weak, the QPO can occur in the model with = 10 4 s. This is illustrated in Figure 7. In the first few hundred simulation days, the QPO appears in a similar fashion as the strong-drag models. As the deep jet grows increasingly westerly, the QPO period increases and the pressure at which the jet migration ends decreases over time. Finally after about 1500 days, the QPO is completely ceased when the deep westward jet approaches a steady state with a speed ∼ 300 m s −1 .
To qualitatively understand varieties of the QPO with different parameters, we briefly recap the fundamental dynamical mechanisms developed to understand the QBO on Earth. The classic framework developed by Lindzen & Holton (1968), Holton & Lindzen (1972) and Plumb (1977) considers a 1-D height-dependent model for equatorial zonal-mean zonal winds with parameterized wave-mean-flow interactions. Waves propagate upward through a background with vertical jet shear, and preferentially dissipate and break near their critical levels (the level where the jet speed is equal to wave zonal phase speeds). Angular momentum of waves is deposited to the jet at where they dissipate. In the absence of large viscosity and counteracting force, the jet is accelerated. Because the peak zonal acceleration of the jet occurs at altitudes below the jet peak, the jet migrates downward over time. As this occurs, the wave-dissipation level moves downward further accompanied with the migration of jet, and this feedback ensures a continuous jet migration. The essence of the QBO can be understood with a simple and elegant framework that considers two discrete upward propagating gravity waves with identical amplitudes and equal but opposite phase speed (Plumb 1977). Radiative damping is applied to wave amplitudes as they propagate vertically, and a minimal viscosity is applied to the jets. Suppose that initially a eastward jet is at a lower altitude than a westward jet, and the eastward propagating wave is trapped below the eastward jet, but the westward propagating wave can propagate freely until reaching the westward jet above. The two waves drive the downward migration of their corresponding jets. At some point a jet approaches the bottom and create a strong shear, which then leads to viscous destruction of that jet. The corresponding wave then freely propagate to higher levels, where wave damping gradually build up a new jet. The process repeats, leading to a periodic downward migration of alternating jets. The period of the oscillation crucially depends on the wave flux. Although latitudinal complications must be considered to explain some features of the QBO (Andrews et al. 1987), the essence of this simple picture generally applies. The rate of radiative damping may have two crucial effects in controlling the QPO behaviors in our models. First, radiative damping may affect the wave momentum flux that can reach the critical levels by damping the wave amplitudes during their propagation. Strong radiative damping with a rate approaching the Doppler-shifted wave frequency may reduce the sensitivity of wave dissipation to the background mean flow. 1 Radiative damping also influences the type of the waves that actually drive QPOs by preferentially damping out slowly propagating waves and preserving fast waves. Second, radiative damping serves as a source of diffusion on the zonal jets. To see this, we write the thermal wind equation Eq. (5) in the equatorial plane (Holton 1973, Chapter 12.6 In the thermal-wind limit which holds reasonably well for zonalmean flows, significant vertical shear requires substantial meridional temperature differences. If the radiative damping is strong, this meridional temperature anomaly can be difficult to maintain over a timescale comparable to the jet migration timescale. In this case, the radiative damping will tend to trigger a residual circulation to oppose the driving force of the QPO. Putting the above arguments together, in hotter atmospheres, the eddy flux responsible for the QPOs may be sufficiently reduced, and the "radiative diffusion" on the zonal jets is so strong that it may statistically balance the vertical wave forcing. Thus, the equatorial regions in hot atmospheres can easily achieve steady states compared to cooler atmospheres regardless the storm timescale. If wave flux is small, a moderate rate of radiative damping can diffuse away the QPO forcing, which is why models with much smaller thermal perturbations amplitude lack a QPO even with a low temperature of eff = 500 K. The storm timescale likely has significant impacts on the wave properties. Despite that limited discretized modes are in the thermal perturbations, waves spanning a wide range of frequency and wavelength are triggered (see Figure 14 in Showman et al. 2019). This indicates that processes generating these waves include nonlinear effects. It may be possible that these processes with a shorter tend to trigger waves whose spectra spans a wider range of frequency and therefore the wave flux may be greater. On the other hand, the excited wave spectra may be more limited by a longer . A shorter is then expected to result in a shorter oscillation periods than those with a longer due to a larger excited wave flux, as shown in Figure  6. Quantitative demonstration of the above speculation is out of our scope and merits further studies.
The development of a strong jet near the wave-generation level affects wave propagation. Because perturbations are applied in the rotational frame, waves are triggered in the rotational frame but immediately experience the jet. When the jet is strong, waves with zonal phase velocities the same sign as the jet encounter a critical layer right after they form and never propagate out. This leaves only waves in the other phase direction propagating freely to the upper atmosphere. When the jet is not strong (compared to the zonal phase speeds of waves), the influence is weaker and does not totally suppress all waves in one direction. The strong deep jet breaks the east-west symmetry of wave forcing necessary for the QPOs. In Figure 7, the increasing QPO period suggests that wave fluxes decrease as the deep jet increases, indicating the increasing ability of the deep jet on blocking part of waves. As the deep jet approaches the same speed as that of the upper downward migrating westward jet, almost all west-propagating waves that are associated with driving the QPO are trapped, and the QPO is then ceased.
Finally we note that details of QPOs in our models are somewhat sensitive to vertical resolution, similar to Showman et al. (2019) and earth GCMs without gravity wave parameterization and tuning (e.g., Nissen et al. 2000). Our goal here is not to precisely predict properties of QPOs that might occur in real BDs and EGPs but to explore their first-order behaviours over a wide parameter space.

TRANSPORT OF PASSIVE TRACERS
In this section we show that the local overturning circulation generates vertical mixing of passive tracers and discuss the dependence of mixing on different tracer sources and sinks. Based on the GCM results, we estimate the vertical diffusion coefficients, zz , for different tracers if the global-mean vertical transport is approximated as diffusion. Finally we discuss how to analytically estimate the zz at conditions relevant for BDs and isolated EGPs without running a GCM. Figure 8. Colored panels: Time averaged, zonal-mean profiles of relative cloud mixing ratio (in logarithm scale) as a function of latitude and pressure for different particle sizes. Only mixing ratio larger than 10 −8 is plotted, and the blank regions have mixing ratio less than 10 −8 . The black dash lines are the condensation level, below which the tracers are in the form of vapor. The model has an effective temperature of 1000 K, a rotation period 15 hours, a surface gravity 1000 ms −2 , a storm timescale = 10 5 s, a drag timescale drag = 10 5 s, and a forcing amplitude 5 × 10 −4 K s −1 . Lower right panel: time and globally averaged relative cloud mixing ratios as a function of pressure for particles with different size. Solid lines with different colors (from black to blue) correspond to particles with different size from 0.05 to 5 m.

Cloud tracers
Cloud particles can be vertically transported in the stratified atmosphere against sedimentation via local overturning circulations driven by the thermal perturbations. We present representative results using models with eff = 1000 K. Models with different eff but the same forcing pattern show qualitatively similar results because the settling speeds depend primarily on assumed cloud particle size and gravity but only moderately on temperature. In a primary model, the following setup is used: a surface gravity of 1000 ms −2 , a storm timescale = 10 5 s, a drag timescale drag = 10 5 s, and a forcing amplitude = 5 × 10 −4 Ks −1 . However, here we use a longer rotation period of 15 hours and a higher horizontal resolution (C128) to better resolve the fine structures within the local Rossby deformation radius. This is because we found extra sensitivity of tracer transport to the horizontal resolution compared to jet formation, and very fine structures within a deformation radius are needed to be resolved in order to achieve convergent results in terms of tracer transport. In the following models, the condensation pressure cond is set at 10 bars, below which the tracers are in a vapor form and above which they are subjected to sedimentation. We examine particle radius of = 0.05, 0.1, 0.2, 0.5, 1.0, 2.0 and 5.0 m. With a given circulation, small particles are easily mixed up to high altitudes relative to the condensation level due to low sedimentation speeds, while large particles sink down quickly and form vertically compact cloud layers. This is shown in Figure 8, in which the colored contour panels show time-averaged zonal-mean mixing ratio of cloud tracers as a function of latitude and pressure (on a logarithmic scale) for particles with different size. Dashed lines are the condensation level. The global-mean mixing ratio as a function of pressure for each particle size is in the lower right panel. The 0.05-m particle can extend roughly 4 to 5 scale heights from the condensation level, while the mixing ratio of 5-m particle decrease significantly within about one scale height above the condensation level. The cloud thickness smoothly decreases from small to large particle size. Near the cloud top, cloud abundance drops rapidly with increasing altitude due to the sensitivity of settling velocity to pressure at low pressures 2 . The zonal-mean mixing ratio of clouds exhibits interesting meridional 2 At low pressure the ratio of molecular mean free path to the size of cloud  variation of the cloud top, in which clouds are usually mixed higher at mid-high latitudes than at low latitudes. This may be due to the weaker local overturning circulation at low latitudes, which is related to the smaller temperature variation at low latitudes (see Figure 2). Pressures of global-mean cloud tops (here loosely defined as the pressure lower than which the global-mean tracer abundance is below particles is larger than unity, and the terminal velocity in this regime is inversely proportional to the pressure.  Figure 11. The mean meridional transport, mean vertical transport, eddy meridional transport and eddy vertical transport terms for the zonal-mean tracer budget (see Eq. 7) as a function of latitude at 1 bar for a cloud tracer with a size of 0.2 m. The result is time-averaged after the model reaching a statistical equilibrium.
1 × 10 −8 ) roughly corresponds to where the particle settling speed exceeds the vertical wind speed. Figure 9 shows the particle settling speeds as a function of pressure for different particle sizes. Profiles for the isobaric RMS vertical wind speed and the mean value of the top 0.1% vertical speeds are also included for a comparison. For small particles with sizes 0.2 m, conjunctions between the RMS vertical wind speed and the particle settling speeds are very reasonable to describe the cloud top pressures (see a comparison between Figure 9 the lower right panel in Figure 8). For larger particles with sizes 2 m, the the comparisons using the RMS vertical velocity obviously underestimate the cloud-top altitudes; and that in the case with 5 m even indicates that no clouds at all above 10 bars. Instead, the comparisons using the top 0.1% vertical speeds are more reasonable.
Cloud layers with smaller sizes are nearly horizontally homogeneous near the cloud base but show increasing horizontal inhomogeneity at altitudes closer to the cloud top. Figure 10 depicts snapshots of the cloud mixing ratio at different pressure levels for a particle size 0.1 m. At relative low pressures, the cloud pattern exhibits strong patchiness with typical shape and length scales strongly related to the regional-scale turbulence structure. However, larger particles ( 0.5 m) show patchiness even close to the condensation level due to their large settling speed. Despite that the degree of patchiness differs, the shapes of cloud patterns are consistent between different pressures, indicating vertically coherent tracer transport from the cloud base up to the cloud top. The patchiness naturally originates from the vertical coherence of the overturning circulation. The atmosphere at sufficiently low pressures is cloud-free because clouds must eventually settle out due to the increasing settling velocity with decreasing pressure. The overturning circulation transport the cloud-free air to high pressure with only moderate "contamination" from the surrounding cloudy area. If the cloud is thick as those with small sizes, the lateral mixing has sufficient time to reduce the horizontal variation between the cloud-free downdraft and cloudy updraft. Otherwise for those with large sizes, the inhomogeneity is easily achieved by the overturning circulation.
The cloud patterns are dominated by local-scale overturning circulations rather than the zonal-mean meridional circulation. This is easily visualized in Figure 10 and we quantify it here. The circulation can be decomposed to zonal-mean and eddy (relative to the zonal Figure 12. Results from the model with a weak drag of drag = 10 7 s and other parameters the same as that shown in Figure 8. Top three panels: timeaveraged zonal-mean tracer abundances for cloud particles with sizes of 0.05, 0.1 and 0.2 m. Bottom panel: instantaneous zonal-mean zonal wind as a function of latitude and pressure. These results are obtained after the model reaching a statistical equilibrium. mean) components, and the zonal-mean transport of a tracer is written in pressure coordinates as follows: where represents a zonal mean of quantity , is the planetary radius, is latitude, is the zonal-mean meridional velocity, the zonal-mean vertical velocity, and are the eddies, and represents the sinks. The vertical transport of clouds are primarily driven by the eddy vertical winds (represented by the term ( )/ ) rather than the zonal-mean circulation. Figure 11 shows the above several terms as a function of latitude at 1 bar for a cloud tracer with a size of 0.2 m. The result is time-averaged after the model reaching a statistical equilibrium. The vertical eddy transport dominate the overall local zonal-mean tracer transport, and the mean vertical and eddy horizontal transport play secondary roles. When we integrate over latitudes to obtain the global mean vertical transport, only the vertical eddy transport is essential, which is balanced by the gravitational settling. The vertical eddy transport show minimum near the equator, corresponding to the smaller tracers profiles near the equator.
We also perform a model with a weak-drag model of drag = 10 7 s and other parameters the same as the above. Figure 12 shows the time-averaged zonal-mean tracer mixing ratio for particles with sizes of 0.05, 0.1 and 0.2 m in the upper three panels, and the instantaneous zonal-mean zonal wind in the bottom panel. Strong westward jets develop near the equator and eastward jets are at high latitudes with speeds exceeding ∼ 300 ms −1 . The cloud thickness is roughly flat in between ∼ ±30 • latitudes and then decreases poleward. The latitudinal change of cloud thickness occurs at where the zonal-mean zonal wind changes from westward to eastward. This differs from the strong-drag case, in which cloud thickness monotonically increases poleward. A careful comparison between Figure 8 and 12 suggests that the overall vertical mixing is weaker in the weak-drag model compared to the strong-drag model, and is primarily contributed by the differences poleward of ∼ ±30 • . The strong barotropic midlatitude jets do not directly generate vertical tracer transport but instead weakens the horizontal temperature variation by increasing horizontal advection. The reduce of horizontal temperature variation then weakens the local overturning circulation and the vertical mixing of tracers. We also perform a set of similar experiments in models with a cooler effective temperature of eff = 500 K. Still, we see an overall weaker mixing in the weak drag model where strong zonal jets develop. This is not surprising because mixing is driven by the local overturning circulation. The strong jets tend to reduce the local temperature variation by advection, and these weak-drag models show less vigorous local eddy velocities and so less vertical mixing.
1D atmospheric models often treat global-mean tracer transport as a diffusion process with a prescribed profile of the diffusion coefficient zz (e.g., Ackerman & Marley 2001;Moses 2014;Charnay et al. 2018;Woitke et al. 2020). In the stratospheres of BDs and isolated young EGPs, the source of zz has been interpreted as small-scale internal gravity wave mixing (e.g., Freytag et al. 2010). In this study, we have shown that large-scale circulation serves as an alternative tracer transport mechanism. To help guiding 1D models on what reasonable zz profiles should be from large-scale dynamics, we may derive the vertical diffusion coefficients for tracers from our GCM results assuming that the global-mean transport is described as a diffusion process. The basic procedure is to diagnose the isobaric tracer flux in the GCM and relate that to the vertical mean tracer gradient (Parmentier et al. 2013;Zhang & Showman 2018a): where the brackets represent global and time averaging, = − log( / * ) is the log-pressure coordinate, is the scale height, * is a reference pressure and is the vertical velocity in -coordinates.
The diffusion coefficients derived from our GCM as a function of pressure are shown in the top panel of Figure 13 for different particle size. The profiles at pressures with tracer abundances < 10 −8 are truncated. First of all, the zz is different for different particle sizes. The zz is on the order of 10 2 m 2 s −1 near the cloud base and generally decreases with decreasing pressure; for small particles the zz is on the order of 10 m 2 s −1 at 1 bar. The decrease is likely due to the decreasing vertical velocity with increasing altitudes (Figure 9). The zz rapidly become very small when approaching the cloud-top pressures. At pressures much larger than the cloud-top pressures, the zz for different particle sizes seem to be independent of particle size. This is especially the case for particles with 0.05 to 0.5 m particles at between about 0.7 to 10 bars. The insensitivity of zz to particle size when settling speed is small compared to the flow vertical speed is similar to that found in hot-Jupiter GCMs (Parmentier et al. 2013;Komacek et al. 2019). The zz derived based on the global-mean tracer flux using Eq. (8) is about an order of magnitude smaller than that using a relation zz ≈ where is a scale height (the latter is shown in the bottom panel of Figure 13). The latter has been widely accepted in many 1D models or even interpretations of 3D models, but its validity has never been justified. We show that it leads to systematic overestimation of the zz compared to that based on true tracer fluxes. This agrees to the conclusion in Parmentier et al. (2013) and Komacek et al. (2019) for hot Jupiter atmospheres.
We also perform models with slightly different parameters than the one shown in Figure 13, and their derived zz profiles are shown in Figure 14. In two models, we use half and twice forcing amplitudes as the original model but with other parameters the same, and their zz profiles are shown in the left and middle panels of Figure 14. As expected, the zz generally increases with increasing forcing amplitudes, but the dependence of zz to forcing amplitudes seems to be slightly steeper than a linear relation. In another model, we set the condensation pressure cond to 3 bars (instead of 10 bars in the original model) with other parameters the same. The resulting zz profiles are shown in the right panel, and they generally agree well with those shown in Figure 13. Overall, these experiments show a common property that when particle settling speed is negligible compared to atmospheric vertical speed, the zz seems independent of particle size.

Chemical tracers
Chemical tracers with equilibrium profiles varying with pressure are subjected to vertical mixing and can potentially be in disequilibrium depending on the chemical timescales relative to the dynamical timescales. Using the strong-drag GCM with the same parameters as that for the passive clouds, we first show the mixing of hypothetical chemical tracers with several chemical timescales of 10 4 , 3.16 × 10 4 , 10 5 , 3.16 × 10 5 ,10 6 , 3.16 × 10 6 and 10 7 s that are independent of pressure. The numerical implementation of idealized chemical tracers is described in Section 2.2. The equilibrium profile is the same for all tracers and is shown as the dashed line in the top panel of Figure  15. The time-averaged global-mean tracer profiles as a function of pressure are also shown as solid lines. Tracers with chem 10 5 s are close to their prescribed equilibrium at almost all pressures. Tracers with chem 10 5 s exhibit obvious deviation to their equilibrium profiles. However, their global-mean abundances still decrease with decreasing pressure. This is somewhat different to chemical mixing using convection resolving models in Bordwell et al. (2018), in which the mean tracers are quenched, i.e., above a certain level the abundance is independent of height and is determined by the value near a level where the quenching occurs. The horizontal tracer patterns are considerably inhomogeneous, showing isobaric tracer variations the same magnitude as the isobaric mean values. The spacial distributions of the variations are similar to those shown for cloud tracers ( Figure 10) and are primarily driven by local overturning circulations. Similarly, these horizontal structures are coherent at different pressures.
We performed a similar exercise to diagnose the vertical diffusion coefficient zz for the global-mean chemical tracers. The bottom panel of Figure 15 shows the resulting GCM-derived zz profiles for these tracers. The zz profiles with different chemical timescales show little difference at large pressures near 10 bars but moderate differences up to a factor of several at pressures lower than a few bars. The zz is on the order of 10 2 m 2 s −1 near 10 bars then decreases with decreasing pressure (less than 10 m 2 s −1 at pressures < 1 bar) partly due to the decrease of vertical wind speed. There seems to be a trend that the zz is lower for tracers with shorter chemical timescales, and this is theoretically expected as will be shown below.
We also show a model with chemical timescales increases with decreasing pressure. The detailed implementation is described in Section 2.2. The chemical timescale at 10 bars is fixed at chem,bot = 10 3 s and that at the top boundary adopt several values as chem,top = 10 4 , 3.16 × 10 4 , 10 5 , 3.16 × 10 5 ,10 6 , 3.16 × 10 6 and 10 7 s. These tracers show much smaller deviations to their equilibrium profiles (not shown) compared to those with constant chemical timescales. The resulting zz profiles for the global-mean transport are shown in Figure 16 as a function of pressure. Their profiles are qualitative similar to but meanwhile are systematic lower than those with constant chemical timescales (bottom panel of Figure 15).

Analytic estimation of
Building upon Holton (1986), Zhang & Showman (2018a) presented a theory to derive the vertical diffusion coefficient if the global-mean tracer transport is described as a diffusion process. In a regime in which the tracers are relatively short-lived with equilibrium profiles depending on pressure along, the can be approximated as where W is the magnitude of vertical wind speed. is a horizontal mixing timescale characterized by ℎ / where ℎ is a characteristic horizontal length scale over which distinctive tracer columns are mixed. Komacek et al. (2019) subsequently derived similar formula for vertical mixing in hot Jupiters' atmospheres.
Here we apply this theory to explain the GCM-derived zz values, and provide estimates of the theoretical at conditions appropriate for BDs and young EGPs without running GCMs. Our goal is to  Figure 13. Middle: similar to the left panel but from a model with twice forcing amplitude as that shown in Figure 13. Right: GCM derived zz for particles with different sizes from a model with a condensation level at 3 bars but other parameters the same as that shown in Figure 13. find appropriate scaling relations for the vertical wind speed and horizontal mixing timescale given a set of planetary parameters. Eq. (9) can be applied to chemical tracers straightforwardly. However, the diffusion coefficient for cloud tracers, which differ to chemical tracers in the sink and source, was not discussed in Zhang & Showman (2018a) and Komacek et al. (2019). Now we derive an expression of zz for cloud tracers. Following Zhang & Showman (2018a), we write the tracer equation for clouds equivalent to Eq. (10) in Zhang where is the isobaric deviation to the global mean, term is the terminal velocity in coordinates and is the scale height. Multiplying Eq. (10) with and globally averaging, we obtain Our goal is to derive a relation of to the vertical mean tracer gradient, and then can be estimated without knowing the distribution of tracers a priori. When term W, the last term in the left-hand side of Eq. (11) is negligible and the relation of and is straightforward. When term is non negligible, clouds usually exhibit significant horizontal inhomogeneity due to the coherent vertical circulation. In this situation, we argue that ∼ W . 3 Using this relation, Eq. (11), the equivalent eddy diffusion approximation Eq. (8) and 2 ≈ W 2 , the for clouds is At levels where term W (and therefore term and term are also likely much smaller than −1 ), is independent of particle size. At levels where the settling speed approaches the vertical wind speed, decreases rapidly and we expect that vertical cloud transport is ceased above those levels. These properties agree well with the zz for clouds derived from the GCM. Now we seek scalings for the vertical wind speed W and horizontal mixing timescale . Vertical mixing at mid-to-high latitudes is representative for cases with either strong or weak drags, therefore we adopt the quasi-geostrophic dynamical properties there. Thermal wind balance Eq. (5) provides a reasonable estimate of the horizontal eddy wind speed, which reads in order of magnitude ∼ where Δ is a characteristic isobaric temperature variation and ∼ Ω. Temperature perturbations generated by convective penetration are expected to have vertical wavelengths on the order of the scale height, therefore eddies driven by these perturbation have ln ∼ 1. From the GCMs, we note that the eddy horizontal structures are mostly around the size of a deformation radius, and therefore ℎ can be approximated as ℎ ∼ /Ω. Then, the characteristic horizontal mixing timescale is Applying the typical conditions of Δ ∼ 60 to 100 K from our GCMs, we obtain of a few hundred m s −1 , consistent with horizontal eddy wind speeds at mid latitudes (see Figure 2). The resulting horizontal mixing timescale is typically on the order of 10 4 s. For general conditions of BDs and isolated young EGPs, the balance between convective perturbation and radiative damping may results in isobaric temperature variation Δ between 50 to 200 K (see the detailed argument in the appendix of Showman et al. 2019), which is the typical range adopted in our analytic estimation.
Vertical motions are a result of horizontal convergence, and this can be decomposed into that by the ageostrophic motions and the geostrophic motions. The former scales as ageo / ∼ W/ and the latter scales as geo / ∼ geo / , where is the planetary radius . We may place an upper limit on the characteristic vertical wind speed due to convergence of the agostrophic wind using W/ ∼ / , where = / is the Rossby number and is much less than 1 at mid latitudes. The convergence due to ageostrophic motions is typically much greater than that by geostrophic motions in condition relevant here, and the scacling for the vertical velocity can be written as where we argue again that the relevant horizontal lengthscale is the global mixing. As the correlation between vertical velocity and abundance anomaly is expected to be positive, one type of columns is upwelling with a vertical velocity W and an eddy abundance + and the other is downwelling with a vertical velocity −W and an eddy abundance − + . The average flux is = ( W + − W (− + ) )/2. As clouds are quite inhomogeneous, we assume + ∼ , and thus, we expect that ∼ W .  Figure 17. Analytically derived zz profiles for clouds with different particle sizes using Eqs. (12), (13) and (14), and adopting Δ = 80 K, Ω = 1.16 × 10 −4 s −1 and a mean temperature of 1100 K. The settling speed is calculated using a T-P profile of eff = 1000 K.
deformation radius and ∼ Ω. It can be shown that where = / and = 1 for an isothermal atmosphere and = 0 for an adiabatic atmosphere. For conciseness, we use a fixed value for for a given T-P profile rather than evaluating locally at each pressure levels. This is because eddies span certain vertical depths and a somewhat averaged is more appropriate. Here we adopted a value = 0.8.
With the above scalings of Eqs. (13) and (14), settling speed of particles given a T-P profile and the known chemical timescales, we are in the position to estimate the theoretical for both chemical tracers using Eq. (9) and cloud tracers using Eq. (12). Note that we do not intend to represent variations of vertical wind speeds and horizontal mixing as a function of height, and the resulting vertical dependence of analytic zz comes from the variation of tracer sink and source. Here we apply Δ = 80 K, Ω = 1.16 × 10 −4 s −1 and a mean temperature of 1100 K relevant to GCMs in this section. For chemical tracers with local chemical timescales of 10 3 , 10 4 , 10 5 , 10 6 and 10 7 s, the analytic is about 2, 17, 49, 61 and 62 m 2 s −1 , respectively. The analytic captures the order of magnitude behavior of the GCM results but cannot capture the pressure dependent variation of the GCM-derived even in cases with constant chem . The strong dependence of zz on chem suggested that the variation of GCM-derived zz shown in Figures  15 and 16 is contributed by the variation of chem .
For cloud tracers, Figure 17 shows the analytic profiles as a function of pressure for particle sizes of 0.05, 0.1, 0.2, 0.5 and 1 m. Clouds with sizes of 2 and 5 m have settling speeds greater than the estimated vertical wind speed, and therefore cannot be applied. For small particles, the zz is on the order between 10 to 10 2 m 2 s −1 at high pressures but rapidly decreases when it approaches cloud top levels. The analytic profiles again capture the overall behavior as well, including the agreement on the order of magnitude value and the behavior that the zz approaches very small values near the level where setting speed roughly match the mean vertical wind speed. The theory of by global circulation takes vast idealization of the actual dynamical system and its inaccuracy when comparing to GCM results is not surprising (Zhang & Showman 2018a). Nevertheless, the theory captures some key properties and the order-of-magnitude estimation of as a function of various physical parameters and provides a guidance for parameterization of mixing for 1D models.

Effects of dynamics on the synthetic lightcurves
Recent long-term lightcurve observations of several BDs provide a valuable window to probe jets, waves, large-scale vortices and their time evolution in these atmospheres (Apai et al. 2017;Apai et al. 2021). The injected temperature perturbation patterns by themselves alone produce lightcurve variability even without dynamics. The question is then how does dynamics shape the lightcurve on top of the forcing pattern? Here we show simulated lightcurves viewed at the equator for models with a 5-hour rotation period, eff = 1000 K, with drag = 10 7 , 10 5 s and without dynamics in the top panel of Figure 18. The rotational modulations are obvious in all cases, showing sinusoidal-like shapes with a period corresponding to rotation period of 5 hours. The variability amplitudes are small (∼ 1%), and those with dynamics show lower amplitudes than that without dynamics. Long-time evolution of the lightcurves results from temporal variability in the temperature spatial patterns. We perform periodogram to the lightcurves that cover over 100 rotation periods and their nomalized power spectrum are shown in the middle and bottom panels of Figure 18. Properties of the lightcurve from the strong-drag model show small difference to that without dynamics. This is because with a strong drag, temperature patterns are close to the forcing patterns. The spectral power of the weak-drag model is overall quite similar but shows some quantitative differences to that without dynamics: the relative power of certain peaks near the rotation period is different, and perhaps more notably, the peak at 2.5 hour is exceptionally strong. The weak-drag model develops a global-scale, zonally propagating Rossby wave in each hemisphere and they induce certain temperature variations. The two waves have roughly the same zonal phase velocities but their zonal phases differ. Their effect on the lightcurve is as if there is a feature with zonal wavenumber 2, which is likely responsible for the peak at 2.5 hour in the power spectrum.

Connections to the interior and constraints on the jets
Large-scale winds are weak in the metallic interiors of BDs and young EGPs where MHD effects is important to break large-scale zonal jets, which have been suggested by theoretical arguments and numerical simulations of convection in rotating, electrically conducting spherical shells (e.g., Busse 2002;Liu et al. 2008). Application to typical BDs showing that the typical wind speed are of ∼ 0.1 − 0.3 m s −1 in the metallic regions ). In the molecular envelope overlying the metallic region, the thermal wind balance likely holds well for rapidly rotating BDs and EGPs. Due to the rapid rotation, zonal flows tend to be organized into columns parallel to the rotation axis (Busse 1976). As winds are weak at the bottom of the molecular envelope, the characteristic zonal jet speeds near the top of the convective zone depend on the characteristic large-scale horizontal density differences along the rotation axis. Scaling theories and numerical simulations suggest that this density difference is small in the fully convective region (e.g., Stevenson 1979;Showman et al. 2011).  suggest that the typical offequatorial zonal wind speed near the top of the convective zone is expected to be a few m s −1 (if the jet width is close to the planetary  On the other hand, the strong off-equatorial zonal jets of Jupiter likely penetrate to a significant depth in the molecular envelope . Recent studies have suggested that the presence of stably stratified layers above the metallic region may explain the deep penetration of the off-equatorial zonal jets (e.g., Christensen et al. 2020;Gastine & Wicht 2021). These stratification may be a result of compositional separation of hydrogen and helium. For BDs and hot EGPs, their interior may be too hot for helium-hydrogen separation to occur in their interior (Guillot 2005;. If the the interior is fully convective (and therefore is closed to be barotropic), theoretically we would expect that off-equatorial jets near the top of convective zone would be at most a few tens of m s −1 . Other possible mechanisms generating stratification or breaking the thermal wind balance in the interior may potentially alter this expectation. Given that the off-equatorial jet structures in our shallow models are largely barotropic, the above consideration may set a somewhat strong dissipation to the zonal jets in the observable layers of BDs and young EGPs generated by shallow mechanisms.
The equatorial jet needs special consideration. Near the equator, cylinders concentric with the rotation axis do not intersect the layer with MHD drag, and the equatorial jet may not suffer from the drag (Schneider & Liu 2009). The latitudinal region of this "drag-free" zone depends on the depth over which the MHD drag becomes sig-nificant. The depth may be very shallow for BDs and young EGPs as their interiors are hot and thermal ionization could provide significant conductivity even before hydrogen becomes metallic. To crudely illustrate this, we consider the ionization of potassium alone and estimate effects of the Lorentz drag on the horizontal winds in the kinematic limit as in Perna et al. (2010). Assuming a solar metallicity and a surface magnetic field strength of 1 kG (which may be appropriate for some field BDs, e.g., Kao et al. 2018), for three models shown in Figure 1, the maximum characteristic drag timescales at 100 bars are on the order of 10 −1 , 10 0 and 5×10 3 s. For a Jupiter-like field strength of 3 G, the drag timescales are about 5 × 10 3 , 2 × 10 5 and 6 × 10 8 s. These estimates indicate that, with a magnetic field strength reasonable for these objects, the MHD effects extend significantly to near the surface for L and T dwarfs. Assuming the surface is at about 1 bar, the depth at, for example, 100 bars, is on the order of 5 × 10 −4 of the radius assuming a Jovian size, leading to the meridional width of the equatorial "drag-free" zone of about ±2 • within the equator. Using the width of the drag-free equatorial zone and the argument of absolute vorticity homogenization (which assumes that jet is shallow and cannot violate the barotropic stability), one may obtain an estimate of the equatorial jet speed (Schneider & Liu 2009). Applying this argument, we have 2 /2 where is the distance between the equator and the latitude where the equatorial jet profile drops to zero. For a 5-hour rotating object with a Jovian size, = 1 × 10 −11 m −1 s −1 , and we obtain ≈ 24 to 500 m s −1 if assuming ∼ 2.2×10 6 to 10 7 m (a depth-to-radius ratio of 5×10 −4 to 5 × 10 −3 ). Of course, there is no guarantee that the jet profile will follow the absolute vorticity homogenization configuration-the jet speed can be well below these values if the absolute vorticity is less homogenized or exceed those if the jet violate the barotropic stability criteria.
Lightcurves in IR and radio wavelengths of a nearby T dwarf revealed evidence of fast eastward traveling features with a speed up to ∼ 600 m s −1 (Allers et al. 2020). Models in Tan & Showman (2021b) showed that zonal waves triggered by cloud radiative feedbacks may explain this net eastward signal. Other possibility includes advection by eastward zonal jets. Using both types of models and careful comparisons with observations, it may be possible to disentangle between the jet and wave mechanisms.

Other implications
Thermal perturbations injected by interactions between convection and the overlying stratified layer provide the minimal circulation in the stratified layers that easily advect sub-micron cloud particles several scale heights above the condensation level. Cloud decks of submicron particles are expected to be homogeneous near the cloud base and inhomogeneous towards the cloud top. However, the cloud decks of larger particles are vertically thin and significantly inhomogeneous even in the cloud base. L dwarf atmospheres are characterized by thick clouds composed of mostly sub-micron particles (Marocco et al. 2014;Hiranaka et al. 2016;Burningham et al. 2017), while T dwarf atmospheres are characterized by large particles and vertically thin layers (e.g., Allard et al. 2001;Burrows et al. 2006;Saumon & Marley 2008;Charnay et al. 2018). The emergence of significant cloud patchiness and thinning has been proposed to explain the rapid L to T transition (Knapp et al. 2004;Marley et al. 2010), and our circulation models and cloud transport provide a natural dynamical mechanism as to why these cloud behaviors would be expected when cloud particles become large. Chemical disequilibrium has been inferred in many field L and T dwarfs. The needed vertical diffusion coefficient zz of chemistry mixing in 1D models to match the ob-served spectrum spans a wide range between 1 to 10 4 m 2 s −1 (10 4 to 10 8 cm 2 s −1 ), and most models need the mixing above the convective layers (e.g., Saumon et al. 2006Saumon et al. , 2007Stephens et al. 2009;Moses et al. 2016;Leggett et al. 2016Leggett et al. , 2017Miles et al. 2020). The lower end of this range of zz seems much smaller than that derived from direct convective mixing (e.g., Ackerman & Marley 2001), but agrees well to our expectation from large-scale transport in the stratified layers driven by convective perturbations. In the convection resolving models by Freytag et al. (2010), the zz of cloud tracers can reach up to ∼ 10 3 m 2 s −1 in the stratified regions by internal gravity waves, providing the high-end mixing in the stratified layers needed to interpret observations.

Summary
Using an updated general circulation model with a grey radiative transfer scheme, we have investigated the formation and evolution of jet streams and tracer mixing in atmospheres of brown dwarfs and distant young giant planets under globally isotropic thermal perturbations which are thought to be caused by interactions between convection and the overlying stratified atmosphere . We have reached the following conclusions: • Under relatively weak damping conditions, robust zonal jet streams with speeds up to several hundred m s −1 are a natural outcome of interactions of waves and turbulence with the planetary rotation. At a particular atmospheric temperature with a reasonable thermal perturbation rate, weak bottom frictional drags promote jet formation at all latitudes, while strong drags leave isotropic turbulence at mid-to-high latitudes and jet only at the equator and at altitudes above the drag domain. With a fixed thermal perturbation rate and weak drags, low atmospheric temperature (representing weak radiative damping) promotes an overall stronger circulation and zonal jets, while high temperature (strong radiative damping) leads to a much weaker circulation. With both weak frictional and radiative damping, the jet speed increases and the number of jets decreases with increasing thermal perturbation rates. The off-equatorial jets are largely pressure independent but the equatorial jets usually exhibit significant vertical jet shears. There is negligible systematic meridional temperature variation associated with the off-equatorial jet structure.
• Similar to Showman et al. (2019), our models exhibit long-term oscillations in the equatorial zonal jets, consisting of stacked eastward and westward zonal jets that migrate downward over timescales ranging from tens of days to months. The mechanism responsible for the quasi-periodic oscillations (QPOs) is similar to those in stratospheres of Earth, Jupiter and Saturn. Conditions with moderate to low atmospheric temperatures and strong thermal perturbations favor the formation of QPOs, while high temperature (strong radiative damping) inhibit QPOs.
• Circulation driven by the thermal perturbations leads to efficient vertical mixing of tracers including clouds and chemical species in the stratified atmosphere. The mixing is mostly due to local-scale overturning circulation rather than the global-scale mean meridional circulation. Sub-micron cloud particles are easily mixed up to several scale heights above the condensation levels, showing horizontal homogenization near the cloud base but patchiness toward the cloud top. Larger cloud particles show vertically thin layers with significant patchiness even near the cloud base. Chemical tracers with relatively long chemical timescales (> 10 5 s) are driven out of equilibrium throughout the atmosphere. If the global-mean tracer transport is characterized by diffusion, the vertical diffusion coefficient, zz , de-rived from the GCM results is on the order of 1 to 10 2 m 2 s −1 (10 4 to 10 6 cm 2 s −1 ), and its profile differs for different types of tracers. We have derived analytic scaling relations to estimate zz for different types of tracers appropriate for atmospheres of BDs and young isolated EGPs. Our analytic zz profiles qualitatively capture some key features of and on order of magnitude agree well to those derived from GCM results.