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

Xianyu Tan1 and Adam P. Showman2,3,†
1Atmospheric Oceanic and Planetary Physics, Department of Physics, University
of Oxford, OX1 3PU, UK
2Lunar and Planetary Laboratory, University of Arizona, 1629 University
Boulevard, Tucson, AZ 85721, USA
3Department of Atmospheric and Oceanic Sciences, Peking University, Beijing,
People’s Republic of China
${}^{\dagger}{\rm Deceased.}$ E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

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

###### keywords:

hydrodynamics — methods: numerical — planets and satellites: atmospheres —
planets and satellites: gaseous planets — brown dwarfs

††pubyear: 2021††pagerange: Atmospheric circulation of brown dwarfs and
directly imaged exoplanets driven by cloud radiative feedback: global and
equatorial dynamics–B.3

##  1 Introduction

Brown dwarfs (BDs) are objects with masses intermediate between stars and
giant planets, but are analogous to giant planets in temperature, composition
and size (Burrows et al., 2001). Isolated BDs are easier to observe than
exoplanets orbiting bright stars, making them ideal targets to investigate
physical, chemical and dynamical processes in the context of planetary
atmospheres (Apai et al., 2017). Growing observations of BDs in the past
decade have revealed that active atmospheric circulation is common among BDs.
Here we summarize several key types of observations that set direct
constraints on the large-scale circulation:

  * •

Broadband lightcurve time variability at infrared (IR) wavelengths is commonly
observed for field BDs (e.g., Gelino et al., 2002; Clarke et al., 2008;
Artigau et al., 2009; Gillon et al., 2013; Buenzli et al., 2014; Wilson et
al., 2014; Metchev et al., 2015; Leggett et al., 2016; Miles-Páez et al.,
2017; Apai et al., 2017; Manjavacas et al., 2017; Vos et al., 2019; Eriksson
et al., 2019; Bowler et al., 2020; Hitchcock et al., 2020; Vos et al., 2020,
see also recent reviews by Biller, 2017 and Artigau, 2018). Most of them are
likely caused by the rotation of large-scale surface inhomogeneities related
to different cloud opacity and/or temperature that move in and out of view as
the objects rotate. The shapes of some lightcurves are irregular and change
over timescales comparable to or slightly longer than rotational timescales,
indicating rapid evolution of the large-scale spatial patchiness in some
cases. Apai et al. (2017) summarized a few types of irregularities of the
lightcurve variability and suggested a possible solution by invoking
differential zonal propagation of waves.

  * •

Multi-wavelength near-IR lightcurve variability provides additional
information on the vertical structure of the surface patchiness because
different wavelengths probe different atmospheric pressures (Buenzli et al.,
2012; Radigan et al., 2012; Biller et al., 2013; Apai et al., 2013; Yang et
al., 2015, 2016; Lew et al., 2016; Zhou et al., 2018, 2020; Lew et al., 2020).
The amplitude of the variability changes with wavelength, showing smaller
amplitude near the water absorption bands and larger amplitude near spectral
windows. These are typically consistent with the scenario that the surface
patchiness results from a combination of different cloud thicknesses rather
than a combination of a single cloud type and a completely clear-sky region
(Buenzli et al., 2012; Apai et al., 2013; Yang et al., 2015; Lew et al., 2016,
2020). Alternative scenarios invoking variation of temperature in clear-sky
atmospheres have also been proposed (Robinson & Marley, 2014; Morley et al.,
2014; Tremblin et al., 2020). The phases of lightcurve rotational variations
can be different at different wavelengths, and the extent of this difference
varies with spectral types (Yang et al., 2016).

  * •

Viewing-angle dependence of variability and near-IR color of field BDs has
been suggested by Vos et al. (2017), Vos et al. (2018) and Vos et al. (2020),
who showed that BDs viewed more equator-on tend to have higher variability
amplitudes and redder near-IR colors than those viewed more pole-on. These
results indicate the possible presence of systematic equator-to-pole
differences of cloud properties, temperature and/or chemical composition.

  * •

Doppler imaging could provide perhaps the most direct constraints on the
presence of large-scale surface patchiness and their morphology. This has been
done for Luhman 16B (the spectral type T1 component in the BD binary system
Lhuman 16AB) and yielded surface flux differences between the patchiness of a
few hundred K in terms of effective temperature (Crossfield et al., 2014).

  * •

Simultaneous tracking of near-IR and radio variability, in which the former
measures the period at which atmospheric features rotate in and out of the
view and the latter likely reflects the rotation period of the interior, could
track possible differential rotation between the atmosphere and the interior.
Recently, Allers et al. (2020) applied this technique to a nearby T6 brown
dwarf and showed that the period of the near-IR variability is slightly
shorter than that of the radio emission, suggesting that the dominant
atmospheric features travel eastward relative to the interior with a zonal
speed of a few hundred$\;{\rm m\;s^{-1}}$.

  * •

Net near-IR polarization could be due to scattering by cloud particles
together with either the presence of oblateness due to fast rotation or cloud
inhomogeneity (e.g., Sengupta & Marley, 2009; Marley & Sengupta, 2011).
Recently, Millar-Blanchaer et al. (2020) unambiguously detected net near-IR
polarization as well as its time variability on Luhman 16B. The time-mean
polarization indicates a latitudinally dependent cloud structure, and its time
variation indicates longitudinal cloud patchiness on top of the zonal cloud
structure.

Directly imaged extrasolar giant planets (EGPs) are mostly young, hot and
relatively distant from their host stars such that they receive negligible
external stellar irradiation compared to their interior heat flux. Their
spectrum and near-IR colors show similarities to the dusty field BDs (e.g.,
Currie et al., 2011; Barman et al., 2011; Marley et al., 2012; Rajan et al.,
2015; Chauvin et al., 2017; Stolker et al., 2020). Near-IR lightcurve
variability has been observed on directly imaged EGPs and planetary-mass,
free-floating objects (Biller et al., 2015; Zhou et al., 2016; Vos et al.,
2018; Biller et al., 2018; Schneider et al., 2018; Manjavacas et al., 2019;
Miles-Páez et al., 2019). From a meteorological point of view, directly imaged
EGPs resemble low-gravity versions of BDs and fall into the same category as
field BDs in terms of atmospheric dynamical regime.

These observations motivate the investigation of global atmospheric dynamics
of BDs and directly imaged EGPs, their fundamental properties and effects on
cloud formation and chemistry. There have been several studies investigating
atmospheric dynamics appropriate for these objects (Freytag et al., 2010;
Zhang & Showman, 2014; Tan & Showman, 2017; Showman et al., 2019; Tan &
Showman, 2020, and see recent reviews by Showman et al., 2020 and Zhang,
2020). Cloud radiative feedback has been proposed as a robust and novel
mechanism to drive spontaneous atmospheric variability and dynamics in these
atmospheres (Tan & Showman, 2019, 2020). This mechanism works as follows:
imagine two adjacent atmospheric columns, one of which has a thick cloud layer
and the other has a relatively thin cloud layer. The two columns experience
different radiative heating and cooling rates due to the different cloud
opacity, which then generate isobaric temperature differences. This drives an
overturning flow that maintains the patchy clouds against settling,
potentially sustaining the circulation. In a previous work, we confirmed the
viability of this mechanism and explored dynamical properties of the
circulation using local three-dimensional models that assume a Cartesian
geometry and a constant Coriolis parameter $f$ across the domain (Tan &
Showman, 2020). We demonstrated the importance of rotation in regulating the
large-scale atmospheric dynamics and cloud formation, showing that the typical
horizontal length scales of vortices and cloud patterns are inversely
proportional to the Coriolis parameter $f$ if other parameters are held fixed,
and that the mean cloud vertical extent decreases with decreasing rotation
period.

The latitudinal variation of the Coriolis parameter in global geometry, the
so-called $\beta$ effect where $\beta=df/dy$ and $y$ is distance increasing
northward, introduces additional dynamical behaviors compared to those in the
$f$-plane models (Tan & Showman, 2020). For example, the $\beta$ effect can
introduce horizontal anisotropy in the turbulence (Rhines, 1975; Vallis &
Maltrud, 1993) as well as large-scale atmospheric waves (Holton & Hakim,
2012), both of which play crucial roles in the global circulation of giant
planets (see reviews by Ingersoll et al., 2004; Vasavada & Showman, 2005; Del
Genio et al., 2009; Showman et al., 2018). In addition, the regional models in
Tan & Showman (2020) occupy only a limited fraction of the surface area of BDs
and EGPs, and therefore qualitative comparisons between model outputs and
observed lightcurve variability were lacking.

In this study, we extend the modeling framework of Tan & Showman (2020) to a
global geometry to investigate the global atmospheric circulation driven by
the cloud radiative feedback and the resulting lightcurve variability. We show
that when the bottom frictional dissipation is relatively strong, horizontally
isotropic storms and turbulence are prevalent at mid-to-high latitudes while
zonally propagating waves are present at low latitudes near the observable
layer. The differential propagation of equatorially trapped waves induces
short-term evolution of simulated lightcurves, analogous to the wave beating
effects described in Apai et al. (2017). Eastward propagating equatorial
Kelvin waves are sometimes dominant, causing a slightly shorter rotation
period of the atmospheric features relative to the underlying planetary
rotation period, which agrees well with observational results by Allers et al.
(2020). When the bottom dissipation is weak, strong and broad zonal jets
develop and modify wave propagation and lightcurve variability. We find
systematic equator-to-pole differences of clouds and temperature structures
due to the latitudinal variation of the Coriolis parameter $f$, supporting
recent observations by Vos et al. (2017) and Vos et al. (2020). Large-scale
equatorial disturbances may help to explain the nature of longitudinal
variation in Doppler mapping (Crossfield et al., 2014) and time-varying
polarization (Millar-Blanchaer et al., 2020) of Luhman 16B.

The paper is organized as follows. Section 2 briefly introduces the global
numerical model. Section 3 describes results of models with different rotation
periods and drag timescales as well as their resulting lightcurve variability.
In Section 4 we discuss implications of the results and draw conclusions.

##  2 Model

The general circulation model (GCM) used in this study is fully described in
Tan & Showman (2020), and here we apply the GCM to a global domain. Key
elements of the model are summarized below. We solve the standard 3D
hydrostatic primitive equations using an atmospheric GCM, the MITgcm (Adcroft
et al., 2004). Two tracer equations which represent the mass mixing ratio of
condensable vapor and clouds are integrated simultaneously. We assume the
ideal gas law for the equation of state. The deep layers in our models have
reached the convective zone which is typically at pressures larger than a few
bars for typical L and T dwarfs (Burrows et al., 2001). Effects of rapid
convective mixing on both entropy and tracers are parameterized using a simple
convective adjustment scheme that instantaneously homogenizes entropy and
tracers in the convectively unstable regions.

A Rayleigh drag is applied to the horizontal winds at pressures higher than 5
bars to crudely represent interactions between the modeled atmosphere and the
quiescent interior where large-scale flows are likely retarded due to
significant magnetohydrodynamic dissipation. Rapid convective mixing of
specific entropy and fast rotation may lead to the Taylor–Proudman effect that
could be efficient in slowing down large-scale winds in the shallow outer
layer (see the detailed argument in Schneider & Liu, 2009). Because of the
much hotter interior of L and T dwarfs and likely strong magnetic fields of
BDs ($\sim$kG, Kao et al., 2016, 2018), the “quiescent" region inside BDs is
expected to extend to a much shallower depth than that of Jupiter. Therefore
the drag is applied uniformly in the horizontal direction. The form and
strength of the drag are highly uncertain because of the unknown nature of
interactions between the interior and the shallow outer layer. Yet no
theoretical study systematically investigate how such interactions should be
parameterized in GCMs of gaseous planets. Nevertheless, this simple drag
serves to dissipate flows in our bottom model layers, and similar drag schemes
have been used in previous studies of Jupiter (Schneider & Liu, 2009) and hot
Jupiters (e.g., Liu & Showman, 2013; Tan & Komacek, 2019; Carone et al.,
2020).

The drag is the strongest at the bottom boundary and is characterized by a
drag timescale $\tau_{\mathrm{drag}}$, then it linearly decreases with
decreasing pressure until it reaches zero at 5 bars. There is no drag at
pressures lower than 5 bars. As in Tan & Showman (2020), we set a relatively
strong drag timescale $\tau_{\mathrm{drag}}=10^{5}$ s for the major suit of
models. We have tested models with stronger drags
($\tau_{\mathrm{drag}}=10^{4}$ and $10^{3}$ s), and they are quantitatively
similar to those with $\tau_{\mathrm{drag}}=10^{5}$ s. Therefore,
$\tau_{\mathrm{drag}}=10^{5}$ s is chosen to represent dynamics in the
“strong-drag" regime. For $\tau_{\mathrm{drag}}$ sufficiently larger than
$10^{5}$ s, dynamics is in the weak-drag regime. The drag is said to be
“strong" in a practical sense that the resulting zonal-mean zonal flows near
the bottom layer are rather weak (with velocities much smaller than $100\;{\rm
m\;s^{-1}}$) based on our modeling results discussed in Section 3.2, whereas
the drag is said to be “weak" when the resulting zonal flows near the bottom
layer are strong (with velocities much greater than $100\;{\rm m\;s^{-1}}$).
Using scaling argument in the convective layers, Showman & Kaspi (2013) argue
that the characteristic speeds of large-scale zonal flows near the top of
convective zone range from only a few to tens of $\;{\rm m\;s^{-1}}$ depending
on the number of jets. We therefore expect that the strong-drag regime might
be appropriate for BDs and hot directly imaged EGPs. Nevertheless, in section
3.2, we show results with weaker drags of $\tau_{\mathrm{drag}}=10^{6}$ and
$10^{7}$ s to explore possible circulation patterns for two reasons. First, it
is theoretically motivated to understand dynamics in different regimes in
spite of its applicability in realistic situations. Second, recent gravity
measurements of Jupiter and Saturn suggest that significant meridional density
gradients (and therefore the vertical wind shears associated with the zonal
jets) exist deep in the convective zone of Jupiter and Saturn (Kaspi et al.,
2018; Guillot et al., 2018; Iess et al., 2019). This indicates that strong
organized zonal flows near the top of convective zone remain possible for BDs.

Atmospheric motions are driven by the horizontal pressure gradient that is
rooted from horizontal differential radiative heating and cooling. In
calculating the radiative flux, we assume a gray atmosphere with a single
broad thermal band for simplicity and computational efficiency. The gas
opacity in our model atmosphere is
$\kappa_{\rm{gas}}=\max(\kappa_{\rm{R,g}},\kappa_{\rm{min}})$, where
$\kappa_{\rm{R,g}}$ is the Rosseland-mean opacity from Freedman et al. (2014)
assuming solar composition. The Rosseland-mean opacity gives a good estimation
of radiative flux in the optically thick limit. The atmosphere above about 1
bar is optically thin by the gas opacity alone, and there is no good choice a
priori for the opacity in the grey approximation. A minimal opacity
$\kappa_{\rm{min}}$ is then imposed in the gas opacity. We assume $\kappa_{\rm
min}=10^{-3}\;{\rm m^{2}\;kg^{-1}}$, which is consistent with that used in
semi-grey models in Guillot (2010) for hot giant planets. The total
atmospheric opacity $\kappa$ is simply the sum of the gas and cloud opacities
$\kappa=\kappa_{\rm{gas}}+\kappa_{c}$, the latter of which is determined by
the cloud mixing ratio as a function of time and location.

The cloud formation scheme is the same as that in Tan & Showman (2020). Cloud
forms when the mixing ratio of condensable vapor exceeds the prescribed
saturation mixing ratio $q_{s}$. Otherwise, evaporation occurs when the vapor
mixing ratio is less than $q_{s}$. The saturation mixing ratio $q_{s}$ is
assumed to depend on pressure alone: $q_{s}$ is $q_{\rm{deep}}$ at a
condensation pressure $p_{\rm cond}$ and then rapidly decreases when pressure
is less than $p_{\rm cond}$. At pressures larger than $p_{\rm cond}$, $q_{s}$
is assumed to be arbitrarily large such that no condensation would occur. The
condensation pressure $p_{\rm cond}$ is assumed to be 0.5 bar, roughly
consistent with conditions of typical mid-L dwarfs in which clouds form at and
above $\sim 1$ bar (e.g., Tsuji, 2002; Burrows et al., 2006). The deep mixing
ratio $q_{\rm{deep}}$ is assumed to be $2\times 10^{-4}~{}\rm{kg/kg}$. At
pressures higher than 5 bars, the vapor field is relaxed towards $q_{\rm
deep}$ over a timescale of $10^{3}$ s. The assumed $q_{\rm deep}$ in this
study is less than that of silicate vapor in atmospheres with solar abundance.
This value is tuned to generate cloud opacity that greatly exceeds the gas
opacity, but not to trigger convection within the cloud-forming region.
Convection within the cloud forming region complicates the dynamics as it
would introduce instability that does not totally rely on resolved large-scale
flow (Tan & Showman, 2019). We aim at exploring large-scale cloud-driven
dynamics in a cleanest possible context, therefore we do not chose a realistic
$q_{\rm deep}$.

We assume a constant cloud particle number per dry air mass $\mathcal{N}_{c}$
(in unit of $\rm{kg^{-1}}$) throughout the atmosphere. Cloud particles are
assumed to be spherical and have a single size locally in each grid cell, and
the particle size is then determined via a given $\mathcal{N}_{c}$ and time-
and location-dependent amount of condensate (using the scheme of Tan &
Showman, 2020). We assume $\mathcal{N}_{c}=10^{11}~{}\rm{kg^{-1}}$, which
results in typical cloud particle size around 0.5 $\mu$m given the assumed
$q_{\rm deep}$, consistent with that expected for L dwarfs (Saumon & Marley,
2008; Burningham et al., 2017). Optical properties of cloud particles,
including the extinction coefficient, scattering coefficient and asymmetry
parameter are averaged using the Rosseland-mean definition. Tables containing
these parameters as functions of temperature and pressure are pre-calculated
for the GCM to read in during the integration. We use enstatite (${\rm
MgSiO_{3}}$) to represent properties of the cloud species, including a density
$\rho_{c}=3190\;{\rm kg\;m^{-3}}$ and the refractive index. Enstatite is one
of the most representative cloud species in atmospheres of L and early T
dwarfs, and is appropriate for our atmospheric conditions (Ackerman & Marley,
2001; Saumon & Marley, 2008; Marley & Robinson, 2015). However, the choice of
cloud species is not critical in this study as long as the cloud opacity
exceeds the gas opacity.

The global models in this study have different horizontal resolution depending
on the rotation period—the shorter the rotation period, the smaller the
deformation radius and thus higher resolution is needed. Our horizontal
resolution of the cubed-sphere grid is C96 (equivalent to 384$\times$192 grid
points in longitude and latitude), C128 (512$\times$256), C192
(768$\times$384) and C256 (1024$\times$512) for rotation period of 20, 10, 5
and 2.5 hours, respectively. These horizontal resolutions are among the
highest published so far for exoplanet and brown dwarf global models111If the
horizontal resolution is insufficient to fully resolve dynamics within the
deformation radius of local regions (especially near the polar regions),
storms will cease to exist and there would not be dynamics in these regions..
The radius of the object is assumed to be $7\times 10^{7}$ m, similar to the
Jovian radius $R_{J}$. Although radius of real BDs varies from more than
$2R_{J}$ to slightly less than $R_{J}$ depending on their masses and ages, our
choice of radius should be representative for field BDs. A standard fourth-
order Shapiro filter is applied in the horizontal velocity and temperature
fields to maintain numerical stability (Shapiro, 1970). The pressure domain is
between 10 bars and $10^{-3}$ bar, which is discretized evenly into 55 layers
in log pressure. The model pressure domain is deep enough to reach the
convective zone, where we argue that rapid convective mixing leads to small
horizontal temperature gradient. The upper boundary is high enough for models
to safely capture dynamics associated with clouds. The temperature at the
bottom boundary (10 bars) is fixed at 2600 K, resulting in atmospheric
temperatures comparable to those of L dwarfs. We adopt physical parameters
relevant for BDs, including the specific heat $c_{p}=1.3\times
10^{4}~{}\rm{Jkg^{-1}K^{-1}}$, the specific gas constant
$R=3714~{}\rm{Jkg^{-1}K^{-1}}$, and a surface gravity $g=1000~{}\rm{ms^{-2}}$.
Models with $\tau_{\mathrm{drag}}=10^{5}$ s equilibrated after $\sim 100$
simulation days. After equilibration, we integrated them for additional 300
days, and the time-mean statistical results were obtained with outputs of the
last 200 days. Models with $\tau_{\mathrm{drag}}=10^{6}$ and $10^{7}$ s
equilibrated after $\sim 200$ and $\sim 1000$ simulation days, respectively,
and their statistical results were obtained similarly to those with
$\tau_{\mathrm{drag}}=10^{5}$ s.

##  3 Results

###  3.1 Results with varying rotation and strong bottom drag

####  3.1.1 Temperature and cloud patterns

Figure 1: Instantaneous global temperature maps at 0.23 bar on the left column
and the corresponding horizontal cloud mixing ratio at 0.23 bar on the right
column. These results are from models with four rotation periods of 20, 10, 5
and 2.5 hours (from the top to the bottom row) and with a drag timescale
$\tau_{\rm{drag}}=10^{5}$ s at the bottom boundary. Horizontal thin dotted
lines in the maps are constant-latitude lines with a spacing of $15^{\circ}$,
and vertical thin dotted lines are constant-longitude lines with a spacing of
$30^{\circ}$ Figure 2: Instantaneous global map of outgoing top-of-atmosphere
thermal flux corresponding to those shown in Figure 1.

Results of global models with rotation periods of 20, 10, 5 and 2.5 hours and
with a drag timescale $\tau_{\rm{drag}}=10^{5}$ s (note that the frictional
drag is applied only at pressures greater than 5 bars) are shown in Figure 1.
Figure 2 shows the corresponding outgoing top-of-atmosphere thermal flux.
These results are obtained after the models reach statistical equilibrium.
Pressures of thermally radiating levels are sensitive to cloud top
pressures—the radiating pressure could be lower than 0.05 bar in cloudy
regions and could be about 3 bars in clear-sky regions. Dynamics at 0.23 bar
is representative for layers near and slightly below the cloud forming region.
Quantities at slightly higher or lower pressures are qualitatively very
similar to those at 0.23 bar.

There is a lack of strong zonal jets, in the sense that the zonal-mean zonal
winds are much weaker than the eddy velocities, in all models shown in Figure
1 and 2 due to dissipation of kinetic energy by the strong bottom drag. The
circulation in the cloud forming region is dominated by active vortices,
turbulence and waves. At 0.23 bar, typical isobaric temperature differences
are above 100 K, and local horizontal wind speeds can exceed 1000$\;{\rm
m\;s^{-1}}$. Mid-to-high latitudes are filled with cyclones222Cyclones have
relative vorticity the same sign as the local planetary rotation, whereas
anticyclones have the opposite sign of relative vorticity. which are less
cloudy than their surroundings, and anticyclones which are associated with
thick clouds. In the cloudy areas, less thermal radiation escapes to space due
to large cloud opacity, and this cloud greenhouse effect warms up the lower
cloud forming region. In the cloud-free areas, more radiation can escape to
space from the hotter, deeper region due to the lack of cloud opacity,
efficiently cooling down the atmospheric column. This patchy-cloud
configuration is self-maintained by the circulation driven by clouds
themselves (Tan & Showman, 2020). Positive buoyancy is generated in the cloudy
regions, and it maintains the cloudiness against cloud gravitational settling
by transport of condensable vapor upward to the condensation level; whereas
negative buoyancy in the cloud-free regions generates downwelling that advects
cloud-free air from above, maintaining the IR cooling. In rapidly rotating
conditions, the tendency of geostrophy implies that the warm, cloudy regions
are anticyclonic and cool, cloud-free regions are cyclonic. The spatial
patterns of the outgoing top-of-atmosphere thermal flux are highly correlated
with that of cloud patchiness. Other than the vortices, turbulence and
transient waves with smaller horizontal length scales and higher oscillation
frequencies are also present. The basic findings from these global models
agree well with those of box models with a fixed Coriolis parameter $f$ across
the domain (Tan & Showman, 2020).

At a given rotation period, typical sizes of vortices are generally larger at
lower latitudes and smaller at higher latitudes. This is more prominent in
rapidly rotating models with rotation periods of 5 and 2.5 hours. The overall
sizes of vortices are larger with longer rotation period. This is because the
typical sizes of vortices driven by cloud radiative feedback are close to the
Rossby deformation radius $L_{d}=c_{g}/|f|$ (Tan & Showman, 2020), where
$c_{g}$ is the phase speed of gravity waves, $f=2\Omega\sin\phi$, $\Omega$ is
planetary rotation rate and $\phi$ is latitude. The dependence of $L_{d}$ on
$\Omega$ and $\phi$ naturally leads to the spatial and rotation-period
dependence of vortex sizes. Individual vortices are short lived and
experiences merging, disruption, shearing and straining over typical
timescales of tens of hours. There is no systematic migration of individual
vortices after they form. This is likely because the atmosphere is filled with
vortices, and individual vortex is significantly disrupted by adjacent
vortices before they can migrate over a long distance.

Figure 3: Time- and zonal-average cloud mass mixing ratio as a function of
latitude and pressure for models with four rotation periods of 20, 10, 5 and
2.5 hours (from the left to right) and with a drag timescale
$\tau_{\rm{drag}}=10^{5}$ s at the bottom. The white regions have
$\log_{10}q_{c}<-7$. Figure 4: Time- and zonal-average temperature-pressure
profiles at selected latitudes from the model with a rotation period of 2.5
hours and with a drag timescale $\tau_{\rm{drag}}=10^{5}$ s at the bottom.
Figure 5: Hovm$\ddot{\rm o}$ller diagrams (longitude-time sequence) showing
the time evolution of outgoing thermal flux at the equator for models with
four rotation periods of 20, 10, 5 and 2.5 hours (from the left to right) and
a drag timescale $\tau_{\rm{drag}}=10^{5}$ s. Fluxes shown in these panels are
sampled every 0.3 day.

As shown by Tan & Showman (2020), greater $|f|$ leads to overall thinner
clouds. Qualitatively, stronger rotation leads to a larger degree of
geostrophic balance in the flows which exerts greater suppression to the
vertical velocity, making the flows less efficient to vertically transport
tracers against cloud gravitational settling. Therefore, the variation of the
Coriolis parameter $f$ from the equator to the poles as well as the different
rotation period have significant consequences for the vertical extent of
clouds. In Figure 3 we show the time- and zonal-mean cloud mixing ratio as a
function of pressure and latitude for models with four rotation periods. With
a given rotation period, the vertical extent of clouds is largest at the
equator and decreases poleward. The cloud thicknesses at the equator are
almost the same among models with different rotation periods. At the poles,
the zonal-mean cloud thickness is smaller when the rotation period is shorter.
For rapidly rotating models, $f$ spans a larger range and is responsible for
the greater change of the zonal-mean cloud thickness. As a result of the
equator-to-pole cloud thickness variation, the time- and zonal-mean
temperature-pressure (T-P) structures exhibit a systematic equator-to-pole
variation due to the cloud radiative feedback. Figure 4 shows the time-
averaged zonal-mean T-P profiles at selected latitudes for the model with a
rotation period of 2.5 hours. The isobaric temperature variation can reach
more than 100 K from the equator to the poles. On top of this systematic
equator-to-pole zonal-mean variation, instantaneous differences in the T-P
structures associated with the cloudy and clear-sky regions are also present,
similar to that shown in Tan & Showman (2020).

####  3.1.2 Equatorial wave properties

Dynamics near the equator shows qualitative differences to that at mid-to-high
latitudes, which can be seen from two features. The first is the morphology of
the eddies. Eddies near the equator are likely more elongated in the zonal
direction, which is in contrast to the horizontally isotropic vortices at mid-
to-high latitudes. This feature is more prominent in models with rotation
periods of 5 and 2.5 hours. The second difference is the time evolution of the
eddies. Equatorial eddies exhibit systematic eastward or westward propagation,
while trajectories of mid-to-high-latitude vortices are more akin to random
walks and show no systematic migration along certain directions. Figure 5
shows the Hovm$\ddot{\rm o}$ller diagrams (longitude-time sequence) of
outgoing thermal flux at the equator as a function of longitude and time for
models with different rotation periods. All models exhibit obvious
characteristic eastward propagating patterns (those moving towards the upper
right of the panels) with zonal speeds on the order of a few hundred $\;{\rm
m\;s^{-1}}$. On the contrary, Figure 17 in Appendix A similarly shows the
Hovm$\ddot{\rm o}$ller diagrams at $45^{\circ}$ latitude where stochastically
evolving vortices dominate, and there is not evidence of eastward or westward
propagation of the eddies. These propagating eddies at low latitudes have much
larger zonal wavelengths than the mid-to-high-latitude vortices (comparing
Figures 5 and 17), showing dominant features characterized by zonal
wavenumbers of a few. Westward propagation of eddies is also present with a
slower phase speed. The time evolution of the zonally propagating eddies is
not always coherent, and often shows complications as they evolve. For
example, some existing perturbations disappear when they propagate while some
new perturbations are generated. The propagating eddies at low latitudes are
likely equatorially trapped waves, and they exert significant effects on the
short-term evolution of lightcurve variability as will be discussed in Section
3.3.

Based on both the morphology and time evolution of the eddies, latitudinal
boundaries that separate the two groups of atmospheric motions seem to emerge.
The boundary is closer to the equator with shorter rotation period. This
boundary is likely related to the equatorial deformation radius $L_{\rm
eq}=\sqrt{c_{g}/\beta}$. Poleward of this distance from the equator, there is
sufficient room for mature vortices to develop and the dynamics is dominated
by vortices. Within this scale from the equator, dynamics are shaped by
equatorially trapped waves that are coupled to the cloud radiative effect.
Using $c_{g}\sim 2000\;{\rm m\;s^{-1}}$ obtained from fitting theoretical off-
equatorial deformation radius to the measured sizes of dominant vortices as a
function of $f$ in the constant$-f$ models of Tan & Showman (2020), the
equatorial deformation radius $L_{\rm eq}$ extends to about $24^{\circ}$,
$17^{\circ}$, $12^{\circ}$ and $9^{\circ}$ away from the equator for models
with rotation periods of 20, 10, 5 and 2.5 hours, respectively. This is
roughly consistent with that shown in Figure 1 and 2 and their time evolution
(not shown).

Figure 6: Wavenumber-frequency relative power spectrum of the symmetric
components around the equator for models with four rotation periods of 20, 10,
5 and 2.5 hours and with a drag timescale $\tau_{\mathrm{drag}}=10^{5}$ s.
Solid lines are three sets of $n=1$ equatorially trapped Rossby waves with
three equivalent depths of $h_{e}=300$, 150 and 50 m (equivalent to Kelvin
wave phase speeds of 548, 387 and 224$\;{\rm m\;s^{-1}}$). Here $n$ is the
meridional mode number used in classic equatorial wave theory. Dashed lines
are three sets of Kelvin waves with the same equilibrium depths. The
horizontal axis represents zonal wavenumber with negative values representing
westward propagation and positive values representing eastward propagation.
The vertical axis is frequency in a unit of cycle per day (CPD). The relative
power spectrum shown here is non dimensional. Fluxes used in this analysis are
sampled every 0.3 day.  Figure 7: Wavenumber-frequency relative power spectrum
of the anti-symmetric components around the equator for the same data sets
used in Figure 6. Dotted lines are three sets of MRG waves with three
equivalent depths of $h_{e}=400$, 200 and 100 m (equivalent to having Kelvin
wave phase speeds of 632, 447 and 316$\;{\rm m\;s^{-1}}$). Fluxes used in this
analysis are sampled every 0.1 day.

Now we turn to characterize properties of equatorial waves that are present in
our simulations. We perform a spectral analysis at the equator in the
wavenumber-frequency domain, similar to that performed in Wheeler & Kiladis
(1999) and Showman et al. (2019). The brief procedure is the following. We
perform two-dimensional fast Fourier transforms on the outgoing thermal flux
as a function of longitude and time near the equator to obtain the raw Fourier
coefficients in the wavenumber-frequency space. These raw coefficients are
then heavily smoothed in the wavenumber-frequency space to generate the
background spectrum. Finally, the raw coefficients are divided by the
background spectrum to obtain the relative power spectrum, in which
significant signals will appear to have values greater than 1. The eddy fields
are further decomposed into symmetric and anti-symmetric components about the
equator because it helps to clarify the wave properties (Wheeler & Kiladis,
1999; Kiladis et al., 2009). The relative power spectrum of models with
different rotation periods are shown in Figure 6 for the symmetric components
and Figure 7 for the anti-symmetric components.

Our space-time spectral analysis demonstrates the robust existence of groups
of zonally propagating waves. In the symmetric components, a group of eastward
waves with zonal wavenumbers between 2 to 10 are present in all models. Weaker
signals of westward waves with zonal wavenumbers between $-10$ and $-2$ exist
in models with rotation periods of 5 and 2.5 hours but not in those with
rotation periods of 20 and 10 hours. In the anti-symmetric components, models
with rotation periods of 20 and 10 hours exhibit evidence of both eastward and
westward waves with zonal wavenumbers in between about $-5$ to $5$, and those
with 5 and 2.5 hours host westward waves with zonal wavenumbers between $-10$
to $1$.

We compare the signals in Figure 6 and 7 to dispersion relations of adiabatic,
freely propagating equatorial waves derived from the shallow-water system
(Matsuno, 1966). These waves are characterized by an equivalent depth $h_{e}$,
such that the gravity wave speed (which is the same as the Kelvin wave speed)
$c_{g}=\sqrt{gh_{e}}$ where $g$ is the surface gravity.333Linear equatorial
wave theories in the continuously stratified atmospheres often decompose the
atmosphere into a set of shallow-water systems with different equilibrium
depths corresponding to different vertical modes of the continuously
stratified atmosphere (see discussion in, e.g., Kiladis et al., 2009; Tsai et
al., 2014). This is partly why shallow-water equatorial waves are often used
as a baseline for comparisons to either observations or simulated atmospheres
(see a review by Kiladis et al., 2009).  Although waves in our simulations are
highly diabatic and tightly coupled to radiative effects of clouds, this
serves as a baseline comparison. Solid lines in Figure 6 represent the
dispersion relations of adiabatic free equatorial Rossby waves with a
meridional mode number $n=1$ and with three equivalent depths of 300, 150 and
50 m (equivalent to gravity wave phase speeds of 548, 387 and 224$\;{\rm
m\;s^{-1}}$). Dashed lines in Figure 6 represent dispersion relations of
Kelvin waves with the same equivalent depths. In the anti-symmetric components
shown in Figure 7, $n=0$ westward mixed Rossby-gravity (MRG) waves and $n=0$
and eastward inertia gravity (EIG) waves with three equivalent depths of 400,
200 and 100 m (gravity wave phase speed of 632, 447 and 316$\;{\rm
m\;s^{-1}}$) are plotted as dotted lines.

In the symmetric components (Figure 6), evidence of Kelvin waves are quite
strong for all models as the spectral powers follow the Kelvin wave dispersion
relations (dashed lines) reasonably well. Amplitudes of the Kelvin waves are
stronger among low zonal wavenumbers between $2-6$. The wave speeds of our
simulated Kelvin waves are somewhat dispersive. Kelvin modes with lower zonal
wavenumbers (longer zonal wavelengths) typically have slower phase speeds than
those with higher zonal wavenumbers. A phase speed of $\sim 550\;{\rm
m\;s^{-1}}$ brackets the upper limit of phase speeds in all models. The phase
speeds of waves with a zonal wavenumber of 2 in cases with rotation period of
10 and 5 hours are much slower than $220\;{\rm m\;s^{-1}}$. There is no
evidence of equatorial Rossby waves in models with rotation periods of 20 and
10 hours, and the evidence is tentative in models with 5 and 2.5 hours.
However, even in the case with 2.5 hours, the westward branch with relatively
high frequencies $>0.5$ cycle per day (CPD) obviously deviates from the
theoretical Rossby-wave dispersion relation. There are no wave signals at
frequencies higher than those shown in Figure 6 in the symmetric components,
suggesting the absence of high-frequency inertia gravity waves.

In the anti-symmetric components (Figure 7), MRG and $n=0$ EIG waves are
evident in cases with rotation periods of 20, 10 and 5 hours. They have either
eastward or westward phase velocities with speeds larger than that of the
Kelvin waves. In the case of 2.5 hours, although the wave signals lie in
between the dispersion relations of theoretical MRG waves assuming different
equilibrium depth, the signals show somewhat constant frequency between zonal
wavenumber -7 to -3, making the interpretation less obvious. Nevertheless,
these diagrams (Figure 6 and 7) quantify the propagation of waves and help to
recognize the wave types.

Figure 8:  Panel (a): a snapshot of the horizontal structure of eastwardly
propagating, symmetric wave at 0.23 bar with a zonal wavenumber 4 and a
frequency of about 0.15 CPD in the simulation with a rotation period of 10
hours. Thick solid lines are positive geopotential anomalies from $0.5\times
10^{4}$ to $3.5\times 10^{4}\;{\rm m^{2}s^{-2}}$ and thick dashed lines are
negative geopotential anomalies from $-3.5\times 10^{4}$ to $-0.5\times
10^{4}\;{\rm m^{2}s^{-2}}$ (with H marking the high-geopotential center and L
marking low-geopotential center). Arrows represent wind vectors, and colours
represent cloud mixing ratio associated with this wave. Panel (b): an example
of an adiabatic free Kelvin wave as a function of non-dimensional zonal and
meridional distances. Similar to panel (a), arrows represent wind vectors and
thick lines represent height perturbations in the shallow-water system. Panel
(c): An analytic forced-damped stationary wave solution forced by a zonal-
wavenumber-1 stationary forcing centered around (0,0). Light colours
represents larger values in the height field and arrows are eddy wind vectors.
This is adopted from Showman & Polvani (2011), which the details of the
solution is referred to.

Because eastward Kelvin waves seem dominant in both the equatorial time
evolution of out-going thermal flux (Figure 5) and the spectral analysis
(Figure 6), we further illustrate their properties in the physical domain by
extracting from the total dataset through filtering in the wavenumber-
frequency domain, similar to the procedure in Wheeler & Kiladis (1999). At a
given latitude, we choose a signal with a zonal wavenumber and a frequency in
the wavenumber-frequency domain, then project it back to the longitude-time
domain. We repeat this process for all latitudes near the equator and obtain
the propagating wave in the physical space. Panel (a) in Figure 8 shows the
snapshot of an eastward symmetric wave with a zonal wavenumber of 4 and a
frequency of about 0.15 CPD that shows a strong signal in the model with a
rotation period of 10 hours (the upper right panel in Figure 6). An example of
adiabatic free Kelvin wave is shown in panel (b) as a comparison to our
simulated wave. The adiabatic free Kelvin wave has eddy zonal velocity in
phase with the pressure anomaly, and the disturbances decay away from the
equator. Our simulated wave shows certain similarities to the adiabatic free
Kelvin wave — velocities have much larger zonal components than the meridional
components in most places, and the geopotential anomalies decay away from the
equator. Note that cloud abundances correlate well with the
convergence/divergence of the wind field (see a comparison between vectors and
the colors that represent the cloud abundance).

Meanwhile, in the simulated wave, the eastward zonal velocity pattern exhibits
a moderate eastward phase shift compared to the geopotential anomalies in the
zonal direction (whereas those in the adiabatic free Kelvin wave are in-
phase). Moreover, the simulated wave exhibits a northwest-southeast phase tilt
in the northern hemisphere and a southwest-northeast phase tilt in the
southern hemisphere for the geopotential anomalies. This likely reflects a
tendency that the circulation resembles wave patterns excited by a stationary
heat source that is symmetric about the equator (Matsuno, 1966; Gill, 1980).
An example of such a stationary wave solution subjected to damping on both
velocity and layer thickness in the shallow-water system adopted from Showman
& Polvani (2011) is shown in panel (c). This so-called Matsuno-Gill pattern is
characterized by an eastward shift of the Kelvin component at the equator and
a westward shift of Rossby components off the equator. Therefore it shows the
northwest-southeast phase tilt in the northern hemisphere and a southwest-
northeast phase tilt in the southern hemisphere, as well as that the maximum
positive zonal velocity is east of the height maximum at the equator. In our
simulations, when cloud forms by convergence at the equator, the accompanied
heat source drives the flow towards the Matsuno-Gill shape. However, because
clouds themselves are also advected by flows over a timescale that is
comparable to that required to form the Matsuno-Gill pattern, the advection of
the heat source disrupts the complete formation of the Matsuno-Gill pattern.
This is why our simulated wave still shows quantitative discrepancies to that
in panel (c).

It is likely that our simulated eastward, Kelvin-like waves resemble
properties of both the adiabatic free Kelvin waves and the forced-damped
Matsuno-Gill circulation. Recently, using a quasi-linear moist shallow-water
system, Vallis & Penn (2020) showed that the Matsuno-Gill-like circulation is
excited by latent heating (which itself is coupled to the flow, somewhat
analogous to the cloud radiative heating in our cases) at the equator, and a
slow eastward migration of the whole pattern is driven by the interactions of
the flow and moisture field. It is unclear how much the mechanism shown in
Vallis & Penn (2020) participates in the eastward propagation of our simulated
waves in addition to that related to the free Kelvin wave. Teasing out the
detailed mechanisms in our simulations is beyond the scope of this paper as
they are highly nonlinear. Here we still refer these zonally traveling
disturbances as waves, partly because these disturbances follow wave
dispersion relations reasonably well (as shown in Figures 6 and 7).

The origin of the equatorial waves in our simulations is likely related to
self-excitation by cloud radiative feedbacks. However, even treating them as
waves, mechanisms controlling the wave properties are unclear and little
previous work exists. Linear stability analyses on the cloud radiatively
coupled dynamics was carried out for inertia gravity waves in a $f$-plane
system and a quasi-geostrophic system by Gierasch et al. (1973), and they
showed that unstable modes are possible as a result of cloud radiative
instability. Here, we extend their theory to the equatorial $\beta$-plane
which is appropriate for equatorially trapped waves. The derivations are shown
in Appendix B. In the linear theory, we find that unstable modes are possible
for a set of Kelvin and $n=0$ MRG modes, suggesting a source of kinetic energy
on the equatorial waves shown in our models. However, the unstable modes from
the theory do not propagate. All other modes that resemble the classic
adiabatic free equatorial waves discovered by Matsuno (1966), including
propagating Kelvin, Rossby, MRG, eastward and westward inertia gravity modes,
show damping on the eddy amplitudes due to thermal radiation. The linear
theory is not a total failure because it still predicts the kinetic energy
sources for most equatorial waves in our simulations, but it cannot explain
the propagation of the simulated waves. This situation is similar to that in
the non-rotating two-dimensional system shown in Tan & Showman (2020).

The wave speeds shown in Figure 6 and 7 are significantly slower than expected
from adiabatic waves with long vertical wavelengths in conditions appropriate
for our simulated atmospheres, the latter of which is characterized by a much
larger dry equivalent depth $h_{e,{\rm dry}}=c^{2}_{g,{\rm dry}}/g\sim 4000$
m, where $c_{g,{\rm dry}}\approx 2NH\sim 2000\;{\rm m\;s^{-1}}$, $N$ is the
Brunt-Vaisala frequency and $H$ is the scale height. This is similar to
equatorial waves in Earth’s troposphere which are affected by latent heat
released from moist convection (Wheeler & Kiladis, 1999). To understand the
reduced phase speed of tropical waves by diabatic effects, the following
idealized framework has been proposed (see a review by Kiladis et al., 2009).
Suppose that the large-scale diabatic heating and cooling $Q$ is included in
the linearized thermodynamics system as the following

$\frac{\partial}{\partial t}\left(\frac{\partial\phi}{\partial
z}\right)+wN^{2}=Q,$ (1)

where $\phi$ is the geopotential (note that $\partial\phi/\partial z$ is
proportional to the temperature perturbation), $z=-H\log(p/p_{s})$ is the log-
pressure coordinate, $p$ is pressure, $p_{s}$ is a reference pressure, and $w$
is vertical velocity in log-pressure coordinates. Furthermore, if the heating
and cooling are proportional to the vertical velocity such that $Q=\alpha
wN^{2}$ (where $\alpha$ is an arbitrary constant), then Equation (1) becomes

$\frac{\partial}{\partial t}\left(\frac{\partial\phi}{\partial
z}\right)+w(1-\alpha)N^{2}=0.$ (2)

If $\alpha$ is positive and less than 1, diabatic heating and cooling
effectively reduce the stability of the atmosphere and therefore reduce the
wave phase speeds. For Earth’s tropical waves, the problem reduces to
determine theoretical values based on parameterized moist convection for
$\alpha$ to explain the observed wave speeds (e.g., Neelin & Held, 1987;
Emanuel et al., 1994; Haertel & Kiladis, 2004), but the problem has not been
completely solved (Kiladis et al., 2009). For waves coupled purely with cloud
radiative feedback, the relation $Q=\alpha wN^{2}$ likely holds as well. In
regions with upwelling, vapor is advected above the condensation level and
clouds form, which generates warming near the condensation level due to cloud
radiative effect. In regions with downwelling, cloud-free air is advected
downward and clears out the region, enhancing the IR flux to space and
inducing cooling. To quantitatively confirm the positive correlation between
upwelling and heating, or downwelling and cooling, we calculate the cospectral
power density for the quantity $wQ$ in the cloud forming layer, which is
simply $2\mathbb{R}(w_{k}Q^{\ast}_{k})$ where $w_{k}$ and $Q_{k}$ are the
coefficients at wavenumber $k$ space for vertical velocity and heating rate,
and $Q^{\ast}_{k}$ is the conjugate of $Q_{k}$. A detailed description of
similar exercises in constant$-f$ models is referred to Tan & Showman (2020).
The $wQ$ cospectral power density shows positive values for almost all zonal
wavenumbers, proving that $Q=\alpha wN^{2}$ with a positive $\alpha$ holds
well in our models. Therefore, we expect the above idealized framework may be
used to qualitatively understand the reduced wave speeds seen in our
simulations (however it does not address questions as why and how certain wave
modes are excited but others are not). As for the detailed determination of
$\alpha$, we defer to a future study.

###  3.2 Results with weaker bottom drag

Figure 9: Time-averaged and zonal-mean zonal velocity as a function of
pressure and latitude for models with a rotation period of 5 hours and
different drag timescales $\tau_{\rm{drag}}=10^{5}$, $10^{6}$ and $10^{7}$ s
from the top to bottom row.

We examine properties of the circulation when the bottom frictional drag
becomes weaker by performing two additional models with drag timescales of
$\tau_{\mathrm{drag}}=10^{6}$ s and $10^{7}$ s and a rotation period of 5
hours. Note that in all models, the frictional drag is applied only at
pressures greater than 5 bars, well below the cloud condensation level of 0.5
bar. Figure 9 shows the time-averaged zonal-mean zonal wind as a function of
pressure and latitude for models with drag timescales of
$\tau_{\rm{drag}}=10^{5}$, $10^{6}$ and $10^{7}$ s from the top to the bottom
panel. When the drag is strong ($\tau_{\mathrm{drag}}=10^{5}$ s), the zonal-
mean zonal wind is much weaker than local horizontal wind speeds of vortices
and turbulence. The zonal-mean zonal flow at low latitudes exhibits an
interesting vertical wind shear, with a westward mean flow centering at around
0.5 bar and an eastward mean flow around 0.01 bar. The equatorial eastward jet
centering around 0.01 bar corresponds to a superrotation, which requires up-
gradient angular momentum transport to the equator by eddies. Our analysis
(not shown) suggests that horizontal eddy momentum transport by transient
waves is responsible for this local superrotation. This might be somewhat
analog to those proposed for superrotation in solar-system bodies, such as
Venus, Titan, tropospheres of Jupiter and Saturn (see a recent review by
Imamura et al., 2020). In the model with $\tau_{\mathrm{drag}}=10^{6}$ s, a
broad equatorial westward jet and high-latitude eastward jets with speeds
$\sim 400\;{\rm m\;s^{-1}}$ emerge. At low latitudes, the vertical wind shear
is similar to that with $\tau_{\mathrm{drag}}=10^{5}$ s despite the overall
equatorial jet velocity being westward. The jet speeds are comparable to the
horizontal mean eddy velocity, but the vortices are still nearly isotropic at
mid-to-high latitudes. Near the equator, the equatorial waves are Doppler
shifted by the westward jet whose speed is comparable to that of wave
propagation. In the model with $\tau_{\mathrm{drag}}=10^{7}$ s, the jet speeds
become much stronger, reaching about $-2200\;{\rm m\;s^{-1}}$ at the equator
and more than $1500\;{\rm m\;s^{-1}}$ at about $\pm 50^{\circ}$ latitude.
Compared to the jet structure with $\tau_{\mathrm{drag}}=10^{6}$ s, the
meridional width of the equatorial jet with $\tau_{\mathrm{drag}}=10^{7}$ s is
very similar and only the jet speed increases. However, at high latitudes, the
cores of eastward jets are closer to the equator than those with
$\tau_{\mathrm{drag}}=10^{6}$ s. There is a strong barotropic (pressure-
independent) component for zonal jets with $\tau_{\mathrm{drag}}=10^{7}$ s.

The RMS horizontal eddy velocity in the cloud forming regions is about
$400-550\;{\rm m\;s^{-1}}$, which is much larger than the jet speeds in the
model with $\tau_{\mathrm{drag}}=10^{5}$ s. Even in the model with
$\tau_{\mathrm{drag}}=10^{6}$ s, the jet speeds are only comparable to the
eddy velocities. In the following analysis, we focus primarily on the model
with $\tau_{\mathrm{drag}}=10^{7}$ s in which the jet speeds well exceed the
eddy velocities and the dynamics exhibits obvious differences to those with
$\tau_{\mathrm{drag}}=10^{5}$ s.

Figure 10: Instantaneous results from the model with a rotation period of 5
hours and a drag timescale $\tau_{\rm{drag}}=10^{7}$ s after the model reached
a statistical equilibrium state. Panel (a) shows the cloud mixing ratio at
0.23 bar; panel (b) shows temperature at 0.23 bar; and panel (c) shows the
absolute vorticity at 0.23 bar. Note that strong eastward jets are located at
$\pm 48^{\circ}$ latitude.

Cloud formation is affected by the presence of strong zonal jets in the case
with $\tau_{\mathrm{drag}}=10^{7}$ s. Wind shears at the flanks of the mid-
latitude jets create strong cyclonic regions poleward of the jets and strong
anti-cyclonic regions equatorward of the jets. Anti-cyclones which are
associated with cloud formation tend to be more vulnerable in the strong
cyclonic zones, in the sense that they are sheared apart and destroyed more
easily than in other regions. On the other hand, formation of cyclones which
are associated with thin clouds are more prevalent in the cyclonic zones. In
addition, vortices tend to migrate towards regions with background absolute
vorticity closer to that of vortices (e.g., Scott, 2011; O’Neill et al.,
2015). There is a tendency that cyclones formed near the mid-latitude jets
migrate to the poleward flanks of the jets where it is cyclonic, whereas anti-
cyclones tend to migrate to the equatorward flanks of the jets. These
behaviors are observed in time evolution of storms in the simulation with
$\tau_{\mathrm{drag}}=10^{7}$ s (not shown). All these imply that, clouds
associated with anti-cyclones are relatively depleted in the polar flanks of
the mid-latitude jets and are more abundant in the equatorial flanks of the
jets. The vortex behaviors influenced by the strong zonal jets described above
may be one of the mechanisms responsible for the cloud mixing ratio at 0.23
bar in panel (a) of Figure 10. There are also zonal-mean upwelling equatorward
of the mid-latitude jets where it is cloudy and zonal-mean downwelling
poleward of the mid-latitude jets where it is less cloudy (not shown). This
mean meridional circulation may or may not be a mechanism driving the
meridional cloud gradient. It could simply be radiatively driven given the
existing meridional cloud gradient. Shapes of vortices are deformed due to the
significant horizontal wind shear, especially near latitudes between
$15^{\circ}-45^{\circ}$ and $50^{\circ}-70^{\circ}$ (see Figure 10).
Interestingly, there is a strong polar cyclone at each pole with the center of
the cyclone not aligned with the pole. Clouds are also strongly depleted in
the cyclones, and the same reasons mentioned above for cloud depletion
poleward of the mid-latitude jets might also be relevant. An instantaneous
temperature map at 0.23 bar is shown in panel (b) and instantaneous absolute
vorticity $f+\zeta$ is shown in panel (c) of Figure 10, where
$\zeta=\mathbf{k}\cdot\nabla_{p}\times\mathbf{v}$ is the relative vorticity,
$\mathbf{v}$ is horizontal velocity vector, $\mathbf{k}$ is the local upward
unit vector on the sphere and $\nabla_{p}$ is the horizontal gradient in
pressure coordinates. Due to the radiative effects of clouds, strong cyclonic
regions where clouds are relatively depleted are systematically colder than
anti-cyclonic regions.

The development of robust zonal jets with a weaker bottom frictional drag is
consistent with that found in Showman et al. (2019), who studied atmospheric
circulation of brown dwarfs using horizontally isotropic, randomly evolving
thermal forcing. Storms and vortices excite Rossby waves, and their
generation, propagation and wave breaking interact with the mean flow, driving
zonal jets (Dritschel & McIntyre, 2008). Without efficient removal of kinetic
energy by the strong bottom frictional drag, zonal jets can be pumped up and
maintained by wave-mean-flow interactions. Some features of the jets in our
simulations are interesting. First, the jets are quite meridionally broad,
with the equatorial westward jet extending to $\pm 40^{\circ}$ latitude and
subsequent eastward jets at about $\pm 70^{\circ}$ latitude for the case with
$\tau_{\mathrm{drag}}=10^{6}$ s and at about $\pm 48^{\circ}$ latitude for the
case with $\tau_{\mathrm{drag}}=10^{7}$ s. As a comparison, Jupiter has about
7 subtropical zonal jets in each hemisphere (e.g., Ingersoll et al., 2004).
Classic turbulence-driven jet theory predicts that the meridional jet spacing
is related to wind speed via the Rhines scale $\pi\sqrt{U/\beta}$ where $U$ is
a characteristic wind speed, indicating that the number of jet on a sphere is
roughly given by $N_{\rm jet}\sim\sqrt{2\Omega a/U}$ where $\Omega$ is the
rotation rate and $a$ is planetary radius (see reviews by, e.g., Vasavada &
Showman, 2005; Dritschel & McIntyre, 2008; Showman et al., 2010). Given a
typical wind speed of several hundreds of $\;{\rm m\;s^{-1}}$ in the case with
$\tau_{\rm{drag}}=10^{6}$ s and $\sim 2000\;{\rm m\;s^{-1}}$ in the case with
$\tau_{\rm{drag}}=10^{7}$ s, one would expect a number of zonal jets of $\sim
9$ in the case with $\tau_{\rm{drag}}=10^{6}$ s and $\sim 5$ in the case with
$\tau_{\rm{drag}}=10^{7}$ s, which are obviously greater than that seen in our
simulations. Second, the equatorial jets are westward in our simulations,
whereas those in Jupiter and Saturn are eastward.

Figure 11: Left panel: time- and vertically averaged zonal-mean zonal wind as
a function of latitude for the model with a rotation period of 5 hours and a
drag timescale $\tau_{\rm{drag}}=10^{7}$ s after the model reached a
statistical equilibrium state. Right panel: zonal-mean absolute vorticity
($f+\overline{\zeta}$) corresponding to the zonal-mean zonal wind in the left
panel and planetary vorticity $f$.

The broader-than-expected jets are probably related to the efficient
horizontal mixing of potential vorticity (PV) caused by strong vortices, which
tends to break the PV staircase and smooth out the jet structure. One
interpretation of the Rhines jet scaling is that the jets are a natural result
of PV staircases. Eastward jets correspond to the PV discontinuity whereas the
westward jets correspond to the PV homogenization, and jet speed is determined
by the meridional width of the PV staircase (e.g., Dritschel & McIntyre, 2008;
Scott & Dritschel, 2012). But if the magnitude of vorticity sources is much
larger than the background PV discontinuity, mixing of PV between the
staircases leads to destruction of the PV staircase, therefore the jet scaling
does not apply well in this situation (Scott & Dritschel, 2012). In strongly
forced and damped cases of Showman et al. (2019), the PV structure is
disrupted by eddy vorticity and the resulting zonal jets are less sharp and
meridionally broader than those that are weakly forced and damped. In our
case, the vorticity mixing is even more extreme. We find that in our
simulations, the vorticity associated with mid-to-high-latitude vortices are
typically much stronger than either the background planetary vorticity or that
associated with the zonal jets. These vortices occurs somewhat randomly in
space and time—their occurrence is not particularly constrained by the
presence of zonal jets—and they can easily disrupt the zonal-mean PV
structure. In panel (c) of Figure 10, we can see that the eddy relative
vorticity is far greater than the background planetary vorticity in the case
with $\tau_{\mathrm{drag}}=10^{7}$ s. Figure 11 shows the vertically averaged
zonal-mean zonal wind in the left panel and its corresponding absolute
vorticity (together with the planetary vorticity) in the right panel. Because
the barotropic component dominates the jet structure in the case with
$\tau_{\mathrm{drag}}=10^{7}$ s (panel (c) in Figure 9), the absolute
vorticity is a good representation of PV. The eastward jets correspond to the
sharp gradient of absolute vorticity at $\pm 48^{\circ}$ latitude. There is a
lack of other absolute vorticity staircase in other regions. Meanwhile,
absolute vorticity between $\pm 40^{\circ}$ latitude tends to be homogenized
(although still far from being completely well mixed). In some cases, the
cross-equator homogenization of absolute vorticity can lead to a strong
equatorial westward jet (Dunkerton & Scott, 2008). Such a tendency may be
responsible for the broad, strong equatorial westward jets in our simulations.

Figure 12:  Wavenumber-frequency relative power spectrum of the symmetric
component around the equator (on the left panel) and of the anti-symmetric
component (on the right panel) for the model with a rotation period of 5 hours
and $\tau_{\mathrm{drag}}=10^{7}$ s. Doppler shifts in the power spectrum
associated with a uniform zonal-mean zonal wind $-2130\;\;{\rm m\;s^{-1}}$ are
removed, and most signals shown in these panels are likely associated with
waves propagating relative to the zonal-mean zonal wind. On the left panel,
solid lines are theoretical equatorial Rossby waves and dashed lines are
theoretical Kelvin waves, same as those shown in Figure 6 for a rotation
period of 5 hours. On the right panel, dotted lines are theoretical MRG waves,
same as those shown in Figure 7 for a rotation period of 5 hours.

There is a local absolute vorticity maximum northward of the jet at
$48^{\circ}$ latitude and a local minimum southward of the jet, indicating
that the jet violates the barotropic stability criteria. This is associated
with the accumulation of cyclones northward of the jet and anti-cyclones
southward of the jet. A positive feedback likely takes place in this case: the
jet structure triggers accumulation of vortices, and the accumulation of those
vortices further strengthens the jet structure that already promotes the
accumulation. Similar influences of the migration and accumulation of vortices
on the jet structure has been proposed in previous modeling studies of Jovian
atmospheric dynamics (Thomson & McIntyre, 2016). A balance may be reached
between the accumulation of vortices and the instability associated with the
jet that tends to restore the jet structure towards stability.

In box simulations with constant Coriolis parameter $f$ performed in Tan &
Showman (2020), we showed that when the bottom frictional drag is weak
($\tau_{\rm drag}=10^{7}$ s), kinetic energy accumulates and a pair of a
cyclone and a anticyclone forms with sizes comparable to the simulated domain
size. In the global domain, the presence of $\beta$-effect excites Rossby
waves, and eddy-mean-flow interactions channel the kinetic energy to the zonal
direction, forming zonal jets instead of ever larger vortices. Our simulations
demonstrate the importance of the $\beta$ effect on the formation of zonal
jets.

Finally, equatorially trapped waves also exist in the weak-drag models but
their propagation is influenced by the equatorial jets. The influence is most
prominent in the case with $\tau_{\mathrm{drag}}=10^{7}$ s, in which
equatorial waves are Doppler shifted by the equatorial jet with a zonal-mean
zonal wind $\sim-2200\;\;{\rm m\;s^{-1}}$. In the wavenumber-frequency
analysis for this case, we remove the Doppler shifts associated with the jet
for which we assume a uniform zonal-mean zonal wind of $-2130\;\;{\rm
m\;s^{-1}}$ according to Figure 11. The results are shown in Figure 12. Part
of the signals shown in these panels are likely associated with waves
propagating relative to the zonal-mean zonal wind. We find evidence of
eastward Kelvin waves and westward Rossby waves in the symmetric component
around the equator. A family of MRG waves is also indicated in the analysis of
the anti-symmetric component although their signals slightly deviate from the
theoretical ones. However, in reality, the equatorial jet has horizontal
shears, whose influences on waves cannot be removed by simply assuming a
uniform zonal-mean zonal wind in the analysis shown in Figure 12. Such
influences may be responsible for the abnormally strong westward Rossby wave
signal (which is only tentative in the strong-drag cases shown in Figure 6) in
the left panel, as well as for signals around zonal wavenumber 3-5 and
frequency 0.1-0.2 CPD in the right panel, which is absent in the strong drag
cases (Figure 7). Indeed, when we slightly change the assumed uniform zonal-
mean zonal wind in the wavenumber-frequency analysis, the position or strength
of the above “abnormal" signals is sensitive to the assumed zonal-mean zonal
wind. But the signals associated with the Kelvin waves and MRG waves are only
moderately affected, indicating that Kelvin waves and MRG waves seem robust.

###  3.3 Simulated lightcurves

Figure 13:  The top four panels (a to d) show normalized lightcurves as a
function of rotation period for models with four different rotation periods of
20, 10, 5 and 2.5 hours (from the top to bottom panels) and with a drag
timescale $\tau_{\rm{drag}}=10^{5}$ s. In each panel, the black line
represents a viewing angle of $0^{\circ}$ (equator-on), and the red line
represents a viewing angle of $45^{\circ}$, and the blue line represents a
viewing angle of $90^{\circ}$ (northern-pole-on). The bottom two panels show
the normalized lightcurves for models with a rotation period of 5 hours and
$\tau_{\mathrm{drag}}=10^{6}$ s in panel e and $\tau_{\mathrm{drag}}=10^{7}$ s
in panel f (in which a strong westward equatorial jet develops).

Due to the significant inhomogeneity in the outgoing thermal flux, lightcurve
variability is expected from our simulations. Figure 13 shows the simulated
lightcurves normalized to their time-mean values as a function of time that is
normalized by the rotation periods of the models. Panels a to d are
lightcurves from models with $\tau_{\mathrm{drag}}=10^{5}$ s and rotation
periods of 20, 10, 5 and 2.5 hours, respectively. Panels e and f are
lightcurves from models with a rotation period of 5 hours and
$\tau_{\mathrm{drag}}=10^{6}$ and $10^{7}$ s, respectively. Black, red and
blue lines in each panel represent viewing angles of $0^{\circ}$, $45^{\circ}$
and $90^{\circ}$ relative to the equator, respectively—$0^{\circ}$ means an
equator-on view and $90^{\circ}$ means a northern-pole-on view. The former
maximizes the rotational modulation of the lightcurve whereas in the latter
there are no rotational effects. The atmospheres in our models are
statistically symmetric between the northern and southern hemisphere, and
therefore simulated lightcurves viewed from two hemispheres are qualitatively
similar.

The simulated lightcurves exhibit several important characteristics. For all
models, the amplitudes of lightcurve variability are maximized when viewed
equator-on, whereas they are minimized when viewed pole-on. Part of the reason
is because the horizontal length scales of storms are larger at low latitudes,
which causes larger flux perturbations when the object is viewed equator-on.
This is consistent with the viewing angle dependence of near-IR flux
variability amplitude found by Vos et al. (2017). The peak-to-peak amplitude
of the equator-on lightcurve is typically a few percent and decreases with
decreasing rotation period—the amplitude is almost 4% for the case with
rotation period of 20 and 10 hours, and slightly more than 2% for the case
with 5-hour rotation, and finally slightly more than 1% for the case with
2.5-hour rotation. These characteristics are further summarized in Figure 14,
which shows the peak-to-peak amplitudes of normalized lightcurves with
different viewing angles.

Figure 14: Peak-to-peak variability amplitudes for normalized lightcurves with
different viewing angles from models with different rotation periods and
$\tau_{\mathrm{drag}}=10^{5}$ s. A subset of lightcurves are shown in Figure
13.

Lightcurves with rotation periods of 10, 5 and 2.5 hours exhibit clear
periodicity related to the rotation period at least during certain times in
their evolution. However, at some points the rotational periodicity can be
complicated by the time evolution of the surface patchiness. For instance, in
the case with a rotation period of 5 hours, the 5-hour periodicity is clear
before the 12th rotation period, but transition to double peaks within a
rotation period after that. By eye, the shape of lightcurves for a rotation
period of 20 hours is more complicated by short-term irregularities. This is
perhaps because the rotation period is long and the typical sizes of vortices
are relatively large, which together indicate that evolution of individual
vortices (typically over a timescale of several tens of hours) may
significantly impact the lightcurve evolution over rotational timescales.

Figure 15: Normalized lightcurves with an equator-on viewing angle for models
with four different rotation periods of 20, 10, 5 and 2.5 hours (from the top
to bottom panels) and with a drag timescale $\tau_{\rm{drag}}=10^{5}$ s. In
each panel, the thick grey line is the full-model lightcurve. The red line is
from thermal flux outputs wherein the equatorial region is held fixed at time
zero and stays invariant, while the mid-to-high-latitude regions are not
modified. On the contrary, the blue line is from outputs wherein the mid-to-
high-latitude regions are held fixed at time zero. The equatorial regions here
are defined as between $\pm 25^{\circ}$ for the model with a rotation period
of 20 hours, and between $\pm 20^{\circ}$ for 10 hours, and between $\pm
15^{\circ}$ for both 5 and 2.5 hours.

Viewed from pole-on, lightcurves are much smoother and show much longer
periods of variation. The amplitude of pole-on lightcurves decreases with
decreasing rotation period. This lightcurve variability is caused by the
statistical fluctuations of the ensemble of storms as discussed in Tan &
Showman (2020). The fewer storms projected in the disk, the larger the effects
that the evolution of an individual storm can have on the the lightcurve.
Sizes of storms near the poles inversely decrease with decreasing rotation
period, therefore the above trend found in pole-on lightcurves emerges. Even
though smaller than those viewed equator-on, the peak-to-peak amplitudes of
pole-on lightcurves can still reach almost 2% for cases with 20 and 10 hours
and 1% for the case with 5 hours, which are detectable given sensitivities of
current instruments (e.g., Wilson et al., 2014; Metchev et al., 2015). This
may contribute to some long-term variations in observed lightcurves that are
not easily explained by rotation modulation (e.g., some samples can be seen in
Metchev et al., 2015).

Weak-drag models with $\tau_{\mathrm{drag}}=10^{6}$ and $10^{7}$ s similarly
exhibit significantly time-varying waves and vortices, and some properties of
their lightcurves are similar to those with $\tau_{\mathrm{drag}}=10^{5}$ s,
including peak-to-peak normalized lightcurve amplitude of a few percent,
certain irregular time variability, viewing-angle dependence of variability
amplitude, and complications on the periodicity due to evolution of cloud
patchiness. However, because weak-drag models develop meridionally broad,
strong zonal jets, zonal advection of clouds by the jets can modify the
periodicity of the lightcurves. As shown in Figure 9, the equatorial jets are
both westward with speeds of a few hundred $\;{\rm m\;s^{-1}}$ and about
$2000\;{\rm m\;s^{-1}}$ in models with $\tau_{\mathrm{drag}}=10^{6}$ and
$10^{7}$ s, respectively. Although Kelvin waves have eastward phase speeds
relative to the mean flow, in the rotational frame they travel to the west due
to the strong westward equatorial jet. Their lightcurves, especially viewed
equator-on, are expected to show periodicity longer than the rotational period
of the model. This is most prominent in the case with
$\tau_{\mathrm{drag}}=10^{7}$ s in panel f of Figure 13, wherein the
lightcurve shows only 18 periods over 20 rotational periods.

We divide contributions to the lightcurve variability by surface
inhomogeneities into two groups of dynamical features: the zonally propagating
equatorial waves and mid-to-high-latitude vortices (the latter do not migrate
along the zonal direction). We diagnose the effects of the two groups in the
lightcurve by the following process. First, in the model outputs, we
artificially hold the equatorial region to be static. A synthetic lightcurve
is generated based on this configuration, and the time evolution of the shape
of the variability is caused only by evolution of mid-to-high-latitude
vortices. Then, we artificially hold the mid-to-high-latitude regions to be
static but do not hold the equatorial region static. The resulting time
evolution of the shape of the variability represents only effects of the
propagating equatorial waves. Figure 15 shows these experiments for four
models with $\tau_{\mathrm{drag}}=10^{5}$ s, in which the thick grey lines are
the original full lightcurves; the red lines are lightcurves wherein the
equatorial regions are held fixed from time zero; and the blue lines are
lightcurves wherein mid-to-high latitudes are held fixed from time zero. All
cases are viewed equator-on. The equatorial region is defined as in between
$\pm 25^{\circ}$ latitude for the case with rotation period of 20 hours,
between $\pm 20^{\circ}$ latitude with 10 hours, and between $\pm 15^{\circ}$
latitude with 5 and 2.5 hours. These latitudinal bands are chosen to safely
include the equatorial trapped waves. Slightly changes of these latitudes do
not influence our results and conclusions.

In general, both the equatorial waves and mid-to-high-latitude vortices
contribute to the short-term evolution of the full lightcurve
variability—removing the time evolution of either component results in
significant changes in the lightcurves. However, equatorial waves impact the
time evolution of lightcurve shapes in more critical ways. First, the wave
beating effect that causes the change of variability amplitude with time is
much weaker when the equatorial regions are held fixed. This is most obvious
in the case with a rotation period of 5 hours, in which the red curve (for
which equatorial region is held time invariant) has almost a constant
amplitude, and the splitting to double sub-peaks shown in the original
lightcurve does not occur. Second, there are significant phase differences in
the lightcurve variability between the full lightcurve and the lightcurve with
equatorial regions fixed. This is obvious in all cases and we take the case
with 2.5-hour rotation as an example: the red curve starts to show phase
differences relative to the thick grey curve at a time of about the 7th
rotation period; towards the end the red and grey curves show an almost
$180^{\circ}$ phase difference. The blue curve, for which mid-to-high-latitude
regions are held time invariant, mostly only show differences in the local
peaks and troughs of the variability, but not in the long-term evolution of
the amplitude and the phase of the variability. Our diagnostic analysis
suggests that equatorial waves have major impacts on the lightcurve
variability and the time evolution of the shape of the variability, and mid-
to-high-latitude vortices contribute to the local features of the lightcurves.
This is in good agreement with the analysis of long-term observed lightcurves
of a few BDs by Apai et al. (2017), in which they found that waves can explain
the major evolution of lightcurves, and “spots” are needed to fit the
remaining local inconsistency between data and the wave model.

Figure 16: Periodograms of normalized lightcurves for models with four
rotation periods of 20, 10, 5 and 2.5 hours and with a drag timescale
$\tau_{\mathrm{drag}}=10^{5}$ s (panel a to d). Periodogram for the model with
a rotation period of 5 hours and $\tau_{\mathrm{drag}}=10^{6}$ s is shown in
panel e, and for that with $\tau_{\mathrm{drag}}=10^{7}$ s is in panel f. In
panel c, additional dotted lines are from output samples at different times
much later than that for the solid blue line.

Both eastward and westward waves are present at low latitudes, and the
eastward Kelvin waves are sometimes more dominant as visually shown in Figure
5 and in the wavenumber-frequency analysis (Figure 6 and 7). This implies a
faster rotation of the equatorial features than the intrinsic planetary
rotation, which may lead to a slightly shorter period of the lightcurve
variability than the underlying rotation period. This effect has been clearly
shown in the case with 2.5-hour rotation in Figure 15, and now we quantify
this using periodogram analyses of the simulated lightcurves shown in Figure
13. Figure 16 shows power densities (in arbitrary unit) by the periodogram
analysis for normalized lightcurves viewed equator-on for our models.

Several features are interesting in the periodograms for models with
$\tau_{\mathrm{drag}}=10^{5}$ s. First, some cases exhibit slight shifts of
the major peaks to a period slightly shorter than the rotation period. In
particular, cases with 20 and 2.5 hours show obvious shifts. This is because
of the eastwardly propagating Kelvin waves with phase speeds of a few hundred
$\;{\rm m\;s^{-1}}$. Second, the power densities exhibit major peaks very
close to (in some cases exactly the same as) the underlying rotation period of
the models. However, these peaks are broadened, indicating that the dynamical
features, either eastward or westward with various phase speeds, could
contaminate the periodicity of the lightcurve variability. Third, there is
usually a secondary peak of the power density at a period of approximately
half of the rotation period. This may be caused by the equatorial features
with a zonal wavenumber 2 that appears twice to the “observer" as the objects
rotate once. Finally, our models have shown that different equatorial waves
that travel in different zonal velocities may contribute to the flux
patchiness differently at different times (Figure 5). This could cause time
variability on the properties of lightcurves when the lightcurves are sampled
over a finite period (as the real observations do). This is illustrated for
the case with 5-hour rotation and $\tau_{\mathrm{drag}}=10^{5}$ s in panel c
of Figure 16. We show two additional periodograms of lightcurves that sample
the model outputs at different times, and each of them also samples about 20
rotational periods as the original one. Their relative shapes differ slightly,
and importantly, their major power density peaks can be at periods both longer
and shorter than the underlying rotation period. This suggests that not only
these waves can cause differences between the “observed" lightcurve
periodicity and the underlying rotation period, but also the degree of this
deviation may vary with time.

For models with weaker frictional drags, the development of strong equatorial
westward jets induce shifts in the period of variability. Panel e and f in
Figure 16 show periodograms for the case with $\tau_{\mathrm{drag}}=10^{6}$
and $10^{7}$ s, respectively. The first-order feature of them is that their
major peaks shift to longer periods compared to the rotation period, and the
shift is stronger in the case with $\tau_{\mathrm{drag}}=10^{7}$ s.

##  4 Discussion and Conclusions

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

Recently, Vos et al. (2017), Vos et al. (2018) and Vos et al. (2020) suggested
a viewing angle dependence of the observed near-IR colors and variability
amplitude for a large sample of field BDs. Our dynamical models provide
support for their observational results. Larger variability amplitude when
viewed more equator-on is a natural result of the equatorial maximum of
rotational variability along with our finding that cloud patches reach maximum
sizes at low latitudes (Figure 1, 2, 13 and 14). The vertical extent of zonal-
mean cloud mixing ratio is higher at lower latitudes (Figure 3), which could
be responsible for the redder near-IR colors when viewed more equator-on. The
near-IR colors show a wide scatter in the color-magnitude diagram for mid-to-
late L dwarfs (see color-magnitude diagrams for a large sample of BDs in,
e.g., Faherty et al., 2016; Liu et al., 2016). Due to the latitudinal
variation of the cloud thickness at a given rotation period and the dependence
of global-average cloud thickness on varying rotation period, the different
viewing angles and the variation of rotation periods of the field BDs might
contribute to the scatter of observed near-IR colors. Millar-Blanchaer et al.
(2020) showed that assuming two broad zonal bands with different cloud
properties in each hemisphere, they can reproduce time-averaged polarization
measured for the nearby BD Luhman 16B. Our models do not show clear zonally
banded cloud structure like those in Jupiter and Saturn, but exhibit smooth
equator-to-pole cloud thickness variations. It is worthwhile to explore how
can polarization be produced by the smooth equator-to-pole cloud thickness
variation with different viewing angles and compare it to the measured
polarization.

Amplitudes of our simulated lightcurves typically range from 0.5 percent to a
few percent depending on the rotation period and viewing angle, consistent
with those found in the majority of observed lightcurves (Buenzli et al.,
2014; Radigan et al., 2014; Radigan, 2014; Wilson et al., 2014; Metchev et
al., 2015; Vos et al., 2018). We did not reproduce variability with amplitude
much greater than a few percent. Yet, several field BDs and free-floating low-
mass objects have shown large variability with peak-to-peak amplitudes greater
than 10 percent (e.g., Buenzli et al., 2012; Apai et al., 2013; Buenzli et
al., 2015; Lew et al., 2016; Apai et al., 2017; Zhou et al., 2020). Either
cloud patches with sizes larger than those in our simulations or a greater
horizontal contrast of outgoing thermal flux are required to explain those
very large varability amplitudes. Our models only occupy a very limited
parameter space. Further exploring parameter space, including broad range of
surface gravity, atmospheric temperature and varying cloud particle size, will
yield richer dynamical behaviors.

Our models are highly idealized in the sense that radiative transfer and cloud
formation (including the cloud microphysics and the treatment of sub-grid-
scale structure) are vastly simplified in order to emphasize the role of
atmospheric dynamics. Our modeled cloud structures capture the first-order
behavior of cloud formation that has been shown by previous cloud formation
models of BDs—a sharp cloud base typically emerges around the condensation
level and the cloud mixing ratio smoothly decreases with decreasing pressure
due to mixing (see reviews by, e.g., Helling & Fomins, 2013; Marley &
Robinson, 2015; Helling, 2019). Our cloud scheme neglects temperature feedback
on cloud formation, which could influence locations of the cloud base, and it
certainly does not capture all the sophisticated microphysical processes.
Future endeavours should include better representation of radiative transfer,
cloud microphysics and parameterization of sub-grid cloud structure.

Finally, we summarize our major findings in this study as follows:

  * •

Vigorous atmospheric circulation can be triggered and maintained by cloud
radiative feedback in conditions appropriate for BDs and directly imaged EGPs.
When the bottom frictional drag is strong, zonal flows in the deep layers of
our models are weak (with speeds much smaller than 100$\;{\rm m\;s^{-1}}$). In
the observable layer where clouds form, mid-to-high latitudes are dominated by
isotropic vortices, with thick clouds forming in anticyclones and thin clouds
or clear sky in cyclones. This is consistent with the results of previous
$f$-plane models (Tan & Showman, 2020). At low latitudes, large-scale zonally
propagating waves with both eastward and westward phase speeds are the
dominant dynamical feature. At a given rotation period, sizes of storms and
vortices are typically larger at lower latitudes than those at higher
latitudes. The overall sizes of storms are larger when the rotation period is
longer.

  * •

Lightcurves from our simulations have amplitudes from 0.5 percent to several
percent, consistent with the majority of observed lightcurves. For a given
rotation period, the lightcurve amplitude decreases with increasing viewing
angle (0∘ means equator-on and 90∘ means pole-on), while it typically
increases with increasing rotation period at a given viewing angle. When the
bottom drag is strong, zonally propagating waves at low latitudes have typical
phase speeds of a few hundred $\;{\rm m\;s^{-1}}$. They can cause short-term
evolution of lightcurves via wave beating effects, qualitatively similar to
the observed lightcurves of some field BDs as summarized in Apai et al.
(2017). The eastward Kelvin waves can cause the equatorial flux inhomogeneity
rotating faster than the underlying planetary rotation, which may help to
explain the observed shorter rotation period of atmospheric features than that
of the interior (Allers et al., 2020). Isotropic storms and vortices at mid-
to-high latitudes also contribute to the lightcurve variability.

  * •

Clouds are generally mixed to higher altitudes near the equator than at high
latitudes due to the stronger effect of rotation at high latitudes. This
supports the observed redder IR colors for objects viewed more equator on (Vos
et al., 2017; Vos et al., 2020).

  * •

We expect that the strong-drag models might be appropriate for BDs and
directly imaged EGPs because efficient convective mixing in the interior is
expected to suppress strong zonal flows near the top of convective zone. But
we still perform experiments with weaker drags to explore dynamics in the
weak-drag regime. We find that robust zonal jets with speeds from several
hundred to more than 2000$\;{\rm m\;s^{-1}}$ can form in our weak-drag models,
with typically a broad westward equatorial jet and a high-latitude eastward
jet in each hemisphere. Similar to the weak-drag cases, vortices form at mid-
to-high latitudes and equatorially trapped waves form at low latitudes. Both
the zonal propagation of equatorial waves and the spacial distribution of
vortices are affected by the presence of strong jets. Simulated lightcurves
show longer periodicity than the underlying rotation period due to the strong
westward equatorial jets. The meridionally broad jet structure may be related
to the efficient potential vorticity mixing associated with intense eddies.

  * •

The origin of the equatorially propagating waves in our simulations is likely
related to the self-excitation by cloud radiative feedbacks. Physical
properties of the equatorially symmetric eastward waves resembles properties
of both adiabatic free Kelvin waves and forced-damped waves triggered by a
stationary equatorial heat source. Linear stability theory of waves coupled
with cloud radiative feedback may help to explain the origin of our simulated
waves but fails to explain their propagation.

## Acknowledgements

We thank Tad Komacek and Xi Zhang for comments on the draft and Ray
Pierrehumbert for discussion. X.T. acknowledges support from the European
community through the ERC advanced grant EXOCONDENSE (PI: R.T. Pierrehumbert).
This work was completed with resources provided by the department of Physics
at University of Oxford and the Lunar and Planetary Laboratory at University
of Arizona.

## Data availability

The data underlying this article will be shared on reasonable request to the
corresponding author.

## References

  * Ackerman & Marley (2001) Ackerman A. S., Marley M. S., 2001, ApJ, 556, 872
  * Adcroft et al. (2004) Adcroft A., Campin J.-M., Hill C., Marshall J., 2004, Monthly Weather Review, 132, 2845
  * Allers et al. (2020) Allers K. N., Vos J. M., Biller B. A., Williams P. K., 2020, Science, 368, 169
  * Apai et al. (2013) Apai D., Radigan J., Buenzli E., Burrows A., Reid I. N., Jayawardhana R., 2013, ApJ, 768, 121
  * Apai et al. (2017) Apai D., et al., 2017, Science, 357, 683
  * Artigau (2018) Artigau É., 2018, preprint,  (arXiv:1803.07672)
  * Artigau et al. (2009) Artigau É., Bouchard S., Doyon R., Lafrenière D., 2009, The Astrophysical Journal, 701, 1534
  * Barman et al. (2011) Barman T. S., Macintosh B., Konopacky Q. M., Marois C., 2011, ApJ, 733, 65
  * Biller (2017) Biller B., 2017, Astronomical Review, 13, 1
  * Biller et al. (2013) Biller B. A., et al., 2013, ApJ, 778, L10
  * Biller et al. (2015) Biller B. A., et al., 2015, The Astrophysical Journal Letters, 813, L23
  * Biller et al. (2018) Biller B. A., et al., 2018, AJ, 155, 95
  * Bowler et al. (2020) Bowler B. P., Zhou Y., Morley C. V., Kataria T., Bryan M. L., Benneke B., Batygin K., 2020, ApJ, 893, L30
  * Buenzli et al. (2012) Buenzli E., et al., 2012, ApJ, 760, L31
  * Buenzli et al. (2014) Buenzli E., Apai D., Radigan J., Reid I. N., Flateau D., 2014, The Astrophysical Journal, 782, 77
  * Buenzli et al. (2015) Buenzli E., Marley M. S., Apai D., Saumon D., Biller B. A., Crossfield I. J., Radigan J., 2015, The Astrophysical Journal, 812, 163
  * Burningham et al. (2017) Burningham B., Marley M. S., Line M. R., Lupu R., Visscher C., Morley C., Saumon D., Freedman R., 2017, Monthly Notices of the Royal Astronomical Society, p. stx1246
  * Burrows et al. (2001) Burrows A., Hubbard W. B., Lunine J. I., Liebert J., 2001, Reviews of Modern Physics, 73, 719
  * Burrows et al. (2006) Burrows A., Sudarsky D., Hubeny I., 2006, ApJ, 640, 1063
  * Carone et al. (2020) Carone L., et al., 2020, Monthly Notices of the Royal Astronomical Society, 496, 3582
  * Chauvin et al. (2017) Chauvin G., et al., 2017, Astronomy & Astrophysics, 605, L9
  * Clarke et al. (2008) Clarke F. J., Hodgkin S. T., Oppenheimer B. R., Robertson J., Haubois X., 2008, MNRAS, 386, 2009
  * Crossfield et al. (2014) Crossfield I. J. M., et al., 2014, Nature, 505, 654
  * Currie et al. (2011) Currie T., et al., 2011, The Astrophysical Journal, 729, 128
  * Del Genio et al. (2009) Del Genio A. D., Achterberg R. K., Baines K. H., Flasar F. M., Read P. L., Sánchez-Lavega A., Showman A. P., 2009, in , Saturn from Cassini-Huygens. Springer, pp 113–159
  * Dritschel & McIntyre (2008) Dritschel D., McIntyre M., 2008, Journal of the Atmospheric Sciences, 65, 855
  * Dunkerton & Scott (2008) Dunkerton T. J., Scott R. K., 2008, Journal of the Atmospheric Sciences, 65, 1105
  * Emanuel et al. (1994) Emanuel K. A., David Neelin J., Bretherton C. S., 1994, Quarterly Journal of the Royal Meteorological Society, 120, 1111
  * Eriksson et al. (2019) Eriksson S. C., Janson M., Calissendorff P., 2019, Astronomy & Astrophysics, 629, A145
  * Faherty et al. (2016) Faherty J. K., et al., 2016, The Astrophysical Journal Supplement Series, 225, 10
  * Freedman et al. (2014) Freedman R. S., Lustig-Yaeger J., Fortney J. J., Lupu R. E., Marley M. S., Lodders K., 2014, The Astrophysical Journal Supplement Series, 214, 25
  * Freytag et al. (2010) Freytag B., Allard F., Ludwig H.-G., Homeier D., Steffen M., 2010, A&A, 513, A19
  * Gelino et al. (2002) Gelino C. R., Marley M. S., Holtzman J. A., Ackerman A. S., Lodders K., 2002, The Astrophysical Journal, 577, 433
  * Gierasch et al. (1973) Gierasch P. J., Ingersoll A. P., Williams R. T., 1973, Icarus, 19, 473
  * Gill (1980) Gill A. E., 1980, Quarterly Journal of the Royal Meteorological Society, 106, 447
  * Gillon et al. (2013) Gillon M., Triaud A., Jehin E., Delrez L., Opitom C., Magain P., Lendl M., Queloz D., 2013, Astronomy & Astrophysics, 555, L5
  * Guillot (2010) Guillot T., 2010, A&A, 520, A27
  * Guillot et al. (2018) Guillot T., et al., 2018, Nature, 555, 227
  * Haertel & Kiladis (2004) Haertel P. T., Kiladis G. N., 2004, Journal of the atmospheric sciences, 61, 2707
  * Hayashi (1970) Hayashi Y., 1970, Journal of the Meteorological Society of Japan. Ser. II, 48, 140
  * Helling (2019) Helling C., 2019, Annual Review of Earth and Planetary Sciences, 47, 583
  * Helling & Fomins (2013) Helling C., Fomins A., 2013, Philosophical Transactions of the Royal Society of London A: Mathematical, Physical and Engineering Sciences, 371, 20110581
  * Hitchcock et al. (2020) Hitchcock J. A., et al., 2020, arXiv e-prints, p. arXiv:2005.06906
  * Holton & Hakim (2012) Holton J. R., Hakim G. J., 2012, An introduction to dynamic meteorology. Vol. 88, Academic press
  * Iess et al. (2019) Iess L., et al., 2019, Science, 364
  * Imamura et al. (2020) Imamura T., Mitchell J., Lebonnois S., Kaspi Y., Showman A. P., Korablev O., 2020, Space Sci. Rev., 216, 87
  * Ingersoll et al. (2004) Ingersoll A. P., et al., 2004
  * Kao et al. (2016) Kao M. M., Hallinan G., Pineda J. S., Escala I., Burgasser A., Bourke S., Stevenson D., 2016, The Astrophysical Journal, 818, 24
  * Kao et al. (2018) Kao M. M., Hallinan G., Pineda J. S., Stevenson D., Burgasser A., 2018, The Astrophysical Journal Supplement Series, 237, 25
  * Kaspi et al. (2018) Kaspi Y., et al., 2018, Nature, 555, 223
  * Kiladis et al. (2009) Kiladis G. N., Wheeler M. C., Haertel P. T., Straub K. H., Roundy P. E., 2009, Reviews of Geophysics, 47
  * Leggett et al. (2016) Leggett S. K., et al., 2016, ApJ, 830, 141
  * Lew et al. (2016) Lew B. W., et al., 2016, The Astrophysical Journal Letters, 829, L32
  * Lew et al. (2020) Lew B. W., et al., 2020, The Astronomical Journal, 159, 125
  * Liu & Showman (2013) Liu B., Showman A. P., 2013, ApJ, 770, 42
  * Liu et al. (2016) Liu M. C., Dupuy T. J., Allers K. N., 2016, The Astrophysical Journal, 833, 96
  * Manjavacas et al. (2017) Manjavacas E., et al., 2017, The Astronomical Journal, 155, 11
  * Manjavacas et al. (2019) Manjavacas E., et al., 2019, The Astrophysical Journal Letters, 875, L15
  * Marley & Robinson (2015) Marley M. S., Robinson T. D., 2015, Annual Review of Astronomy and Astrophysics, 53, 279
  * Marley & Sengupta (2011) Marley M. S., Sengupta S., 2011, MNRAS, 417, 2874
  * Marley et al. (2012) Marley M. S., Saumon D., Cushing M., Ackerman A. S., Fortney J. J., Freedman R., 2012, ApJ, 754, 135
  * Matsuno (1966) Matsuno T., 1966, Journal of the Meteorological Society of Japan. Ser. II, 44, 25
  * Metchev et al. (2015) Metchev S. A., et al., 2015, ApJ, 799, 154
  * Miles-Páez et al. (2017) Miles-Páez P. A., Pallé E., Zapatero Osorio M. R., 2017, MNRAS, 472, 2297
  * Miles-Páez et al. (2019) Miles-Páez P. A., et al., 2019, The Astrophysical Journal, 883, 181
  * Millar-Blanchaer et al. (2020) Millar-Blanchaer M. A., et al., 2020, The Astrophysical Journal, 894, 42
  * Morley et al. (2014) Morley C. V., Marley M. S., Fortney J. J., Lupu R., 2014, The Astrophysical Journal Letters, 789, L14
  * Neelin & Held (1987) Neelin J. D., Held I. M., 1987, Monthly Weather Review, 115, 3
  * O’Neill et al. (2015) O’Neill M. E., Emanuel K. A., Flierl G. R., 2015, Nature Geoscience, 8, 523
  * Radigan (2014) Radigan J., 2014, The Astrophysical Journal, 797, 120
  * Radigan et al. (2012) Radigan J., Jayawardhana R., Lafrenière D., Artigau É., Marley M., Saumon D., 2012, ApJ, 750, 105
  * Radigan et al. (2014) Radigan J., Lafrenière D., Jayawardhana R., Artigau E., 2014, The Astrophysical Journal, 793, 75
  * Rajan et al. (2015) Rajan A., et al., 2015, Monthly Notices of the Royal Astronomical Society, 448, 3775
  * Rhines (1975) Rhines P. B., 1975, Journal of Fluid Mechanics, 69, 417
  * Robinson & Marley (2014) Robinson T. D., Marley M. S., 2014, The Astrophysical Journal, 785, 158
  * Saumon & Marley (2008) Saumon D., Marley M. S., 2008, The Astrophysical Journal, 689, 1327
  * Schneider & Liu (2009) Schneider T., Liu J., 2009, Journal of Atmospheric Sciences, 66, 579
  * Schneider et al. (2018) Schneider A. C., Hardegree-Ullman K. K., Cushing M. C., Kirkpatrick J. D., Shkolnik E. L., 2018, The Astronomical Journal, 155, 238
  * Scott (2011) Scott R., 2011, Geophysical & Astrophysical Fluid Dynamics, 105, 409
  * Scott & Dritschel (2012) Scott R. K., Dritschel D. G., 2012, Journal of Fluid Mechanics, 711, 576
  * Sengupta & Marley (2009) Sengupta S., Marley M. S., 2009, ApJ, 707, 716
  * Shapiro (1970) Shapiro R., 1970, Reviews of geophysics, 8, 359
  * Showman & Kaspi (2013) Showman A. P., Kaspi Y., 2013, ApJ, 776, 85
  * Showman & Polvani (2011) Showman A. P., Polvani L. M., 2011, The Astrophysical Journal, 738, 71
  * Showman et al. (2010) Showman A. P., Cho J. Y.-K., Menou K., 2010, Atmospheric Circulation of Exoplanets. pp 471–516
  * Showman et al. (2018) Showman A. P., Ingersoll A. P., Achterberg R., Kaspi Y., 2018, Saturn in the 21st Century, 20, 295
  * Showman et al. (2019) Showman A. P., Tan X., Zhang X., 2019, ApJ, 883, 4
  * Showman et al. (2020) Showman A. P., Tan X., Parmentier V., 2020, Space Sci. Rev., 216, 139
  * Stolker et al. (2020) Stolker T., et al., 2020, Astronomy & Astrophysics, 635, A182
  * Tan & Komacek (2019) Tan X., Komacek T. D., 2019, The Astrophysical Journal, 886, 26
  * Tan & Showman (2017) Tan X., Showman A. P., 2017, ApJ, 835, 186
  * Tan & Showman (2019) Tan X., Showman A. P., 2019, The Astrophysical Journal, 874, 111
  * Tan & Showman (2020) Tan X., Showman A. P., 2020, arXiv e-prints, p. arXiv:2005.12152
  * Thomson & McIntyre (2016) Thomson S. I., McIntyre M. E., 2016, Journal of the Atmospheric Sciences, 73, 1119
  * Tremblin et al. (2020) Tremblin P., Phillips M., Emery A., Baraffe I., Lew B., Apai D., Biller B., Bonnefoy M., 2020, Astronomy & Astrophysics, 643, A23
  * Tsai et al. (2014) Tsai S.-M., Dobbs-Dixon I., Gu P.-G., 2014, ApJ, 793, 141
  * Tsuji (2002) Tsuji T., 2002, ApJ, 575, 264
  * Vallis & Maltrud (1993) Vallis G. K., Maltrud M. E., 1993, Journal of physical oceanography, 23, 1346
  * Vallis & Penn (2020) Vallis G. K., Penn J., 2020, Quarterly Journal of the Royal Meteorological Society
  * Vasavada & Showman (2005) Vasavada A. R., Showman A. P., 2005, Reports on Progress in Physics, 68, 1935
  * Vos et al. (2017) Vos J. M., Allers K. N., Biller B. A., 2017, The Astrophysical Journal, 842, 78
  * Vos et al. (2018) Vos J. M., Allers K. N., Biller B. A., Liu M. C., Dupuy T. J., Gallimore J. F., Adenuga I. J., Best W. M., 2018, Monthly Notices of the Royal Astronomical Society, 474, 1041
  * Vos et al. (2019) Vos J. M., et al., 2019, Monthly Notices of the Royal Astronomical Society, 483, 480
  * Vos et al. (2020) Vos J. M., et al., 2020, AJ, 160, 38
  * Wheeler & Kiladis (1999) Wheeler M., Kiladis G. N., 1999, Journal of the Atmospheric Sciences, 56, 374
  * Wilson et al. (2014) Wilson P. A., Rajan A., Patience J., 2014, A&A, 566, A111
  * Yang et al. (2015) Yang H., et al., 2015, ApJ, 798, L13
  * Yang et al. (2016) Yang H., et al., 2016, ApJ, 826, 8
  * Zhang (2020) Zhang X., 2020, Research in Astronomy and Astrophysics, 20, 099
  * Zhang & Showman (2014) Zhang X., Showman A. P., 2014, ApJ, 788, L6
  * Zhou et al. (2016) Zhou Y., Apai D., Schneider G. H., Marley M. S., Showman A. P., 2016, ApJ, 818, 176
  * Zhou et al. (2018) Zhou Y., et al., 2018, The Astronomical Journal, 155, 132
  * Zhou et al. (2020) Zhou Y., Bowler B. P., Morley C. V., Apai D., Kataria T., Bryan M. L., Benneke B., 2020, arXiv e-prints, p. arXiv:2004.05168

##  Appendix A Hovm$\ddot{\rm o}$ller diagrams at mid latitudes

In Figure 17, we show the Hovm$\ddot{\rm o}$ller diagrams of outgoing thermal
flux at $45^{\circ}$ latitude as a function of longitude and time for models
with different rotation period. These regions are dominated by stochastically
evolving vortices, and there is no systematic eastward or westward propagation
seen in these diagrams, which is in stark contrast to the equatorial
disturbances shown in Figure 5.

Figure 17: Hovm$\ddot{\rm o}$ller diagrams (longitude-time sequence) showing
the time evolution of outgoing thermal flux at $45^{\circ}$ latitude for
models with four rotation periods of 20, 10, 5 and 2.5 hours (from the left to
right) and a drag timescale $\tau_{\rm{drag}}=10^{5}$ s.

##  Appendix B Cloud radiatively induced equatorial waves

Equatorially trapped waves in our simulations play vital roles in the
evolution of simulated lightcurves, and may be responsible for the puzzling
time evolution of some observed lightcurves (Apai et al., 2017). Given their
central importance, we explore their nature using a linear wave theory that is
coupled to cloud radiative effect. We will show that the linear theory may be
used to explain the energetic sources of our simulated waves, but cannot
capture the propagating nature of the simulated waves.

The analysis of cloud radiatively induced instability was first carried out by
Gierasch et al. (1973) in the quasi-geostrophic system and inertia gravity
wave ($f$-plane) system. In the absence of rotation, 2D hydrostatic gravity
waves have a set of pure unstable growing modes (no propagation) and sets of
attenuating, eastward and westward propagating mode. In a $f-$plane
approximation, inertia gravity waves with small horizontal wavelengths can
host purely unstable modes but unstable modes cease at larger horizontal
scales. In a quasi-geostrophic flow, unstable modes are possible for axially
symmetric flows and flows with nonzero $\beta$ (Gierasch et al., 1973). We now
extend this theory to the equatorial $\beta$ plane.

###  B.1 Cloud radiative effect upon thermodynamics

We first present the thermodynamics related to cloud radiative effect
following Gierasch et al. (1973). Assuming that the change of outgoing thermal
flux is due to brightness temperature deviation that is caused by either
actual temperature variation or the cloud-top altitude variation, we can write
the change of emitted flux as

$\delta(\sigma T_{c}^{4})=4\sigma T_{c}^{3}(\delta T_{c}-\Gamma_{c}\delta
z_{c}),$ (3)

where $T_{c}$ is the brightness temperature, $\sigma$ is the Stefan-Boltzmann
constant, $\delta z_{c}$ is the change of cloud-top altitude, and
$\Gamma_{c}=|\frac{dT_{c}}{dz_{c}}$|. Because of the change of outgoing
radiative flux, the atmospheric column is no longer in equilibrium and so net
heating/cooling must occur. Denoting $Q$ as the heating rate per unit mass,
$\rho$ the gas density and $D$ the cloud thickness, we have

$\int_{0}^{D}Q\rho dz=-\delta(\sigma T_{c}^{4})=4\sigma
T_{c}^{3}(\Gamma_{c}\delta z_{c}-\delta T_{c}).$ (4)

For analytic simplicity, $Q$ is assumed to obey the following linear relation
at any level:

$Q=\frac{4\sigma T_{c}^{3}}{M}(\Gamma_{c}\delta z-\delta T),$ (5)

where $M=\int_{0}^{D}\rho dz$. Moreover, clouds are assumed to be perfectly
advected by the flow, i.e., sedimentation is negligible compared to the
vertical advection, so that the rate of change in cloud-top altitude is simply
the vertical velocity $w$. Taking the time derivative of Eq. (5), we have

$\frac{\partial Q}{\partial t}=\frac{4\sigma
T_{c}^{3}}{M}\left(\Gamma_{c}w-\frac{\partial T}{\partial t}\right).$ (6)

Then we assume a basic state of the atmosphere at rest with small horizontal
temperature perturbations. The linearized equation of thermodynamics in height
coordinates expressed using temperature is written:

$\frac{\partial T}{\partial t}+w\left(\frac{\partial\overline{T}}{\partial
z}+\frac{g}{c_{p}}\right)=\frac{Q}{c_{p}}.$ (7)

Combining with Eq. (6) we have the linearized thermodynamic equation

$\frac{\partial^{2}T}{\partial t^{2}}+\frac{\partial w}{\partial
t}\overline{\Gamma}=\frac{1}{\tau}\left(\Gamma_{c}w-\frac{\partial T}{\partial
t}\right),$ (8)

where $\overline{\Gamma}=\frac{d\overline{T}}{dz}+\frac{g}{c_{p}}$, and
$\frac{1}{\tau}=\frac{4\sigma T_{c}^{3}}{c_{p}M}$. Here $\tau$ represents a
radiative timescale of the cloud-forming atmosphere.

By just the terms involving vertical velocity in Equation (8), $\frac{\partial
w}{\partial t}\overline{\Gamma}=\frac{1}{\tau}\Gamma_{c}w$, an instability is
obvious with a growth rate of
$\lambda_{c}=\frac{\Gamma_{c}}{\overline{\Gamma\tau}}$. Now we consider the
full linearized dynamical equations below.

###  B.2 Rossby, mixed Rossby-gravity and inertia-gravity modes

We start with the linearized dynamical equations with a basic state at rest in
an equatorial $\beta$ plane, where the Coriolis parameter is written $f=\beta
y$. The zonal and meridional angular momentum , hydrostatic balance, and
continuity equations are, respectively

$\frac{\partial u}{\partial t}-\beta yv=-\frac{\partial\phi}{\partial x},$ (9)
$\frac{\partial v}{\partial t}+\beta yu=-\frac{\partial\phi}{\partial y},$
(10) $\frac{\partial\phi}{\partial z}=\frac{gT}{T_{0}},$ (11) $\frac{\partial
u}{\partial x}+\frac{\partial v}{\partial y}+\left(\frac{\partial}{\partial
z}-\frac{1}{H}\right)\omega=0,$ (12)

where $T_{0}$ is a reference temperature, $z=-H\ln(p/p_{0})$ is the vertical
coordinate at log pressure, $H=RT_{0}/g$ is the scale height, and $\phi$ is
the geopotential. Equations (8) to (12) forms a closed set.

We seek wave-like solutions:

$\displaystyle\\{u,v,\phi,T,\omega\\}=$ (13)
$\displaystyle\\{\tilde{u}(y),\tilde{v}(y),\tilde{\phi}(y),\tilde{T}(y),\tilde{\omega}(y)\\}e^{z/(2H)}\sin(\pi
z/D)e^{ik_{x}x+\lambda t},$

where $\tilde{u}(y),\tilde{v}(y),\tilde{\phi}(y),\tilde{T}(y)$ and
$\tilde{\omega}(y)$ are functions of $y$ only. To satisfy the form of
equatorially trapped waves, the boundary condition is typically applied to
ensure that the disturbances vanish when $|y|\rightarrow\infty$. Repeating the
derivations for the classic equatorial waves (e.g., Matsuno, 1966; Holton &
Hakim, 2012), we arrive at an equation for $\tilde{v}(y)$:

$\left[\mathcal{L}^{2}_{d}\left(\frac{\lambda^{2}(\frac{1}{\tau}+\lambda)}{c^{2}_{g}(\lambda_{c}-\lambda)}-k^{2}_{x}+\frac{ik_{x}\beta}{\lambda}\right)-Y^{2}\right]\tilde{v}+\frac{d^{2}\tilde{v}}{dY^{2}}=0,$
(14)

where $c^{2}_{g}=N^{2}/k^{2}_{z}$, $N^{2}=g\overline{\Gamma}/T_{0}$,
$Y=y/\mathcal{L}_{d}$, and $\mathcal{L}_{d}$ is the equatorial deformation
radius modified by the diabatic cloud radiative effect:

$\mathcal{L}^{4}_{d}\equiv\frac{c^{2}_{g}(\lambda-\lambda_{c})}{(\frac{1}{\tau}+\lambda)\beta^{2}}.$
(15)

Note that in the equatorial theory of Hayashi (1970) that considered latent
heating effects, a complex equatorial deformation was also possible. In the
adiabatic limit of $\tau\rightarrow\infty$, the deformation radius recovers
the classic definition $\mathcal{L}_{d}=\sqrt{c_{g}/\beta}$ and the dispersion
recovers the classic dispersion relation of adiabatic, unforced equatorially
trapped waves (Matsuno, 1966). In the limit of $\tau\rightarrow 0$, there is
no solution that satisfies the boundary condition. The following relation has
to be met for the given boundary condition:

$\mathcal{L}^{2}_{d}\left(\frac{\lambda^{2}(\frac{1}{\tau}+\lambda)}{c^{2}_{g}(\lambda_{c}-\lambda)}-k^{2}_{x}+\frac{ik_{x}\beta}{\lambda}\right)=2n+1;\quad\quad
n=0,1,2,...$ (16)

Then, the solution has the form

$\tilde{v}(Y)=\tilde{v}_{0}{\rm H}_{n}(Y)\exp(-Y^{2}/2),$ (17)

where $\tilde{v}_{0}$ is a constant with units of velocity, and ${\rm
H}_{n}(Y)$ designates the $n$th Hermite polynomial. In addition, solutions
satisfy the boundary condition if $\exp(-Y^{2}/2)$ diminishes when
$|y|\rightarrow\infty$, which requires that the real part of
$\mathcal{L}_{d}^{2}$ is positive.

We first seek solutions with real, positive $\lambda$. Possible solutions
should be in the range between 0 and $\lambda_{c}$. Outside this range, there
is no solution that satisfies a real $k_{x}$, which is required to have wave-
like zonal disturbances. Equating the imaginary part in the LHS of Equation
(16) to zero, one obtains

$\lambda^{3}+\frac{1}{\tau}\lambda^{2}+c^{2}_{g}k^{2}_{x}\lambda-\lambda_{c}c^{2}_{g}k^{2}_{x}=0.$
(18)

There is one positive real root that satisfies $0<\lambda<\lambda_{c}$. This
dispersion relation is that of the pure hydrostatic gravity waves, the same as
Equation (28) in Gierasch et al. (1973) when $f=0$. Equating the remaining
real parts in Equation (16), one obtains

$(2n+1)^{2}\left(\lambda^{3}+\frac{1}{\tau}\lambda^{2}\right)+c^{2}_{g}k^{2}_{x}\lambda-\lambda_{c}c^{2}_{g}k^{2}_{x}=0.$
(19)

Only when $n=0$ (the mixed Rossby gravity modes), both equations (18) and (19)
can be simultaneously satisfied, and there is a set of purely unstable, non-
propagating $n=0$ modes. In this case, $\mathcal{L}_{d}^{2}$ is purely
imaginary and fails the strict meridional boundary condition. Nevertheless, by
solving equation (16) numerically, we cannot find other modes that contain a
positive real part of $\lambda$. All other modes have both negative real parts
and imaginary parts, and they resemble the classic equatorial waves—the
propagating Rossby, MRG and inertia-gravity modes but with damping in the wave
amplitudes due to thermal radiation. The dispersion relations for all possible
solutions are plotted as solid curves in Figure 18, in which the top panel
shows the $-\lambda_{i}$ (in the same format as that plotted in classic
literature of equatorial waves, e.g., Holton & Hakim, 2012) and the bottom
panel shows the $\lambda_{r}$, where $\lambda=\lambda_{r}+i\lambda_{i}$.
Despite that the special non-propagating unstable mode in Equation (18) does
not strictly satisfy the lateral boundary condition, at least their
disturbances do not amplify with $|y|\rightarrow\infty$. So we still plot the
dispersion relation as the dotted lines in Figure 18.

Figure 18:  Top panel: dispersion relations of the imaginary components
$-\lambda_{i}$, i.e., the propagating components, for cloud-radiatively-
coupled free equatorial waves. Frequency and zonal wavenumbers have been non-
dimensionalized as $\lambda^{\ast}=\lambda/(\beta c_{g})^{1/2}$ and
$k^{\ast}=k(c_{g}/\beta)^{1/2}$. Different modes are labelled in the right
panel and with different colors and line styles. The dotted line corresponds
to a special case in the Kelvin and $n=0$ unstable modes. Bottom panel:
dispersion relations of the real components $\lambda_{r}$, i.e., the
attenuating or growing components, for equatorial waves. Negative
$\lambda_{r}$ represents attenuating modes and positive $\lambda_{r}$
represents growing modes.

###  B.3 Kelvin modes

The Kelvin mode is a special mode in which the meridional velocity is zero.
The governing equations for the zonal and meridional angular momentum, and
continuity are simplified to

$\frac{\partial u}{\partial t}=-\frac{\partial\phi}{\partial x},$ (20) $\beta
yu=-\frac{\partial\phi}{\partial y},$ (21) $\frac{\partial u}{\partial
x}+\left(\frac{\partial}{\partial z}-\frac{1}{H}\right)\omega=0.$ (22)

We make use of equations (20), (11), (22) and (8) and assume wave-like
solutions of equation (13), then we obtain the following dispersion relation

$\lambda^{3}+\frac{1}{\tau}\lambda^{2}+c^{2}_{g}k^{2}_{x}\lambda-\lambda_{c}c^{2}_{g}k^{2}_{x}=0.$
(23)

This is the same as Equation (18) for $n=0$ modes, and growing modes with
positive $\lambda$ exist. Additional constraints from the lateral boundary
condition should be satisfied. Combining equations (20) and (21), one obtains

$\tilde{u}=\tilde{u}_{0}\exp\left(\beta\frac{k_{x}\lambda_{i}+ik_{x}\lambda_{r}}{2(\lambda_{i}^{2}+\lambda_{r}^{2})}y^{2}\right),$
(24)

where $\tilde{u}_{0}$ is the amplitude of the perturbation zonal velocity at
the equator, and $\lambda$ is written as $\lambda=\lambda_{r}+i\lambda_{i}$.
If one restricts $\tilde{u}$ to vanish with $|y|\rightarrow\infty$,
$k_{x}\lambda_{i}<0$ needs to be satisfied. With this regard, the purely
unstable modes fail to satisfy the boundary condition, and only the eastward
propagating but decaying modes are valid solutions. This reaches the same
conclusion as the $n=0$ unstable modes. The dispersion relation of the
unstable Kelvin modes are also represented as dotted lines in Figure 18.

Our key finding is that there is only one set of growing but non-propagating
modes corresponding to the Kelvin and $n=0$ modes that may be marginally
relevant. This may provide a way to excite the Kelvin waves, $n=0$ MRG waves
and $n=0$ eastward inertia gravity waves seen in our simulations. Other
propagating modes have properties quantitatively similar to the adiabatic free
modes but with damping of their amplitudes due to thermal radiation.