1 Introduction

Giant planets, broadly defined, span an enormous range of objects. Limiting ourselves to substellar bodies comprised primarily of hydrogen,Footnote 1 such bodies nevertheless encompass an impressive diversity. Jupiter and Saturn represent the canonical prototypes, and of course are the best observed due to their proximity to Earth. Outside our solar system, hundreds of extrasolar giant planets (EGPs) have been discovered. Hot Jupiters are the most easily observationally characterized; they orbit extremely close to their stars, at distances of typically \(\sim 0.03\mbox{--}0.1~\mbox{AU}\), receive thousands of times more starlight than Jupiter, and thereby achieve temperatures of 1000 K or more (Showman and Guillot 2002). At dayside temperatures exceeding \(\sim 2200~\mbox{K}\), the molecular constituents of their atmospheres start to dissociate and they are called ultra hot Jupiters (Bell and Cowan 2018; Parmentier et al. 2018). Also amenable to atmospheric characterization are the directly imaged planets—that is, planets that are sufficiently hot and distant from their host stars to be imaged as distinct entities. They are hot not because they are strongly irradiated, but because they are massive and young, so that they still glow from their heat of formation—and therefore also have temperatures of typically \(\sim1000~\mbox{K}\) (Bowler 2016). Intermediate between these extremes are a large population of EGPs that are irradiated, but less so than hot Jupiters, and also which are old enough to have lost much of their internal heat of formation (Guillot et al. 1996); as a result, they exhibit cooler temperatures and are harder to observe. One might usefully define “warm Jupiters” to be those objects with temperatures of 300 to 1000 K (corresponding to distances from a sunlike star of approximately 1 to 0.1 AU) and “cool Jupiters” to be those objects with temperatures less than 300 K (corresponding to orbital distances from a sunlike star exceeding 1 AU). Although the warm and cool Jupiters are harder to observe than hot Jupiters and directly imaged planets, far more have been discovered, and they will be increasingly amenable to observational characterization in the future when more sensitive instruments will become available.

Brown dwarfs are objects thought to have formed like stars but which have insufficient mass to fuse hydrogen (Chabrier et al. 2000; Burrows et al. 2001); they are typically defined as objects of \(\sim10\) to 80 Jupiter masses (the stellar mass limit). Lacking strong internal thermonuclear heat generation, they cool off over time, but their super-Jovian mass implies that even after billions of years they may still exhibit atmospheric temperatures of 1000 K or more (e.g., Burrows et al. 2001). Typically, they are isolated objects, far from any star, which makes them easier to observe than exoplanets. Although their formation mechanisms differ from EGPs, brown dwarfs share many physical similarities to the currently known directly imaged planets; from an atmospheric dynamics point of view, the former may be considered high-mass, high-gravity versions of the latter.

Although far less numerous than known EGPs and field brown dwarfs, a population of brown-dwarf companions to stars has also been discovered. Some such objects orbit sunlike stars in tight, several-day orbits and resemble high-mass, high-gravity versions of hot Jupiters (see a summary in Bayliss et al. 2016). Other brown dwarfs orbit white dwarf stars so closely that the two objects nearly touch, and have orbital periods of just several hours (see a summary in Casewell et al. 2015). In some such systems—the cataclysmic variables—the brown dwarf lies so close to the white dwarf that it continually sheds mass onto the white dwarf, while in other systems, the orbital separation (while still tight) is great enough to prevent mass exchange.

Despite this broad diversity, there is merit in considering these objects together as a class. They share in common the fact that they are all fluid, hydrogen-dominated objects; they have radii similar to Jupiter to within a factor of two; their atmospheres merges continuously into their interiors; and because of the large opacity and low viscosity of hydrogen in the conditions of their interiors, they are generally expected to lose their internal heat by convection, implying convective well-mixed interiors.Footnote 2 Like all planetary atmospheres, giant planets and brown dwarfs should share in common many of the fundamental dynamical, physical, and chemical processes that shape the structure and circulation of atmospheres generally. Considering these objects together therefore provides a unique opportunity to better understand the physical and dynamical processes that operate in atmospheres over a wide range, and to understand how these common processes leads to diverse outcomes for different objects. This “grand challenge” is far easier to complete for giant planets than for terrestrial planets because the atmospheres of giant planets are currently amenable to observational characterization over a wide range, whereas for terrestrial planets, such observational characterization over a comparably wide range is still decades away.

The atmospheric circulation and structure of giant planets is shaped by a variety of factors, including the external irradiation (that is, the radiative flux the atmosphere receives from a nearby star), the internal convective heat flux, the planetary mass (therefore gravity), the rotation rate, and the atmospheric composition, including the overall bulk metallicity, as well as specific elemental ratios such as C/O. Even among known objects, these factors vary over ranges of \(\sim 10^{7}\), \(10^{6}\), \(10^{2}\), \(10^{2}\), and \(10^{2}\), respectively. Thus an enormous diversity is represented. Figure 1 summarizes the subpopulations where key observations constraining the atmospheric circulation have been obtained. Fortuitously, these populations span a wide range in several of the above parameters. Jupiter and Saturn exhibit weak irradiation, weak interior heat flux, and rotate rapidly (a factor that strongly influences the dynamical regime). Hot Jupiters are strongly irradiated, but likely have small interior heat fluxes (perhaps within an order of magnitude of Jupiter), and due to tidal locking, they are expected to have modest rotation rates. Field brown dwarfs and directly imaged giant planets receive negligible irradiation from their stars, but have enormous interior heat fluxes and are rapidly rotating—in many cases faster than Jupiter.

Fig. 1
figure 1

Four distinct subfields of astronomy and planetary science are providing important constraints on atmospheric structure and circulation of giant planets—hot Jupiters; solar system giant planets; brown dwarfs and directly imaged giant planets; and irradiated brown dwarfs. They together span a wide range of physical properties

To highlight the wide parameter space involved, Fig. 2 shows these diverse populations on a parameter space of external irradiation and internal heat flux, which are two of the factors that matter most for driving an atmospheric circulation, and which vary over the widest range across the known population. Jupiter, Saturn, Uranus, and Neptune occupy the “weak forcing” regime near the lower left corner, with external and interior heat fluxes that are small and comparable. Jupiter exhibits internal and irradiation fluxes of \(7.5~\mbox{W}\,\mbox{m}^{-2}\) and \(6.6~\mbox{W}\,\mbox{m}^{-2}\) (Li et al. 2018) absorbed and internal fluxes are 0.27 and \(0.70~\mbox{W}\,\mbox{m}^{-2}\) (Pearl and Conrath 1991). In contrast, hot Jupiters reside in the upper left corner, with external fluxes of \(10^{4}\) to \(10^{7}~\mbox{W}\,\mbox{m}^{-2}\), and internal fluxes that are expected to be far smaller. The hot Jupiters show us how giant planets behave when external forcing dominates. Brown dwarfs embody the opposite extreme, with enormous interior fluxes of \(10^{3}\) to \(10^{6}~\mbox{W}\,\mbox{m}^{-2}\) and typically negligible external irradiation. Lying in the lower right corner of Fig. 2, brown dwarfs yield insights on giant-planet behavior under extreme internal fluxes when external forcing is zero. Finally, the irradiated brown dwarfs represent the “strong forcing” regime in both internal and external fluxes; they occupy the upper right corner of Fig. 2.

Fig. 2
figure 2

Regime diagram of various giant planet populations as a function of interior heat flux (abscissa) and global-mean energy flux received from the star (ordinate). Solar system giant planets lie in the lower left corner, hot Jupiters in the upper left corner, and brown dwarfs in the lower right corner. A handful of highly irradiated brown dwarfs are also providing information and reside in the upper right corner. Observational constraints from these diverse populations place tight constraints on how fundamental theories for the atmospheric circulation scale with these heat fluxes over many orders of magnitude. This paper will provide a tour across this parameter space

Our grand challenge is to understand the atmosphere and interior circulation and structure on giant planets and brown dwarfs, broadly defined. Key questions include the following:

  • What is the nature of the atmospheric circulation—including the distribution and importance of zonal jets, vortices, storms, waves, and turbulence? What are the characteristic wind speeds, temperature variations, length scales, and time variability? How do they depend on parameters? These questions could be viewed as one of characterization, both from observations and careful numerical experiments.

  • How does the circulation work? What are the dynamical mechanisms controlling it? How does the interplay of these mechanisms lead to various outcomes in the behavior of the global circulation across the wide parameter space occupied by giant planets? What is the link between these mechanisms and those mechanisms well-known from study of more familiar atmospheres like Earth and Jupiter? This could be viewed as a question of understanding the behavior characterized in the first point.

  • What is the role of condensation and clouds? How critical are they to driving (or influencing) the circulation and climate, via radiative feedbacks, latent heating, re-distribution of condensable chemical species, or other mechanisms?

  • What is the role of coupling between atmospheric dynamics, radiative transfer and chemistry?

  • How do the magnetic field and the atmospheric circulation interact with each other, particularly in the hottest atmospheres?

  • Can we achieve a unified theory of giant planet atmospheric circulation that explains observations of hot Jupiters, brown dwarfs, and solar system giant planets?

  • Does this knowledge provide insights into the circulation and climate of (less easily observed) smaller planets, including habitable terrestrial planets?

Our aim is to broadly survey the atmospheric circulation across these diverse classes of giant planets and brown dwarfs. We place particular attention on dynamical mechanisms and the way they vary among these populations. This review can be viewed as a guided tour through Fig. 2, starting with hot Jupiters in the upper left corner (Sect. 2), moving downward to the warm Jupiters (Sect. 3), hopping across to the brown dwarfs in the lower-right corner (Sect. 4), and then finishing with irradiated brown dwarfs in the upper-right corner (Sect. 5). Given the existence of many prior reviews of atmospheric dynamics on Jupiter and Saturn, we touch on the dynamics of solar system planets only briefly, as points of reference with hotter giant planets. Our review updates prior reviews on the atmospheric dynamics of exoplanets (Showman et al. 2018b, 2010, 2013b; Heng and Showman 2015) and complements the many excellent reviews covering observations and radiative structure of EGPs (e.g. Deming and Seager 2009, 2017; Seager and Deming 2010; Burrows 2014; Madhusudhan et al. 2014; Madhusudhan 2019; Zhang 2020) and brown dwarfs (Stevenson 1991; Burrows et al. 2001; Helling and Casewell 2014; Marley and Robinson 2015; Biller 2017; Zhang 2020).

2 Hot Jupiters

EGPs orbiting extremely close to their stars—hot Jupiters—comprise the most easily characterizable type of exoplanet, and therefore our understanding of these exoplanets is best developed. As little as 0.03–0.05 AU from their parent star, hot Jupiters have orbital periods of just a few Earth days. The immense stellar irradiation heats their atmospheres to temperatures of \(\sim1000\mbox{--}3000~\mbox{K}\), and they therefore radiate enormous IR heat fluxes to space, promoting direct detection of their thermal emission. Close-in planets such as hot Jupiters are also more likely to transit their stars—the transit probability for a planet on a 0.05 AU, randomly inclined orbit around a sunlike star is \(\sim10\%\), versus 0.1% for a planet at Jupiter’s distance—and when such a transiting planet is detected, it enables the determination of the planet’s radius and allows atmospheric characterization through a wide suite of observation methods.

Hot Jupiters are too close to their stars to be distinctly resolvable from their stars in images—what we observe is the combined light from the planet-star system—and indirect methods are therefore needed to tease apart the planetary light from the starlight. When the planet passes behind its star—an event known as secondary eclipse—only the star is visible. Subtracting the total system flux received during secondary eclipse from that received immediately before and afterward (when both the planet and the star contribute to the combined light) yields a spectrum of the planet, and in particular of the planet’s dayside. Moreover, observing the planet throughout its orbit—as its dayside and nightside rotate in and out of view—allows the measurement of the phase variations of the planet’s outgoing IR flux, and thereby provides inferences on the longitudinal variation of temperature near the photosphere (Fig. 3). In particular, such light curve observations yield day-night temperature differences and offsets of the hot spot from the substellar point. And observations during transit—when the planet lies in front of its star—probe the atmosphere in transmission. Starlight passing through the planet’s atmosphere on its way to Earth is preferentially blocked at wavelengths where the atmosphere is more absorbing; therefore, the planet essentially appears bigger at wavelengths corresponding to atmospheric absorption, and a spectrum of the planet’s size is essentially a transmission spectrum of the planet’s atmosphere at its terminator. For recent reviews of these and other methods for characterizing EGP atmospheres, see Deming and Seager (2017) and Madhusudhan (2019).

Fig. 3
figure 3

Hot Jupiter observations and dynamical insights that can be gather from them: (a) Infrared phase curves of the hot Jupiter HD209458b observed with Spitzer at 4.5 μm. (b) Day-night temperature contrast estimated from phase curve observations (Komacek et al. 2017). (c) Thermal structure at different phases estimated based on the spectral phase curve of WASP-43b observed by HST/WFC3 (Stevenson 2016). (d) Brightness map of HD189733b reconstructed from the 8 μm phase curve observed by Spitzer (Knutson et al. 2007). (e) Doppler measurements at the planetary limb during transit allow direct detection of atmospheric winds on the different limbs for the hot Jupiter HD189733b (Louden and Wheatley 2015)

The atmospheric regime of hot Jupiters differs radically from that of any planet in the solar system. Hot Jupiters are blasted by starlight—they typically receive stellar fluxes of \(10^{5}\mbox{--}10^{6}~\mbox{W}\,\mbox{m}^{-2}\) (Fig. 2), which is \(\gtrsim 10^{4}\) times the flux received by Jupiter and hundreds to thousands of times that received by Earth. Moreover, the close-in orbital distances lead to strong tidal forces that slow their spin; for example, a typical hot Jupiter orbiting at 0.05 AU around a Sun-like star has a synchronization timescale of \(\sim10^{6}\) years, which is orders of magnitude shorter than the typical multi-Gyr system ages (Guillot et al. 1996). Therefore, hot Jupiters are generally presumed to be synchronously rotating, with permanent daysides and nightsides. This property, coupled with the intense irradiation, leads to a unique climate regime of permanent day-night forcing, which is not experienced by any planet in the solar system.

In what follows, we first sketch out the basic dynamical regime of hot Jupiters before proceeding to a detailed discussion from simulations and theory of the atmospheric circulation and the processes that maintain it. We also cover various topics of current interest including coupling of the dynamics to clouds and chemistry.

2.1 Basic Dynamical Arguments

The relative slowness of rotation (compared to Earth, Jupiter, and brown dwarfs) causes several important effects. First, it implies that the dynamical length scale on hot Jupiters will be relatively large, approaching the global scale (Showman and Guillot 2002; Menou et al. 2003; Cho et al. 2003). One measure of this effect is the equatorial deformation radius, given by

$$ L_{\mathrm{eq}}=\left (\frac{NH}{\beta }\right )^{1/2} $$
(1)

where \(N\) is the Brunt-Vaisala frequency, \(H\) is the pressure scale height, and \(\beta \) is the gradient of Coriolis parameter \(f\) with latitude, that is \(\beta = df/dy\), where \(f = 2\Omega \sin {\phi }\), \(\Omega \) is the planetary rotation rate (\(2\pi \) over the rotation period), \(\phi \) is latitude, and \(y\) is northward distance. On a sphere, \(\beta = 2\Omega /a\) at low latitudes, where \(a\) is the planetary radius. The quantity \(NH\) can be thought of as the phase speed of horizontally propagating gravity waves, which for a vertically isothermal atmosphere is just \(R\sqrt{T/c_{p}}\), where \(R\) and \(c_{p}\) are the specific gas constant and specific heat at constant pressure, and \(T\) is temperature. Under typical hot Jupiter conditions, these expressions yield an equatorial deformation radius \(L_{\mathrm{eq}} \sim 5 \times 10^{4}~\mbox{km}\)—more than half a planetary radius.Footnote 3 In contrast, for Earth and Jupiter, respectively, the equatorial deformation radius is about 30% and 11% of the planetary radius, respectively.

Thus one expects that the dynamically “tropical” conditions (i.e., where the Coriolis force is not dominant in the horizontal force balance) that prevail within a deformation radius of the equator—including equatorially trapped waves and the effect they exert in adjusting the planet’s atmosphere—will extend to significantly higher latitudes than they do on Earth and especially Jupiter. Indeed, the equatorial regions comprise a waveguide, of meridional half-width \(L_{\mathrm{eq}}\), for a wide population of tropical baroclinic waves (Andrews et al. 1987; Holton and Hakim 2012). These waves are confined to very low latitudes on Jupiter and brown dwarfs but can extend meridionally to mid-latitudes or farther on typical hot Jupiters.

Outside the equatorial waveguide, the deformation radius is given by

$$ L_{D}=\frac{NH}{f}. $$
(2)

Inserting typical numbers for a hot Jupiter yields \(L_{D} \approx 4 \times 10^{7}~\text{m}\), again about half the planetary radius.

An alternate measure of the role of rotation comes from directly assessing the amplitude of Coriolis forces relative to other forces in the equation of motion. This can be accomplished by using the Rossby number, which represents the ratio of the characteristic amplitude of advection forces per mass, \(U^{2}/L\), to Coriolis accelerations, \(fU\), in the horizontal equation of motion, where \(U\) and \(L\) are the characteristic horizontal wind speed and length scale. Taking the ratio of these forces, we have that the Rossby number is \(Ro = U/fL\). Typical wind speeds in hot Jupiters are expected to be \(\sim1\text{--}3~\text{km}\,\text{s}^{-1}\), and adopting a rotation period of 3 days and a length scale of half a planetary radius, we obtain \(Ro \approx 0.6\text{--}2\) at midlatitudes. With Rossby numbers of order unity, hot Jupiters are thus transitional between regimes where the rotation plays little role (\(Ro \gg 1\)) and where it dominates the dynamics (\(Ro \ll 1\)). As we will discuss later, the large-scale circulation away from the equator on Jupiter, Earth’s atmosphere and oceans, and probably brown dwarfs are in the latter regime, which will lead to differing dynamics between hot Jupiters and these other planets. Note that because the Coriolis parameter goes to zero at the equator, the Rossby number always tends to become large near the equator, and in fact one useful dynamical measure of the “tropics” corresponds to the range of latitudes within which \(Ro \gtrsim 1\) (e.g. Showman et al. 2013b). According to this measure, some hot Jupiters—particularly those with faster wind speeds and/or slower rotation rates—will be “all tropics” worlds where the Rossby number always exceeds one even at the poles; in contrast, other hot Jupiters with faster rotation and/or slower wind speeds may exhibit a transition where the low and midlatitudes comprise the tropics but the poles exhibit a more extratropical (\(Ro \lesssim 1\)) behavior.

The high temperatures of hot Jupiters imply that, in the observable atmosphere, they will experience short radiative time constants. Near the photosphere, the radiative time constant is approximately (Showman and Guillot 2002)

$$ \tau _{\mathrm{rad}} = \frac{pc_{p}}{4g\sigma T^{3}} \sim 10^{5}\left ( \frac{p}{0.3~\text{bar}} \right )\left (\frac{1200~\text{K}}{T}\right )^{3}~\text{s}, $$
(3)

where \(\sigma \) is the Stefan-Boltzmann constant, and \(p\) should be interpreted as a pressure near the IR photosphere. Note that, because of the \(T^{3}\) dependence, the radiative time constant varies significantly across the hot Jupiter population, from \(\sim10^{4}~\text{s}\) for ultra hot Jupiters, to \(10^{5}~\text{s}\) for intermediate-temperature planets like HD 189733b, to \(10^{6}~\text{s}\) for warm Jupiters.

What are the expected wind speeds, to order of magnitude? A simple estimate can be obtained by balancing the pressure-gradient force that drives the flow with the greater of the Coriolis force, advective forces, or drag force (if present) in the horizontal momentum equation. Generally, if the frictional drag is weak, and if \(Ro \ll 1\), then the balance is between the pressure-gradient and Coriolis force. This leads to the well-known thermal-wind equation (Vallis 2006; Holton and Hakim 2012):

$$ f\frac{\partial u}{\partial \ln p} = R\frac{\partial T}{\partial y}, $$
(4)

where \(u\) is zonal wind, \(R\) is the specific gas constant, and \(y\) is northward distance on the sphere. To order of magnitude, this equation yields

$$ U\sim \frac{R\Delta T_{\mathrm{horiz}} \delta \ln p}{fL}, $$
(5)

where \(\Delta T_{\mathrm{horiz}}\) is the characteristic horizontal temperature difference, assumed to extend vertically over a number of scale heights \(\delta \ln p\), \(L\) is the characteristic horizontal length scale, and \(U\) is the difference in characteristic horizontal wind speed between some upper level of interest (e.g., the photosphere) and deeper levels. This should be interpreted as the characteristic eddy speed (e.g., associated with the wind flow from day to night) and not the speed of equatorial superrotation. Alternately, if \(Ro \gtrsim 1\), then the pressure gradient forces are typically balanced by horizontal advection forces.Footnote 4 To order of magnitude, the former is \(R\Delta T_{\mathrm{horiz}} \delta \ln p/L\), and the latter is \(U^{2}/L\), so their balance implies a wind speed

$$ U\sim (R\Delta T_{\mathrm{horiz}} \delta \ln p)^{1/2}. $$
(6)

Since \(Ro \sim 1\) on typical hot Jupiters, the Coriolis and advection forces are comparable, and we would expect that these two expressions would yield similar estimates. Indeed, when we plug in typical numbers (e.g., \(R = 3700~\text{J}\,\text{kg}^{-1}\,\text{K}^{-1}\), \(\Delta T_{\mathrm{horiz}} \approx 400~\mbox{K}\), \(\delta \ln p \approx 3\), \(f \approx 3 \times 10^{-5}\), and adopting \(L\) of a Jupiter radius), we obtain \(U \approx 2~\text{km}\,\text{s}^{-1}\) from both estimates.

2.2 GCM Experiments and Comparison to Observations

General circulation models (GCMs) for hot Jupiters have been developed using a variety of different codes and numerical approaches. Most commonly, these solve the primitive equations of atmospheric dynamics over the full globe, assuming that the circulation is driven by the intense stellar irradiation gradient, under conditions of synchronous rotation. Consistent with the expectations described above, these models typically assume a thermal structure that is deeply stratified throughout the atmosphere. The vertical domain typically extends over many pressure scale heights from a pressure of \(\lesssim 1~\text{mbar}\) at the top, to commonly \(\sim100~\text{bars}\) at the bottom.

These GCMs predict atmospheric flows comprising several key features, including (1) eddy and jet structures of near-global scale; (2) large day-night temperature differences reaching hundreds of K; and, (3) most interestingly, a wide, fast eastward equatorial jet—so-called equatorial superrotation—which straddles the equator, extends to latitudes \(\sim 30^{\circ }\), and achieves zonal-mean zonal wind speeds of typically \(2\text{--}4~\text{km}\,\text{s}^{-1}\) (Fig. 4). Despite the synchronous rotation—which in radiative equilibrium would lead to a temperature field comprising a simple day-night temperature difference with the hottest regions at the substellar point—the dynamics distorts the temperature structure in a complex manner, most prominently by inducing an eastward displacement of the day-side hot spot by tens of degrees longitude from the sub-stellar point. The earliest GCMs capturing these general features preceded observations (Showman and Guillot 2002; Cooper and Showman 2005), and predicted that the large day-night temperature differences and eastward offsets would be detectable in IR lightcurves of these planets. This helped motivate searches for these features in IR light curves. Observations from the Spitzer Space Telescope first confirmed this prediction for the hot Jupiter HD 189733b (Knutson et al. 2007), and subsequent full-orbit IR light curve observations from the Spitzer and Hubble Space Telescopes have detected such an eastward offset on the majority of hot Jupiters that have been observed (Fig. 3; for reviews, see Heng and Showman 2015; Parmentier et al. 2015; Parmentier and Crossfield 2018).

Fig. 4
figure 4

Example GCM simulations of hot Jupiters from a variety of groups in the field. Despite differences in forcing setup, numerics, and other factors, all the models exhibit similar circulation regimes, comprising significant day-night temperature differences, a fast eastward equatorial jet, and an eastward shifted dayside hot spot. All models assume synchronously rotation and conditions appropriate for typical hot Jupiters. Top left shows observations of HD 189733b from Knutson et al. (2007). Simulations of HD 209458b are shown from Showman and Guillot (2002), Heng et al. (2011b), Rauscher and Menou (2012b), Amundsen et al. (2016), Cho et al. (2015). Simulations of HD 189733b from Showman et al. (2009). Simulations of WASP-43b from Mendonça et al. (2016). These seven simulations were performed with totally distinct numerical codes, involving varying approximation of radiative forcing, and using seven independent dynamical cores. For each image, the substellar point is in the center of the panel, except for Mendonça et al. (2018), where the antistellar point is in the center

Over the last 15 years, many additional atmospheric dynamics models have been employed to investigate the global circulation of hot Jupiters; interestingly, most of these models agree reasonably well in their qualitative predictions, including the presence of large day-night temperature differences, equatorial superrotation, and—under appropriate conditions—eastward shifted hot spots (e.g., Menou and Rauscher 2009; Rauscher and Menou 2010, 2012a, 2013; Dobbs-Dixon and Lin 2008; Dobbs-Dixon and Agol 2013; Heng et al. 2011b,a; Mayne et al. 2014, 2017; Showman et al. 2008, 2009, 2013a, 2015; Parmentier et al. 2013, 2016; Kataria et al. 2015, 2016; Lewis et al. 2010, 2013; Cho et al. 2015; Zhang and Showman 2017; Menou 2019; Tan and Komacek 2019; Mendonça 2020). Figure 4 presents representative snapshots from several different groups, which, despite the many differences in model setup, highlight the overall similarity in qualitative regime across these diverse models. The earliest hot-Jupiter GCMs adopted extremely idealized schemes to represent the day-night thermal forcing; more recent work has implemented radiative transfer schemes of varying levels of realism that represent the absorption of incoming starlight and radiation of thermal IR.

Interestingly, the qualitative properties of the circulation in these models—including the emergence of equatorial superrotation—seem to be robust to model assumptions, numerics, and the detailed formulation of the day-night thermal forcing. This contrasts with Jupiter and Saturn, for which circulation models can readily predict either eastward or westwardFootnote 5 equatorial jets depending on the detailed model assumptions (for reviews, see Vasavada et al. 2006; Del Genio et al. 2009; Showman et al. 2018a). This suggests that the mechanism that causes the superrotation on hot Jupiters is extremely robust.

Observed IR light curves of hot Jupiters are now sufficiently accurate that detailed comparison to GCM experiments is a fruitful exercise. Performing such comparisons self-consistently requires a model that represents the heating/cooling with an explicit radiative transfer scheme.Footnote 6 Most GCMs to date represent the angular dependence of radiation using the two-stream approximation. Approaches to treating the opacities come in several flavors. The simplest is a dual-band (double grey) scheme, with one opacity band in the IR and one in the visible (e.g., Heng et al. 2011a; Rauscher and Menou 2012a, 2013; Rauscher and Kempton 2014; Perna et al. 2012; Komacek et al. 2017; Tan and Komacek 2019; Mendonça et al. 2018; Flowers et al. 2019). The advantage of this approach is its simplicity—it is ideal for wide parameter explorations, understanding dynamical processes, and capturing the bulk radiative budget, which is sufficient for many applications.

On the other hand, the gaseous radiative transfer is inherently non-grey, with opacities that vary by orders of magnitude from wavelength to wavelength. For the price of increasing the model complexity, capturing the non-grey behavior affords two advantages. First, it allows a more realistic representation of the radiative heating/cooling, so that the overall thermal structure is accurately simulated. Second, IR spectra and lightcurves at different wavebands probe different pressures, where the temperature differ. As a result, IR spectra and lightcurves indicate that hot Jupiters are inherently non-grey bodies. Accurately capturing the wavelength-dependence of IR spectra and lightcurves requires non-grey radiative transfer. Several hot-Jupiter GCMs have been developed that treat the opacities using the well-known correlated-k method: the SPARC/MITgcm (Showman et al. 2009; Lewis et al. 2010; Kataria et al. 2013), and the Exeter group’s hot-Jupiter implementation of the UK Met Office GCM (Amundsen et al. 2014, 2016). In this approach, opacities are treated by dividing the spectra into typically 10–40 wavelength bins, and then statistically representing the opacity information from \(\sim 10^{4}\text{--}10^{5}\) wavelength points within each bin. This allows accuracy of typically 1% or better in heating rates (Showman et al. 2009; Kataria et al. 2013; Amundsen et al. 2014) while retaining much greater computational efficiency than a line-by-line radiative transfer calculation, which is computationally prohibitive in GCMs. Intermediate approaches include bin methods, which divide the spectrum into several dozen wavelength bins, but treat the opacity as essentially grey within each bin (Dobbs-Dixon et al. 2012; Dobbs-Dixon and Agol 2013). A key point is that none of these differing approaches are inherently better or worse; they should be viewed as distinct tools useful for different applications.

Detailed comparisons between IR lightcurve observations (obtained with Spitzer and Hubble) and GCM simulations including radiative transfer have been performed for a wide range of hot Jupiters, including HD 189733b (Showman et al. 2009; Knutson et al. 2012; Dobbs-Dixon and Agol 2013; Drummond et al. 2018a; Steinrueck et al. 2019; Flowers et al. 2019), HD 209458b (Zellem et al. 2014; Amundsen et al. 2016; Drummond et al. 2018b), WASP-43b (Kataria et al. 2015; Stevenson et al. 2017; Mendonça et al. 2018), WASP-18b (Arcangeli et al. 2019), WASP-19b (Wong et al. 2016), HAT-P-7b (Wong et al. 2016), WASP-103b (Kreidberg et al. 2018), WASP-121b (Parmentier et al. 2018), HAT-P-2b (Lewis et al. 2014), and HD 80606b (Lewis et al. 2017). Most of these studies compute the radiative transfer driving the GCM assuming a solar-composition, cloud-free hydrogen atmosphere that is in local chemical equilibrium, and perform straight-up comparisons of the observed lightcurves to a nominal model, without any special attempts to search for an exact match via model tuning. Nevertheless, a few of the models explore the influence of supersolar metallicities, nonsynchronous rotation, specified atmospheric drag (perhaps representing the effect of Lorentz forces), disequilibrium chemistry, or specified haze distributions on the lightcurves.

Generally speaking, the lightcurves predicted from these models capture the qualitative features of the observations—including large day-night flux differences and phase offsets wherein the fluxes peak before secondary eclipse, as expected due to the eastward offset of the dayside hotspot from the substellar point.Footnote 7 Given the lack of tuning, these aspects of qualitative agreement could be viewed as major successes suggesting the models are in approximately the correct regime. However, the GCM lightcurves also exhibit significant discrepancies from the observations—they tend to overpredict the hotspot offset, overpredict the nightside flux, underpredict the dayside flux, and therefore underpredict the day-night flux difference (i.e., the phase curve amplitude). Figure 5 shows examples for the benchmark planets HD 189733b, HD 209458b, WASP-43b and WASP-18b that highlight these successes and failures (see also Parmentier and Crossfield 2018, for a qualitative comparison). Although the degree of discrepancy varies from planet to planet, the persistence of these aspects of discrepancy across the hot-Jupiter sample suggests that the models may consistently be missing one or more crucial ingredients.

Fig. 5
figure 5

Observed thermal phase curves of benchmark hot Jupiters exoplanets compared to 3D global circulation models outputs. For HD209458b (top) we show the comparison of the Zellem et al. (2014) Spitzer observation at 4.5 μm with the SPARC/MITgcm output (left) and with the UK Met office Unified Model (right). For WASP-43b (middle row) we show the comparison between the SPARC/MITgcm and, from left to right, the Spitzer 4.5 μm phase curve, the Spitzer 3.6 μm phase curve, the HST/WFC3 phase curve and the HST/WFC3 dayside spectrum. For WASP-18b (bottom row left) we show the comparison between the HST/WFC3 spectral phase curve and spectrum at several phases and the SPARC/MITgcm output. Finally for HD189733b we show the four Spitzer phase curve and the related day and nightside spectrum and compare them to the SPARC/MITgcm outputs

As yet, there exists no consensus on what these missing factors may be, although several authors have explored the influence of various factors in an attempt to better match the observed lightcurves. Super-solar metallicities (e.g., 5 or 10 times solar) tend to increase the day-night flux differences and decrease the hotspot offset (Showman et al. 2009; Lewis et al. 2010; Kataria et al. 2015), which improves the agreement for certain planets, but generally do not solve the overall problem. For some planets, the addition of atmospheric drag improves the agreement (e.g., Kreidberg et al. 2018 for WASP-103b, and (Arcangeli et al. 2019) for WASP-18b); however, the Lorentz forces such drag parameterizes are only relevant at hot temperatures (\(\gtrsim2000~\text{K}\)), and so this effect cannot provide a solution for the cooler planets in Fig. 5. Chemical disequilibrium between CH4 and CO was suggested as a possible culprit for HD 189733b and HD 209458b (Knutson et al. 2012; Zellem et al. 2014), but recent GCM experiments rule out this explanation (Steinrueck et al. 2019; Drummond et al. 2018a, 2020). More likely is that the discrepancies result from clouds and/or hazes, although this possibility remains to be quantified in detail.

Measurements of the Doppler shift of the atmospheric winds on the planet’s terminator, as detected during transits, has been used to characterize the wind regime for several hot Jupiters. Using terminator-averaged measurements, Snellen et al. (2010) showed that the high-altitude winds on HD 209458b exhibit a preferential blueshift of \(2~\text{km}\,\text{s}^{-1}\), indicating that the atmospheric circulation comprises net flow toward Earth at high altitude (with the return flow presumably occurring at a deeper, unobserved level). Louden and Wheatley (2015, 2019) used transit ingress and egress measurements of HD 189733b and WASP-49b to discriminate between the leading and trailing limbs, and they showed that both planets exhibit equatorial superrotation (corresponding to a redshift on the leading limb, as the air flows from the planet’s nightside to dayside, and a blueshift on the trailing limb, carrying the return flow from the dayside back to the nightside). For HD 189733b, the observations indicate eastward velocities on the leading and trailing limbs of about \(3\text{--}4~\text{km}\,\text{s}^{-1}\) and \(4\text{--}6~\text{km}\,\text{s}^{-1}\), respectively. For WASP-49b, their study likewise identified atmospheric superrotation, and they were further able to determine that the circulation over both poles flows from day to night.

Motivated by these measurements, several GCMs investigations have explicitly determined the Doppler signatures that would result from the atmospheric circulation, and made predictions for various planets (Miller-Ricci Kempton and Rauscher 2012; Showman et al. 2013a; Rauscher and Kempton 2014; Flowers et al. 2019). These models show that at high altitude (pressures of \(\sim0.01\text{--}1~\text{mbar}\), the approximate levels sensed by the Doppler measurements), the circulation tends to comprise a day-night flow at high latitudes combined with a superrotation at lower latitudes, yielding a net signature where the windflow is redshifted at low latitudes on the leading terminator and blue shifted everywhere else along the terminator. These studies show that the relative contributions of the superrotation and the day-night flow to the limb signature depend on the stellar irradiation, planetary rotation rate, and the exact altitude sensed (which depends on the atmospheric composition), among other factors. The addition of sufficiently strong frictional drag or Lorentz forces can damp the jet, causing the day-to-night flow regime to predominate, although the regimes under which this can occur need further study. Generally, the GCM results agree reasonably well with the observed signatures, lending confidence that superrotation exists on hot Jupiters and that the GCMs are in approximately the correct regime.

2.3 Dynamical Mechanisms: Superrotation

We next turn to the dynamical mechanisms responsible for maintaining the flow, starting with the eastward equatorial jets prominent in hot Jupiter simulations (Fig. 4), as well as the mechanisms responsible for controlling the characteristic day-night temperature differences, wind speeds, and eastward offsets of the hot spots. Significant progress has been made on these topics in the past ten years.

The eastward equatorial jets emerging in circulation models of hot Jupiters are important for several reasons. Being perhaps the dominant aspect of the atmospheric flow, they are critical in explaining observations, because they help to cause the eastward shift of the dayside hot spot, influence the day-night temperature distribution, and cause differing conditions on the leading and trailing terminators. Such eastward equatorial jets are also dynamically interesting—they are an example of what is termed superrotation.Footnote 8 The equator is the region of the planet farthest from the rotation axis, and therefore, an eastward equatorial jet represents a local maximum of angular momentum per unit mass as a function of latitude and height—air at higher latitudes or deeper levels will tend to have smaller values of angular momentum per unit mass due to its closer proximity to the rotation axis. Maintaining such a superrotating jet against friction or other processes requires that the atmospheric circulation transport angular momentum up-gradient from regions where it is small (outside the jet) to where it is large (inside the jet). This is the opposite of diffusion, which causes down-gradient transport. Hide’s theorem (Hide 1969) states that such superrotating flows cannot result from axisymmetric circulations such as angular-momentum-conserving Hadley cells, but rather that the necessary up-gradient momentum transport can only be accomplished by waves and/or turbulence. Indeed, in many branches of fluid dynamics, up-gradient momentum transport is a fairly common outcome of turbulence and wave propagation (indeed, wave propagation is non-local in the sense that the wave can cause a flux of angular momentum that does not depend on the local gradients of the flow, as would occur in a diffusive problem). The puzzle in the present context is to understand the specific dynamical mechanism that causes the superrotation on hot Jupiters and to understand why it appears to be so robust.

The dynamical mechanism for equatorial superrotation on tidally locked planets was first investigated by Showman and Polvani (2010, 2011), who showed that the intense day-night heating pattern on synchronously rotating planets induces a standing pattern of large-scale atmospheric waves, and that these waves have a structure that naturally leads to a flux of angular momentum towards the equator, causing the superrotation.

To see how this works, it is useful to first consider the behavior of freely propagating tropical waves. Importantly, tropical waves tend to be confined to an equatorial waveguide whose half-width is equal to the equatorial Rossby deformation radius (Eq. (1)). Two classes of tropical wave are particularly relevant here—the Kelvin waves and the Rossby waves (see Holton and Hakim 2012; Andrews et al. 1987). The Kelvin waves are essentially gravity (buoyancy) waves in the east-west direction, with broad-scale pressure anomalies that are centered on, and symmetric about, the equator. In the north-south direction, Kelvin waves are geostrophically balanced—meaning the wave-induced pressure gradients are supported by Coriolis forces associated with the wave-induced wind anomalies. This balance implies (1) that isolated Kelvin waves exhibit very weak meridional wind relative to the zonal wind, and (2) that the zonal wind exhibits eastward phase in the regions of pressure maxima, and westward phase in the region of pressure minima. It turns out that trait (2) implies that only eastward-propagating waves are possible, and therefore Kelvin waves exhibit phase and group propagation to the east.

In contrast, Rossby waves are vortical waves whose restoring force relies on the planetary rotation, or more specifically, the fact that the Coriolis parameter is a function of latitude. This effect is referred to as the “\(\beta \) effect,” where \(\beta = df/dy\) is the gradient of the Coriolis parameter with northward distance \(y\). Because of the tendency of fluid parcels to conserve their potential vorticity following the flow,Footnote 9 the latitude variation of the planetary vorticity implies that fluid parcels deflected poleward tend to generate anticyclonic vorticity, whereas those deflected equatorward tend to generate cyclonic vorticity. This vorticity comprises both zonal and meridional motions, and it turns out that the meridional motions associated with this wave-generated vorticity tend to deflect the fluid parcel back to its original latitude. This chain of events thus acts as a restoring force that allows wave propagation. These properties imply that Rossby waves exhibit phase propagation to the west (Holton and Hakim 2012). When equatorially trapped, Rossby waves exhibit pressure patterns that are symmetric about the equator, with maxima that peak off the equator, and velocities skirting the pressure maxima in a vortical pattern. For low-order Rossby waves of long wavelength, this structure tend to fill the equatorial waveguide.

Time variable, stochastic forcing can trigger freely propagating Rossby and Kelvin waves, which in a system with weak damping are able to propagate zonally for long distances within the equatorial waveguide. For example, in Earth’s tropics, cumulus convection tends to trigger both of these wave classes, which can be observed because of their influence on the cloud patterns (Wheeler and Kiladis 1999; Kiladis et al. 2009). Models of stochastic forcing in the tropics likewise lead to these and other wave modes (e.g., Salby and Garcia 1987). A large heating pulse applied at the equator, with a size comparable to a deformation radius, for example, might naturally produce a Kelvin-wave packet that propagates to the east, and a Rossby-wave packet that propagates to the west.

Showman and Polvani (2010, 2011) demonstrated that the intense day-night thermal heating pattern on synchronously rotating exoplanets triggers a global-scale, standing wave pattern, with dynamics analogous to the tropical wave dynamics just described, and which explains the emergence of equatorial superrotation. These solutions, which built on the classic solutions of Matsuno (1966) and Gill (1980), adopted the 1.5-layer shallow-water equations, which represent the flow of an active, constant-density atmosphere (here representing the photosphere levels of a giant planet) that overlies a denser interior (representing the deep atmosphere and interior) that is assumed quiescent. Figure 6 depicts the linear, analytic, \(\beta\)-plane solutions for typical hot Jupiter conditions, and adopting radiative and drag time constants that are equal. In the solutions, the steady, day-night heating pattern (Fig. 6a) leads to standing eddy pattern known as a Matsuno-Gill pattern (Fig. 6b). At low latitudes, it comprises predominantly east-west flows that diverge from a point (marked by an X) lying east of the substellar point. This structure can be interpreted as a standing Kelvin wave; the eastward displacement of the thermal structure at low latitudes results from the fact that Kelvin waves propagate to the east. At higher latitudes it comprises vortical flow, with off-equatorial pressure maxima (anticyclones) in the northern and southern hemispheres of the dayside, and pressure minima (cyclones) in the northern and southern hemispheres of the nightside. This structure can be understood as an equatorially trapped Rossby wave; the westward displacement results from the fact that long-wavelength Rossby waves propagate to the west. Importantly, by distorting the thermal field eastward at low latitude (due to the Kelvin wave component) and westward at higher latitudes (due to the Rossby wave component), this latitudinally varying zonal-phase shift induces a global-scale chevron-shaped eddy pattern in which the pressure contours and eddy velocities are tilted northwest-southeast in the northern hemisphere and southwest-northeast in the southern hemisphere. These phase tilts are clearly visible in Fig. 6b. See Heng and Workman (2014) for a wider range of Matsuno-Gill-type solutions, and Pierrehumbert and Hammond (2019) for a tutorial on the Matsuno-Gill model.

Fig. 6
figure 6

Linear, analytical solutions of the shallow-water equations demonstrating how day-night heating can induce equatorial superrotation on tidally locked exoplanets, from Showman and Polvani (2011). Panel (a) shows the day-night heating pattern; specifically, the colors and contours represent the radiative equilibrium state as a function of eastward distance and northward distance (zero northward distance represents the equator and the center of the plot represents the substellar point; distance is plotted in units of the equatorial deformation radius, which for a hot Jupiter is a significant fraction of a planetary radius). The dayside is thick (hot) and the nightside is thin (cold). Panel (b) depicts the steady-state, linear solution that emerges in response to this heating pattern when no background flow is present. The solution can be interpreted as standing Kelvin and equatorially trapped Rossby waves. Notice that the predominant phase tilts of the eddy velocities are from northwest-southeast in the northern hemisphere and southwest-northeast in the southern hemisphere; this pattern causes transport of angular momentum from midlatitudes to the equator. Panel (c) shows the resulting acceleration of the zonal-mean zonal wind. The acceleration from meridional momentum convergence is in black, the acceleration from vertical momentum exchange is in blue, and the net acceleration is shown in red. In the net, eastward acceleration occurs at the equator, which will cause superrotation to emerge

This tilted pattern of eddy velocities causes an equatorward flux of momentum and is therefore key in driving the equatorial superrotation. Let \(u'\) and \(v'\) be the eddy zonal and meridional velocities, respectively, defined as the deviations of the zonal and meridional winds from their zonal-means. Examination of Fig. 6 shows that, in the northern midlatitudes, fluid parcels moving east (\(u' > 0\)) tend to be moving south (\(v' < 0\)), whereas those moving west (\(u' < 0\)) tend to be moving north (\(v' > 0\)). Thus, the velocities are correlated such that \(\overline{u'v'} < 0\), where the overbar denotes a zonal average. In the southern midlatitudes, the correlation is reversed—fluid parcels moving east tend to be moving north, whereas fluid parcels moving west tend to be moving south, such that \(u'v' > 0\). The quantity \(u'v'\) represents a meridional flux of eastward momentum per unit mass—a positive value implies northward transport of eastward momentum, whereas a negative value implies southward transport. Thus, the Matsuno-Gill pattern (Fig. 6) naturally causes a convergence of momentum onto the equator, which leads to an eastward acceleration that drives equatorial superrotation. In this theory, the meridional width of these Matsuno-Gill solutions—and the equatorial jet that they drive—is comparable to the equatorial deformation radius.

The solutions also show that, at the equator, the standing waves cause a net downward transport of zonal momentum, which plays a critical role in the momentum balance. Because of the eastward phase shift of the Kelvin-wave structure, the equatorial eddy velocity on the nightside is preferentially eastward while that on the dayside is preferentially westward (see Fig. 6b). Nightside cooling causes downward transport of air, which takes its eastward eddy momentum with it. Dayside heating brings up air from the deeper atmosphere with only weak zonal velocity. Thus, there is a net downward transport of eastward momentum out of the atmosphere into the interior. This loss of momentum from the atmosphere causes a westward acceleration at the equator that partially—but not completely—cancels out the eastward acceleration caused by the convergence of momentum onto the equator. (In Fig. 6c, the acceleration due to this vertical exchange is shown in blue, that due to meridional momentum convergence is shown in black, and the net acceleration is shown in red.)

The solution shown in Fig. 6 adopts equal radiative and drag timescales, but the linear behavior differs when these timescales are unequal. In particular, Showman and Polvani (2011) and Heng and Workman (2014) explored a wide range of differing \(\tau _{\mathrm{rad}}\) and \(\tau _{\mathrm{drag}}\); they showed analytically that, as the drag becomes weak (i.e., the drag timescale becomes long), the equatorial waveguide corresponding to these forced-damped solutions becomes confined closer and closer to the equator. Their solutions show that, in the linear limit of \(\tau _{\mathrm{drag}} \rightarrow \infty \), the meridional width of the region with prograde phase tilts shrinks to zero,Footnote 10 and a pattern of reverse phase tilts instead emerges—northeast-southwest in the northern hemisphere and southeast-northwest in the southern hemisphere. Such a pattern would not be expected to converge eddy momentum onto the equator or drive strong superrotation. The presence of these reverse phase tilts are well-understood (Showman and Polvani 2011, Appendix C), but this leads to a puzzle: how can superrotation develop in a 3D hot Jupiter GCM with no explicit frictional drag, when the meridional width of the Matsuno-Gill pattern (and its prograde phase tilts) should shrink to zero?

The answer is provided by nonlinearity. At the high forcing amplitudes relevant to typical hot Jupiters, the dynamics becomes nonlinear, which modifies the details of the wave-mean-flow interactions, both due to the nonlinearity of the eddies themselves, and because a strong zonal flow develops when the forcing amplitude is large, and the zonal flow modifies the eddy structures. Showman and Polvani (2011) showed numerically that when nonlinearity is progressively introduced back into the dynamics, the meridional width of the Matsuno-Gill-type pattern broadens accordingly, even in the limit of weak drag, and at very large amplitudes exhibits a meridional width comparable to the equatorial deformation radius. This effect of nonlinearity on the eddies can be understood by appreciating that nonlinear momentum advection can qualitatively play the role of drag in the multi-way force balance that controls eddy behavior, allowing prograde phase tilts to persist to high latitudes, and preventing the Matsuno-Gill-type behavior from collapsing to the equator (Showman and Polvani 2011; for a visual explanation of how this effect can lead to prograde phase tilts, see Showman et al. 2013a). Mathematically, this effect can be appreciated from the momentum equation governing the 1.5-layer shallow-water system (Showman and Polvani 2011, Eq. (12)):

$$ \frac{d\mathbf{v}}{dt}+g\nabla h+f\mathbf{k}\times \mathbf{v}=- \mathbf{v}\left [\frac{1}{\tau _{\mathrm{drag}}} + \frac{1}{\tau _{\mathrm{rad}}}\left (\frac{h_{\mathrm{eq}}-h}{h}\right ) \mathcal{H}(h_{\mathrm{eq}}-h)\right ] $$
(7)

where \(\mathbf{v}\) is the horizontal velocity vector, \(d/dt\) is the material derivative, \(h\) is layer thickness, \(\mathbf{k}\) is the local vertical unit vector, \(h_{\mathrm{eq}}\) is the local radiative equilibrium thickness, and ℋ is the Heaviside step function, which equals one when its argument is positive and zero otherwise. The entire quantity in square brackets plays the role of one over an “effective drag” time constant. At low amplitude, the quantity \((h_{\mathrm{eq}} - h)/h \ll 1\), and thus the second term in square brackets is not generally important. At high amplitude, however, deviations from radiative equilibrium tend to be substantial,Footnote 11 implying that \((h_{\mathrm{eq}} - h)/h \sim 1\). This implies that, on the dayside, the second term in square brackets tends to have a magnitude of order \(1/\tau_{\mathrm{rad }}\). Thus, under typical hot-Jupiter conditions, the “effective” drag time constant caused by the nonlinearity has comparable magnitude to the radiative time constant—even when the actual drag time constant \(\tau {\mathrm{drag}}\) is infinite. The fact that the highly nonlinear, zero-drag regime tends to naturally exhibit radiative and effective drag timescales that are comparable to each other can explain why these solutions have a Matsuno-Gill-like pattern that extends to high latitudes, just like linear solutions with equal radiative and drag timescales (Fig. 6). This behavior also explains more generally why the standard definition of equatorial deformation radius (Eq. (1)) provides a reasonable fit to the actual jet width in most 3D hot-Jupiter GCMs, despite the fact that these GCMs generally do not include explicit frictional drag in the region where the superrotation forms. Debras et al. (2020) further confirmed the above picture by performing full 3D linear wave calculations using the numerical package developed by Debras et al. (2019).

In sum, high-amplitude, fully nonlinear shallow-water solutions on the sphere show that essentially the same mechanism holds even in the nonlinear regime: although the details of the wave-mean-flow interactions are sensitive to nonlinearity, the qualitative mechanism still occurs even at high amplitude. As mentioned above, the eddy nonlinearity is critical for preventing the Matsuno-Gill solutions from collapsing to the equator in the weak-drag limit (Showman and Polvani 2011). Once the jet builds up, these planetary-scale waves can exhibit finite meridional width even in the linear, weak-drag limit (Hammond and Pierrehumbert 2018; Debras et al. 2020), presumably because they are then propagating relative to the moving airflow (even though they are standing in the reference frame of the planet). As the jet builds up, the concomitant changes in the eddy structure lead to changes in the amplitudes of the horizontal and vertical eddy momentum convergences. Steady state is achieved when they balance (along with drag, if any).

In a major study, Tsai et al. (2014) showed that a deeply stratified, 3D atmosphere exhibits global-scale, Matsuno-Gill-type wave solutions that are very similar to those in shallow water. They adopted radiative and frictional damping times that are equal and invariant with depth. The day-night heating was assumed to occur over a broad vertical layer peaking around 0.1 bar. Tsai et al.’s solutions in the absence of a background flow, shown in the left column of Fig. 7, yield a Matsuno-Gill pattern very similar to the shallow-water solutions identified by Showman & Polvani, comprising superposed equatorial Kelvin and Rossby waves whose differential zonal propagation produces eddies tilting predominantly northwest-southeast (southwest-northeast) in the northern (southern) hemisphere. Just as in the shallow-water system, the eddy structure induces a meridional momentum convergence onto the equator in the forcing layer, as well as a downward momentum transport at the equator from the forcing levels to deeper levels. Within the forcing layer, the net accelerations in the meridional plane are eastward at the equator and westward at higher latitudes (Fig. 7d), qualitatively similar to the pattern of accelerations predicted by the shallow-water solutions (Fig. 6c). The striking similarity between the shallow-water and 3D solutions is not surprising, because, as is standard in tidal theory and many other wave problems (Chapman 1970), this 3D problem is mathematically separable and determined by solving a shallow-water equation for the horizontal structure, coupled to a 1D vertical-structure equation for the vertical eigenfunctions.

Fig. 7
figure 7

Linear, analytic, steady \(\beta \)-plane solutions of the Matsuno-Gill problem in 3D primitive equations, assuming equal radiative and frictional time constants that are constant with depth, from Tsai et al. (2014). Left, middle, and right columns adopt a uniform background wind of \(U=0\), 500, and \(1000~\text{m}\,\text{s}^{-1}\), respectively. Top row shows eddy geopotential and eddy wind structure; antistellar point lies in the middle of each plot. Bottom row shows the total wave-induced acceleration; exact values are arbitrary but for plausible forcing amplitudes range from \(-10^{-6}~\text{m}\,\text{s}^{-2}\) (westward) in dark blue to \(3\times 10^{-6}~\text{m}\,\text{s}^{-2}\) (eastward) in red

In Tsai et al. (2014)’s solutions, below the \(\sim 1\text{-bar}\) level, the waves propagate downward away from the forcing layer. This propagation causes the Rossby gyres to shift westward at greater depth, while the Kelvin wave component shifts eastward with depth. The overall longitudinal/vertical tilts of the solution depend on the relative amplitudes of the Kelvin and Rossby components, which vary with the adopted damping time constant. In particular, the Kelvin wave component dominates when the damping timescales are short, while the Rossby wave component dominates when the damping timescales are long. Despite these complications, the visual appearance of their solutions is that the vertical tilts appear to be modest under damping conditions relevant to hot Jupiters.

Still, analysis of fully nonlinear GCMs shows that the build-up of the equatorial jet strongly modifies the standing waves and disrupts the coherency of the prograde phase tilts, leading to a wave structure in the fully equilibrated state that differs substantially from the classical Matsuno-Gill pattern (Showman and Polvani 2011; Tsai et al. 2014; Showman et al. 2015; Mayne et al. 2017; Mendonça 2020). Because of the disruption of the phase tilts, the net meridional momentum transports—while still equatorward|are weaker than would be predicted from the classic Matsuno-Gill model at the same wave amplitude. Quantitatively understanding how the equatorial jet equilibrates therefore requires an understanding of how the jet modifies the eddies. Several authors have addressed this question.

The addition of a uniform background zonal wind to the Matsuno-Gill model has long been known to modify the standing wave modes (Phlips and Gill 1987) and the momentum fluxes they induce (Arnold et al. 2012). To investigate this effect in the context of hot Jupiters, Tsai et al. (2014) included a uniform background zonal wind \(U\) in their analytic solutions (Fig. 7, middle and right columns). The background zonal flow influences the amplitudes and relative phase offset between the Kelvin and Rossby wave components, with crucial implications for the wave-induced accelerations. Because the group velocity for long Rossby waves is about one-third that of the Kelvin wave, the phase offset of the Rossby-wave component is much more easily influenced by the background winds. In the solution in Fig. 7 background winds ranging from zero to \(1~\text{km}\,\text{s}^{-1}\) hardly change the offset of the Kelvin-wave component, but they cause the Rossby wave gyres to shift from \(\sim 20^{\circ }\) west of the substellar longitude to \(50^{\circ }\) east of it. (Compare Fig. 7a and c; see Tsai et al. 2014, Fig. 8, for a decomposition of the Kelvin and Rossby-wave components of the solutions.) At zonal winds of \(\sim1~\text{km}\,\text{s}^{-1}\), the pressure maxima of the Rossby and Kelvin components are nearly in phase, and as can be seen from Fig. 7, this destroys the strong polewared-westward to equatorward-eastward phase tilts. After a transient increase in the strength of the acceleration at modest \(U\) (Fig. 7e), the zonal acceleration then plummets as the zonal flow increases (Fig. 7f). Still, the net acceleration remains eastward throughout the entire range of \(U\) shown in Fig. 7, suggesting that the jet will continually accelerate to speeds exceeding \(1~\text{km}\,\text{s}^{-1}\). Tsai et al. (2014) showed that, for the particular model parameters of Fig. 7, the net, wave-induced zonal acceleration at the equator reaches zero only when the background wind reaches about \(3~\text{km}\,\text{s}^{-1}\). This intriguing result could help to explain why the equilibrated zonal jet speeds in GCMs are typically \(3\text{--}4~\text{km}\,\text{s}^{-1}\) (Fig. 4), although the meridional and vertical shear of the zonal wind are also likely important.

Hammond and Pierrehumbert (2018) explored the effect of meridional shear of the equatorial jet on the wave structure and momentum fluxes in a linear shallow-water model. They assumed a Gaussian-shaped equatorial jet whose zonal-mean wind is \(\overline{U}=U_{0}e^{-y^{2}/L^{2}_{\mathrm{eq}}}\), qualitatively similar to those in GCM solutions. Their results indicate that, even in the presence of the sheared equatorial jet, the eddies still transport momentum from midlatitudes onto the equator, leading to an eastward acceleration at the equator that can help maintain the superrotation. Furthermore, they showed that the inclusion of the meridionally shearing jet and its associated geopotential anomaly triggers more forced wave modes, which results in a better match between the overall surface height from the linear shallow-water model and the GCM temperature field than those without a jet.

Several authors have further investigated the dynamics of the equatorial superrotation in full 3D GCM experiments. As mentioned above, the equilibrated, standing eddy structure shows significant distortion from the classic Matsuno-Gill pattern. Despite this complication, diagnostics from these GCMs still generally support a picture wherein the day-night forcing leads to a large-scale, standing wave pattern that causes a meridional convergence of eddy angular momentum onto the equator, driving the superrotation, and a downward eddy momentum transport at the equator from above the photosphere to deeper levels (Showman and Polvani 2011; Tsai et al. 2014; Showman et al. 2015; Mayne et al. 2017; Mendonça 2020; Debras et al. 2020). This downward eddy momentum transport helps cause the equatorial jet to slowly penetrate more deeply over time at pressures of a few bars and greater. At mid-to-high latitudes, all the studies agree that the steady state zonal-momentum balance involves not only the eddy flux convergences (both meridional and vertical), but also the Coriolis accelerations and momentum advection associated with the mean-meridional circulation. At the equator, the zonal-momentum balance at and above the photosphere tends to comprise a balance primarily between the (eastward) acceleration due to meridional momentum transport and the (westward) acceleration due to the vertical momentum transport (e.g., Showman and Polvani 2011). Nevertheless, several authors have highlighted the role of vertical momentum transport due to the zonal-mean circulation as well (Mayne et al. 2017; Mendonça 2020). Interestingly, although the primary day-night forcing would appear to comprise a zonal wavenumber-one (diurnal) mode, Mendonça (2020) argues that the wavenumber-two (semidiurnal) mode may be at least as important. The details of all these issues are likely to be sensitive to planetary rotation rate, strength and implementation of radiative forcing, and other modeling aspects, so it would be worth exploring these issues over a wider range of conditions.

To summarize, the mechanism predicts that the equatorial jet has a meridional half-width approximately equal to the Rossby deformation radius. For conditions typical of hot Jupiters, this length scale is comparable to a planetary radius, explaining the broad structure of the equatorial jet in hot-Jupiter circulation models (Fig. 4). The process that drives the jet is, in its essence, a direct wave-mean-flow interaction between the eddies and the planetary-scale mean flow. Turbulent cascades or other eddy-eddy interactions are possible, and would affect the details, but are not essential to the basic mechanism.

2.4 Day-Night Temperature Differences and Recirculation Efficiency

Synchronously rotating exoplanets exhibit permanent daysides and nightsides and are therefore subject to a fascinating climate regime of day-night thermal forcing that has no analogy in the solar system. Lacking any direct irradiation, the only possible means for the nightside to exhibit a temperature exceeding tens of K is for heat to be transported from the dayside to the nightside by dynamical motions in the atmosphere or interior. This raises the intriguing question of what processes control this day-night heat transport, what the expected day-night temperature difference is, and how it should vary with the stellar irradiation, planetary rotation rate, atmospheric composition, and other parameters. This problem is analogous to the important question of what controls the meridional temperature gradients and equator-to-pole temperature differences on planets like Earth, but in the case of hot Jupiters, the heat transport occurs not only meridionally but zonally as well. Because the \(\beta \) effect inhibits heat transport meridionally, whereas no such constraint exists in the zonal direction, the dynamics are therefore likely to be quite different.

This problem is intriguing not only because it touches on basic issues in geophysical fluid dynamics (GFD) but because it is directly amenable to observational characterization: on a synchronously rotating planet, energy balance implies that the energy transported from the dayside to the nightside by the atmosphere/interior circulation is radiated to space on the nightside. Therefore, measurements of the total IR flux radiated to space on the nightside—which can be directly estimated from IR lightcurves—provides an observed measure of the day-night dynamical heat transport.

Many authors have posited that the problem of whether the day-night temperature difference is large or small can be cast as a comparison between two timescales, the radiative timescale (Eq. (3)) and an advection timescale for air to travel across a hemisphere from day to night:

$$ \tau _{\mathrm{adv}} \sim \frac{a}{U}, $$
(8)

where \(a\) is the planetary radius and \(U\) is a characteristic horizontal wind speed. Showman and Guillot (2002) first suggested that hot Jupiters will exhibit large day-night differences when \(\tau _{\mathrm{rad }}\ll \tau_{\mathrm{adv}}\) and small day-night temperature differences when \(\tau _{\mathrm{rad }}\gg \tau _{\mathrm{adv}}\). Many authors have suggested that this timescale comparison helps to explain the dependence of the day-night temperature difference on pressure, opacity, and incident stellar irradiation (Cooper and Showman 2005; Showman et al. 2008; Fortney et al. 2008; Rauscher and Menou 2010; Heng et al. 2011b; Cowan and Agol 2011b,a; Perna et al. 2012).

The radiative timescale scales as \(T^{-3}\) (Eq. (3)), implying that it becomes very short for ultra hot Jupiters. In contrast, order-of-magnitude expressions for the typical horizontal velocities (Eqs. (5)–(6)) suggest that the horizontal velocity should also be greater on strongly irradiated hot Jupiters, but that the dependence is weaker. Thus, one might expect that for hot Jupiters whose photospheres exceed some critical temperature—perhaps of order 2000 K—the atmosphere is close to radiative equilibrium and the day-night temperature difference is large. For hot Jupiters cooler than this value, the day-night temperature difference should become smaller.

Perna et al. (2012) were the first authors to systematically explore the influence of stellar irradiation on day-night temperature structure in GCM simulations. Their GCM adopted a dual-band radiative transfer approach with one broad opacity band in the visible and one in the IR, neglecting scattering, with opacities chosen to yield a range of plausible temperature profiles. In a series of numerical experiments, they varied the stellar irradiation over nearly three orders of magnitude—corresponding to planetary effective temperatures of about 600 to 2500 K—and examined how the circulation responds. Their simulations show that, as qualitatively expected, the day-night temperature differences are large for effective temperatures \(\gtrsim 2000~\mbox{K}\) but becomes small for \(T \lesssim 1500~\mbox{K}\) (Fig. 8). They showed that the transition is consistent with the expectation based on the timescale argument presented above: for planets with effective temperatures less than \(\sim1500~\text{K}\), the advective and radiative timescales (as post-processed from the simulation results) are similar, but at effective temperatures exceeding 1500 K, the radiative timescale becomes shorter than the advective timescale, and the day-night flux differences become large (Fig. 8).

Fig. 8
figure 8

Fractional day-night IR flux differences in GCM simulations of synchronously rotating hot Jupiters as a function of their equilibrium temperature for models with different assumptions. Perna et al. (2012) uses the FMS model with semi-grey radiative transfer. Komacek et al. (2017) use the MITgcm model with semi-grey radiative transfer but slightly different choices of opacities and rotation period than Perna et al. (2012). Parmentier et al. (2020) uses the SPARC/MITgcm with non-grey, temperature dependent opacities (assuming chemical equilibrium without TiO/VO) and rotation periods similar to Perna et al. (2012). Despite their differences, the three models provide similar results when used in similar conditions, showing that the decrease of heat redistribution efficiency with increasing temperature is robust. However, when either drag (orange triangles from Komacek et al. 2017) or the radiative feedback of nightside clouds (diamonds from Parmentier et al. 2020) are included in the models, the heat redistribution becomes much smaller. The difference between physical ingredients is much larger than the difference between modelling frameworks

There are several issues with the timescale comparison between \(\tau_{\mathrm{rad }}\) and \(\tau _{\mathrm{adv }}\), however (Perez-Becker and Showman 2013; Komacek and Showman 2016). First, it was not derived in any self-consistent manner from the governing dynamical equations, but rather represents an ad hoc, albeit attractive, hypothesis. Second, it is not predictive—the advection timescale can only be evaluated once the wind speeds are known, yet the wind speeds themselves depend on the temperature differences one is trying to understand. Third, the criterion obscures any role for other key timescales, including rotational timescale, frictional timescale (if any), and timescale for horizontal and vertical wave propagation. These timescales surely affect the circulation, including the day-night temperature differences, and therefore one wound expect a timescale criterion for day-night temperature differences to depend on them.

Perez-Becker and Showman (2013) and Komacek and Showman (2016) performed idealized numerical experiments and constructed analytic scaling theories for the day-night temperature differences and characteristic wind speeds, with the aim of obtaining a self-consistent understanding of both the day-night temperature differences and characteristic wind speeds (hence advection timescale). The former study adopted the 1.5-layer shallow water model, while the latter adopted the 3D primitive equations, but the setup was otherwise quite similar across both studies, involving an idealized day-night forcing where the thermal structure was relaxed towards a prescribed radiative equilibrium (hot on the dayside and cold on the nightside) over a prescribed radiative timescale \(\tau _{\mathrm{rad }}\). In agreement with the results of Perna et al. (2012), the simulations show that the thermal structure exhibits essentially no day-night variation when \(\tau _{\mathrm{rad }}\) is long, but that the day-night contrast becomes large when \(\tau _{\mathrm{rad }}\) is short (Fig. 9). The transition occurs near \(\tau _{\mathrm{rad}} \sim 10^{5}~\text{s}\). Frictional drag was also included in some models to parameterize the possible effects of Lorentz forces, and these authors found that sufficiently strong drag could also lead to large day-night temperature differences.

Fig. 9
figure 9

Dependence of circulation in synchronously rotating hot Jupiter models over a wide range of radiative time constant. The radiative heating is implemented with an idealized Newtonian cooling scheme that relaxes the circulation toward a specified radiative equilibrium structure. Top row presents shallow-water simulations from Showman et al. (2013a) and Perez-Becker and Showman (2013); from left to right panels, radiative time constant is 100, 10, 1, 0.1, and 0.01 Earth days, respectively. Rotation period is 2.2 days in all models. Plot shows thickness (colorscale) and flow velocity (arrows). Bottom row: 3D primitive equation simulations from Komacek and Showman (2016) showing temperature (colorscale) and flow velocity (arrows) structure at 80 mbar; from left to right panels, radiative time constant is \(10^{7}\), \(10^{6}\), \(10^{5}\), \(10^{4}\), and \(10^{3}~\text{s}\), respectively. Rotation period is 3.5 days in all models. Both sets of models exhibit a transition from small day-night thermal contrast at \(\tau_{\mathrm{rad }} \gtrsim 10^{5}~\text{s}\) to large day-night contrast at \(\tau _{\mathrm{rad }} \lesssim 10^{5}~\text{s}\)

The analytic theory in Perez-Becker and Showman (2013) and Komacek and Showman (2016) is constructed by balancing dominant terms in the dynamical equations. Ginzburg and Sari (2016) also present a related analysis. In the momentum equation, the day-night pressure-gradient force, which emerges from the day-night temperature difference and drives the circulation, is balanced by the greater of the Coriolis, frictional drag, horizontal momentum advection, or vertical momentum advection forces, respectively; this leads to different expressions for the day-night contrast in each regime, along with criteria for the transition between them. For example, when the dominant force balance is between the day-night pressure-gradient force and the Coriolis force, which holds when drag is weak and the rotation is sufficiently fast, the theory predicts that the fractional day-night contrast is

$$ A=\left (1+\frac{\tau _{\mathrm{rad}}}{f\tau ^{2}_{\mathrm{wave}}}\right )^{-1}, $$
(9)

where here \(A\) is the fractional day-night thermal contrast (i.e., the fractional day-night thickness contrast in shallow water and fractional day-night temperature difference in 3D). \(f=2\Omega \sin \phi \) is the Coriolis parameter, \(\Omega \) the planetary rotation rate and \(\phi \) the latitude. The quantity \(\tau _{\mathrm{wave }}=a/NH\) is the characteristic timescale for gravity (or Kelvin) waves of long vertical wavelength to propagate over a zonal distance equal to a planetary radius. The typical horizontal wave propagation speed can be estimated by \(NH\), which is about \(1~\text{km}\,\text{s}^{-1}\) for a hot Jupiter (see Sect. 2.1), which leads to a timescale \(\tau _{\mathrm{wave }}\sim 10^{5}~\text{s}\), where we have adopted a Jupiter radius \(a = 7 \times 10^{7}~\text{m}\). In the limit of large \(\tau _{\mathrm{rad }}\), this equation predicts \(A = 0\), corresponding to no variation from day-to-night; in the limit of small \(\tau_{\mathrm{rad}}\), it predicts \(A=1\), corresponding to a temperature structure in radiative equilibrium, exhibiting large day-night temperature differences. The transition between these limits occurs when

$$ \tau _{\mathrm{rad}} \sim f\tau ^{2}_{\mathrm{wave}}, $$
(10)

where \(f\) is the Coriolis parameter. Evaluating the expression for a typical hot Jupiter value \(f=3\times 10^{-5}~\text{s}^{-1}\) yields a critical radiative time constant \(\tau _{\mathrm{rad }}\sim 3\times 10^{5}~\text{s}\). Thus, under typical hot Jupiter conditions, this theory predicts that day-to-night temperature differences will be small when \(\tau _{\mathrm{rad}} \lesssim 3\times 10^{5}~\text{s}\) and large when \(\tau _{\mathrm{rad }} \gtrsim 3\times 10^{5}~\text{s}\). It can be seen that this explains the transition occurring in Fig. 9. More generally, the theory yields self-consistent predictions of the day-night thermal contrast and wind speeds over the full parameter space of \(\tau _{\mathrm{rad }}\), \(\tau _{\mathrm{rad}}\), forcing amplitude (i.e., the day-night difference in radiative-equilibrium temperature) and other parameters, with no free parameters and no tuning. Perez-Becker and Showman (2013) and Komacek and Showman (2016) showed that the theory matches the simulations reasonably well over this wide range. In contrast, the comparison between radiative and horizontal advection timescales does not, particularly when the forcing amplitude is weak.

Perez-Becker and Showman (2013) and Komacek and Showman (2016) showed that the dynamical mechanism for regulating the day-night temperature differences can be interpreted in terms of a wave-adjustment mechanism: the day-night heating contrast triggers planetary-scale waves—including the Kelvin and Rossby wave discussed previously—which propagate longitudinally around the planet. The horizontal convergence and divergence associated with these waves changes the vertical thickness of fluid columns, inducing vertical motion that pushes isentropes up or down. One can make an analogy to throwing a rock into a pond: the perturbation to the water surface caused by throwing the rock into the pond leads to the propagation of gravity waves whose horizontal convergence or divergence change the thickness; once the waves radiate away, the end result is a flat surface. In the atmosphere, a similar process acts to flatten surfaces of constant entropy (for a review, see Showman et al. 2013b). If the damping is weak (i.e. the damping timescales are long), the waves can easily propagate from dayside to nightside and this adjustment process is efficient, leading to small day-night temperature difference. If the damping is strong (i.e., the damping timescales are short), however, then the waves become damped before they can propagate from day to night. This suppresses the wave regulation mechanism, allowing the temperature difference to approach radiative equilibrium.

What is the relationship of this criterion to the comparison between radiative and advective timescales? Perez-Becker and Showman (2013) and Komacek and Showman (2016) showed that the timescale criterion for transition between large and small fractional day-night temperature differences can be expressed as a comparison between the radiative timescale and an appropriately defined vertical advection timescale \(\tau _{\mathrm{vert}}\). This comparison emerges naturally from the dynamical equations because, on global scales, the vertical entropy advection term generally dominates over the horizontal entropy advection term in the thermodynamic energy equation (see Komacek and Showman 2016 for a discussion). This is essentially the “weak temperature gradient” (WTG) balance well-known from Earth tropical meteorology (for reviews, see Showman et al. 2013b; Pierrehumbert and Hammond 2019). The timescale comparison between \(\tau _{\mathrm{rad }}\) and \(\tau_{\mathrm{vert}}\) subsumes the classic comparison between \(\tau_{\mathrm{rad }}\) and \(\tau _{\mathrm{adv }}\): they are essentially the same in the highly forced limit where the day-night temperature difference is order unity, but they differ greatly in the weak forcing limit, and in that limit, the former comparison explains the simulation results far better than the latter one.

Zhang and Showman (2017) combined these regimes into a single compact expression for the temperature difference, valid across the entire parameter space:

$$ \frac{\Delta T}{\Delta T_{\mathrm{eq}}} = 1- \frac{2}{\alpha +\sqrt{\alpha^{2}+4\gamma ^{2}}} $$
(11)

where the non-dimensional, pressure-dependent constants \(\alpha \) and \(\gamma \) are defined as

$$ \alpha =1+ \frac{\left (\Omega +\frac{1}{\tau _{\mathrm{drag}}}\right ) \tau ^{2}_{\mathrm{wave}}}{\tau _{\mathrm{rad}}\delta \ln p} $$
(12)

and

$$ \gamma = \frac{\tau ^{2}_{\mathrm{wave}}}{\tau _{\mathrm{rad}}\tau _{\mathrm{adv},\mathrm{eq}}\delta \ln p}. $$
(13)

In these expressions, \(\tau _{\mathrm{adv},\mathrm{eq}} = a/U_{\mathrm{eq}}\) is the timescale for advection by a reference “equilibrium cyclostrophic wind”, defined as \(U_{\mathrm{eq}} = \sqrt{\Delta T_{\mathrm{eq}} \delta \ln p/2}\). The quantity \(\delta \ln p\) is the difference in log-pressure between some deep pressure where the flow is assumed to decay to zero (assumed to be 10 bars by Komacek and Showman 2016) and some lower pressure of interest, for example at the IR photosphere.

Atmospheres of ultra-hot Jupiters (here loosely defined as having equilibrium temperature \(\gtrsim 2000~\mbox{K}\)) are sufficiently irradiated such that all molecules, including molecular hydrogen and water are partially dissociated in their dayside atmospheres (see Fig. 10 and Parmentier et al. 2018). Heat is absorbed on the dayside to break the strong H2 bond, and then the atomic hydrogen is transported by atmospheric flows towards cooler regions of the atmosphere, where hydrogen atoms recombine to molecules and release heat. The energy per unit mass associated with the molecule-atom transition of hydrogen (\(2.18 \times 10^{8}~\text{J}\,\text{kg}^{-1}\)) is two orders of magnitude higher than that associated with water phase change, and together with the fact that hydrogen is the main atmospheric constituent, suggests that heat transport by hydrogen dissociation and recombination can be enormous in ultra-hot atmospheres. Bell and Cowan (2018) first studied this effect using an semi-analytic energy transport model, and they found that the heating from hydrogen recombination could significantly reduce the day-night phase-curve amplitude and increase phase curve offsets. Building upon the theory of Komacek and Showman (2016), Komacek and Tan (2018) analytically investigated this effect on the day-night heat transport of ultra-hot Jupiters. They found that the day-night temperature difference increases with increasing temperature until reaching an equilibrium temperature of about 2300 K, after which the day-night temperature difference decreases with increasing temperature. This finding can qualitatively explain the relatively small phase curve amplitudes observed for several ultra-hot Jupiters (Fig. 11). Next, Tan and Komacek (2019) implemented the effects of heating and cooling, change of mean molecular weight and gas thermodynamic properties caused by hydrogen recombination and dissociation into an idealized GCM similar to that in Komacek et al. (2017), and they systematically investigated the influence on the atmospheric circulation of ultra-hot Jupiters. They likewise confirmed that the fractional day-night temperature difference as a function of equilibrium temperature peaks at \(\sim2300~\text{K}\), in good agreement with the prediction of Komacek and Tan (2018). Additionally, wind speeds and direction of the equatorial flows at relatively low pressures are strongly affected by including the hydrogen dissociation and recombination, and this influence is more prominent with higher equilibrium temperature.

Fig. 10
figure 10

Percentage of water that is dissociated at the infrared photosphere of irradiated planets as a function of their dayside temperatures and gravity. Planets on the right hand side of the \(20\%\) of water dissociation iso-contour can be considered as ultra hot Jupiters. The efficiency of thermal dissociation is pressure dependent, hence why planets with higher gravity, tend to be less affected by it. A stellar isochrone for an age of 4 Gyr from Baraffe et al. (2008) is shown as a dashed black line

Fig. 11
figure 11

Predictions for day-night temperature differences from analytic theories and comparison to measurements from lightcurves of various hot Jupiters, from Komacek and Tan (2018). Dotted lines show the predictions from Komacek and Showman (2016), and solid lines show the predictions from Komacek and Tan (2018) accounting for the latent heat caused by dissociation and recombination of hydrogen at high temperature. Orange curve assume there is no explicit frictional drag in the observable atmosphere, while the black curves adopt a spatially constant drag coefficient with a timescale of \(10^{4}~\text{s}\)

2.5 Clouds, Hazes, and Chemistry on Hot Jupiters

There now exist several lines of evidence for clouds and hazes in the atmospheres of hot Jupiters. First, many hot Jupiters exhibit transit (transmission) spectra that are considerably flatter than predicted by cloud-free models; by causing broadband scattering, high-altitude clouds/hazes can flatten spectral features in the observed way. Moreover, in the visible, these transmission spectra commonly exhibit a Rayleigh-like slope attributed to scattering by small aerosol particles (e.g., Sing et al. 2011, 2016; Pont et al. 2013; Stevenson 2016; Iyer et al. 2016; Heng 2016; Barstow et al. 2017; Wakeford et al. 2017). Second, measurements of secondary eclipse at visible wavelengths imply that several hot Jupiters exhibit sufficiently high albedo to require clouds (e.g., Evans et al. 2013; Heng and Demory 2013; Esteves et al. 2013, 2015; Angerhausen et al. 2015). Energy-balance arguments for planets with full-orbit IR lightcurves also imply significant Bond albedos for several hot Jupiters (Schwartz and Cowan 2015). Most spectacularly, Kepler visible lightcurves of several hot Jupiters peak after secondary eclipse, indicating that the dayside bright spot in visible wavelengths is west, rather than east, of the substellar point (Demory et al. 2013; Angerhausen et al. 2015; Esteves et al. 2015; Shporer and Hu 2015). This is interpreted as an inhomogeneous cloud distribution, with a cloud-dominated region west of the substellar point.

Atmospheric circulation modulates the 3D distribution of condensate clouds, but in turn the cloud distribution influences the radiative heating/cooling and therefore the circulation. We consider these issues in turn, followed by brief discussions of cloud-induced variability, hazes, and coupling between chemistry and dynamics.

2.5.1 How Does the Atmospheric Circulation Affect the Cloud Distribution?

Two main processes can shape the cloud distribution in hot Jupiter atmospheres. First, the temperature affects the cloud distribution, leading to a day/night contrast in the cloud abundance. Second, the meridional atmospheric circulation leads to an equator to pole cloud variation that superimposes on the day/night variation.

The temperature effect has been highlighted by the observation of optical phase curves (Demory et al. 2013; Esteves et al. 2013, 2015) from which the distribution of clouds on the dayside atmosphere can be mapped (Shporer 2017). As shown by Shporer and Hu (2015), all the best targets for cloud mapping showed a partial cloud coverage on their dayside, indicating that most hot Jupiters are likely partially cloudy, with clouds present on the western part of the dayside and absent from the substellar point and the hot spot. As seen in Fig. 12, this cloud distribution is anticorrelated with the temperature distribution derived from infrared phase curves which points toward an eastward shift of the hottest part of the atmosphere (see Parmentier et al. 2018, for a review). This indicates that clouds are present where the atmosphere is cold and absent where it is hot. For the ultra hot Jupiters this can be easily explained by the fact that no known species can condense at the high temperatures of their dayside (Helling et al. 2019b). For cooler planets, however, something must prevent the clouds to form in the dayside. One possibility, proposed by Parmentier et al. (2016), is that the most refractory elements that could form clouds on the dayside of the cooler hot Jupiters are trapped in the deep layers of the planet. If that is correct, hotter planets would be expected to have silicate (MgSiO3 and Mg2SiO4), corundum (Al2O3), Iron or perovskite (CaTiO3) clouds whereas planets with an equilibrium temperature cooler than 1500 K would have an atmosphere dominated by sulphide clouds such as manganese sulphide (MnS), zinc sulphide (ZnS) or sodium sulphide (Na2S). Such a transition would be very similar to the transition between L and T brown dwarfs Morley et al. (2012). However, this mechanism operates only if cloud particles can be sequestered below the photosphere for cooler planets, which is a balance between settling and vertical mixing. The presence of such a deep sequestration of condensates has been simulating for silicate clouds by Powell et al. (2018). The study shows that part, but not all, of the cloud material was sequestered in the deep layers of the atmosphere. The formation of larger particle sizes, through the inclusion of a larger number of species such as iron could potentially help sequester a larger amount of silicates and explain the trends seen in the Kepler optical phase curves. Overall, determining more precisely the cloud composition of hot Jupiters would shed light into the deep atmosphere, including the rate of vertical mixing and the deep thermal structure, both related to the deep atmospheric circulation.

Fig. 12
figure 12

The optical phase curve of Kepler-41b observed with the Kepler space telescope (Shporer and Hu 2015) compared with the output from cloudless global circulation and models including the presence of silicate clouds of different particle sizes. The left inset shows the temperature distribution on the dayside of a global circulation model with similar equilibrium temperature (Parmentier et al. 2016). The right inset shows the brightness distribution on the dayside of the same model integrated over the Kepler bandpass. The bright crescent on the left is due to the reflection by clouds and creates the positive shift of the optical phase curves

The atmospheric circulation can also affect the cloud coverage by advecting cloud material from one part of the planet to another one. The level of coupling between flow and particles depends on the pressure. For a given particle, the settling timescale increases with decreasing pressure and, below some pressure, always becomes smaller than any circulation timescale. At 1 mbar in HD209458b, Parmentier et al. (2013) show that particles 1 μm or smaller should be coupled to the flow whereas larger particles should be more strongly affected by the vertical settling. As a consequence, simulations often show that the particles horizontal distribution is more homogeneous for small than for large particles. Particularly, for large particles, the interaction between the atmospheric circulation and the settling particles often lead to a peculiar behaviour at the equator. As seen in Fig. 13 in three out of four models the equator gets depleted in clouds whereas in the Lee et al. (2016) models the equator gets enhanced in clouds. Overall, more investigation is needed to understand the banding of clouds in hot Jupiters.

Fig. 13
figure 13

Cloud particle number density at 1 mbar of four different models using different assumptions. Parmentier et al. (2013) and Komacek et al. (2019) model HD209458b using passive tracers advected by the flow with a settling velocity corresponding to 2.5 μm particles. Whereas the former uses a non-grey radiative transfer scheme including TiO/VO, the second one uses a semi-grey radiative transfer. Lee et al. (2016) and Lines et al. (2018) use the same microphysical cloud model. The former is coupled to the GCM of Dobbs-Dixon and Agol (2013) with the parameters of HD189733b whereas the latter is coupled to the Metoffice ENDgame global circulation model of Mayne et al. (2014) with the parameters of HD209458b without TiO. Both include the radiative feedback of the clouds. All models show that the particle number density is different between the equator and the mid-latitudes, showing that banded clouds can be formed by the circulation. Models disagree, however, as to whether the equator should be depleted or enhanced in cloud particles. The dynamical mechanism responsible for this pattern is not yet fully understood

2.5.2 How Does the Cloud Distribution Affect the Atmospheric Circulation?

In hot Jupiter atmospheres, clouds affect the atmospheric circulation mainly through their optical properties, the latent heat release during their formation being negligible. When present on the dayside, they can cool the planet by reflecting some of the incoming stellar radiation but can also warm the atmosphere by increasing the greenhouse effect. When on the nightside, clouds always tend to warm up the atmosphere. The final effect of clouds on the atmospheric circulation therefore depends on their precise location. Roman and Rauscher (2017) show that nightside clouds tend to warm the planet nightside, leading to smaller day/night contrast on isobars. A global cloud, however, produces a temperature structure similar, but cooler than the one of a clear atmosphere. Although the cloud radiative feedback can strongly affect the day/night contrast and the hot spot offset on isobars, all models found that they do not change qualitatively the general behaviour of the atmospheric circulation such as the presence of a superrotating jet (see Parmentier et al. 2016; Lee et al. 2016; Lines et al. 2018 or Roman and Rauscher 2019).

Recently, Parmentier et al. (2020) carried a systematic study of the effect of non-grey nightside clouds on hot Jupiter atmospheres and their observational consequences using the SPARC/MITgcm. They confirm that nightside clouds dramatically reduce the overall day/night heat transport (see Fig. 8). Furthermore, they highlight an apparent contradiction: whereas on isobars nightside clouds decrease the day/night temperature contrast and increase the hot spot shift on isobars, they usually increase the phase curve amplitude and decrease the phase curve offset. This apparent contradiction is due to the fact that phase curves probe hemispherically averaged flux maps. When sharp variations in flux with longitude are present due to sharp variation of the atmospheric opacities, the center of the brightest hemisphere does not necessarily correspond to the longitude of the brightest point. Overall, this shows that phase curve amplitude and phase curve offsets do not necessarily probe the day/night temperature contrast and hot spot offsets on isobars. This should be kept in mind when interpreting the results from analytical theories that often assume isobaric heat transport.

Importantly, for a given planet, the effect of the clouds on the atmospheric circulation can easily be over or underestimated. Neglecting the wavelength dependence of the cloud opacity (such as Roman and Rauscher 2019), or the scattering properties of the clouds (Lee et al. 2016) can strongly enhanced the heating rates and affect the thermal structure and modelled observation in a much stronger way than in reality (Harada et al. 2019). Similarly, considering a species that is not present in the atmosphere can lead to a strong overestimate of the cloud radiative feedback. As shown by Lines et al. (2019), using two different cloud parametrisation can lead to very different radiative feedback, thermal structure and observational consequences.

2.5.3 Clouds and Atmospheric Variability

Atmospheric variability has been long sought in hot Jupiter infrared observations. Clear sky models of hot Jupiters are predicted to vary by less than \(1\%\) in their global temperature (Showman et al. 2009; Rauscher and Menou 2012b; Komacek and Showman 2020), which would not lead to observational evidence given current observational facilities. Optical phase curves, however, have been shown to vary significantly for two planets (Armstrong et al. 2016; Jackson et al. 2019). Although the source of the time variation is still under debate, and possibly includes magnetic coupling between the winds and the planetary magnetic field (Rogers and Komacek 2014; Rogers 2017) it has been postulated that the presence of clouds could significantly enhance the observable variability (Lines et al. 2018). Indeed, a planet dayside can rapidly evolve from barely cloudy to fully cloudy when the global temperature changes by one to two hundreds degrees (Parmentier et al. 2016). Moreover, on a given bandpass, the relative importance of the reflected light and thermal emission can be extremely sensitive to the cloud abundance and the temperature. Armstrong et al. (2016) postulated that the atmosphere of HAT-P-7b could oscillate between two states, one hot and relatively clear dominated by thermal emission and having a bright spot at the temperature maximum, east of the substellar point and a colder, cloudier state dominated by reflected light and having a bright spot at the cloud maximum, west of the substellar point. Such an oscillation would lead to large variation from positive to negative in the phase curve offset while keeping the total dayside luminosity relatively constant.

2.5.4 Haze

Haze modelling, assuming that methane is a main haze precursor, have shown that hazes can form and explain some features in the transmission spectrum of hot Jupiters Lavvas and Koskinen (2017). However, since all the observed planets do not have large methane concentrations, other pathways might be more important. Particularly, recent experiments at \(\approx 1500~\text{K}\) by (Fleury et al. 2019) have shown that photochemical hazes can form in a pure H2-CO gas at high temperature. Haze and condensation clouds are expected to be affected differently by the atmospheric circulation. Whereas cloud particles form and evaporate at similar temperature and pressures, hazes can have a much higher evaporation temperature than the temperature they formed. As a consequence, whereas clouds are expected to track the temperature maps of hot Jupiters, hazes might provide a more homogeneous aerosol cover. The two could potentially be distinguished by measuring the difference between the east and west limb of the planet (Line et al. 2016; Powell et al. 2019). Hazes would cover both the east and west limbs whereas condensation clouds would only be present on the cooler western limb (Kempton et al. 2017). Recent GCM investigations are starting to test these ideas (Steinrueck et al. 2019).

2.6 Chemistry

If the atmospheres of hot Jupiters were in local chemical equilibrium, the large day/night temperature contrast would naturally lead to large horizontal variation in their chemical composition. However, the consensus is that hot Jupiter atmospheres are not in chemical equilibrium. The reason is that the chemical timescale is an exponential function of both temperature and pressure. As an example, the time it takes to convert carbon monoxide to methane at 0.1 bar is of order of 10 days at 2500 K but of the order of millions of years at 1000 Kelvin (Visscher 2012). For comparison, the horizontal advection timescale is on the order of days. As first shown by Cooper and Showman (2006), atmospheric mixing is expected to homogenise the atmospheric abundances horizontally extremely easily. The exact value of the resulting abundances is, however, subject to the details of the atmospheric circulation. To first order, Agúndez et al. (2014) show that the dayside photospheric abundances set the nightside abundances of the atmosphere, with the dayside abundances being themselves set by the vertical quenching. As a consequence the final atmospheric abundances are determined by the vertical mixing strength on the dayside atmosphere and the pressure and temperature at the quench point, where the vertical mixing timescale equals the chemical timescale. Drummond et al. (2018a) and Drummond et al. (2020) additionally showed that the meridional circulation could also mix the chemical composition latitudinally, meaning that the 3D nature of the atmospheric mixing by the circulation should not be ignored. Finally, differences between chemical schemes (e.g. Moses et al. 2011; Venot et al. 2019) and the simplifications needed to couple the scheme to the atmospheric circulation (Tsai et al. 2018) can also change the expected quenched abundances (see Fig. 13 of Drummond et al. 2018b).

Changing the chemical abundances of the atmosphere through chemical quenching changes the optical properties of the atmosphere, hence the energy balance and the atmospheric circulation. However, as shown by Steinrueck et al. (2019), Drummond et al. (2018a,b), for the methane/carbon monoxide reaction, quenching affects the thermal structure of the atmosphere by \(\pm 100~\text{K}\), which is small compared to the day/night contrast of the atmosphere. Although important when comparing model outputs to observations, chemical quenching does not fundamentally alter the atmospheric circulation.

There are two situations where the chemical abundances of CO, CH4, and H2O are expected to depart from being vertically and horizontally constant. First, at pressures lower than \(\approx 10~\upmu \text{bar}\), photochemistry is expected to significantly deplete the dayside atmosphere from molecules such as methane, water or ammonia while increasing the abundances of other species such as HCN or CO2. Horizontal advection is expected to extend the effect of photochemistry to the limb and part of the nightside atmosphere (see Fig. 14). Importantly, photochemistry can also affect the higher pressures through vertical mixing, resulting on an inverted quenching as seen for the case of HCN in Fig. 14.

Fig. 14
figure 14

Abundances of several molecules at the equator of HD189733b calculated by Agúndez et al. (2014) using a 2D model (longitude/pressure). The abundances are determined by a combination of diffusive vertical mixing, advective horizontal mixing and photochemistry. At pressures larger than 1–10 bars, the abundances are in chemical equilibrium, between 1–10 bars and 100–10 μbar the abundances are determined by a combination of vertical and horizontal quenching. For pressures lower than 100–10 μbar photochemistry starts being important

The second situation where abundances are not constant through the atmosphere is extremely high temperatures, such as occur on ultra-hot Jupiters, when molecular dissociation becomes important. The timescale to thermally dissociate and reform molecules is much faster than the timescale for other chemical reactions (Bell and Cowan 2018; Kitzmann et al. 2018). As a consequence, ultra-hot Jupiters are expected to be in close to chemical equilibrium with respect to the dissociation reactions. The molecular abundances of almost all species, including H2, are expected to vary both vertically and horizontally, CO being the only molecule to not dissociate at these temperatures due to its higher binding energy (Parmentier et al. 2018; Kitzmann et al. 2018). In the hottest planets the dissociation of H2 on the dayside and its recombination on the nightside is expected to significantly enhance the day to night heat transport through latent heat transport (Bell and Cowan 2018; Komacek and Tan 2018; Tan and Komacek 2019; Mansfield et al. 2020; Wong et al. 2019).

2.6.1 Vertical Mixing

Vertical mixing is a fundamental outcome of atmospheric circulation. Understanding mixing is critical to understanding how the atmosphere interacts with clouds and chemistry. It affects the abundance of clouds, determines their particle size distribution and whether particles can sequester important chemical species deep in the atmosphere. However, understanding vertical mixing in hot Jupiter is a challenge. Despite being locally stable to convection, the atmospheres of these planets can transport and mix material with large scale atmospheric dynamics. Using passive tracers in a global circulation model of HD209458b, Parmentier et al. (2013) show that the vertical mixing in hot Jupiters is due to the combined effect of large scale upwelling on the dayside, downwelling on the nightside and localised updraft and downdraft close to the terminators. These mixing patterns are very different from known solar system and brown dwarfs equivalent, and, for example, vertical mixing calculations based on mixing length theory cannot be easily applied to the hot Jupiter case. As an example, for HD209458b Parmentier et al. (2013) estimated the vertical mixing by measuring the vertical flux of passive tracers along isobars of their global circulation model and found a vertical mixing coefficient that is 100 times smaller than the prediction from mixing length theory. Nonetheless, the resulting vertical mixing at the photosphere is much larger than at the photosphere of known solar-system planets (e.g., eddy diffusivity of \(\approx10^{5}~\text{m}^{2}/\text{s}\) at 100 mbar).

Earth stratosphere is probably the best solar-system analogue for vertical mixing in hot Jupiters. As shown by Holton (1986), mixing in the stratosphere is determined by large scale atmospheric motions. The formulation of Holton (1986) served as the base of the analytical models developed by Zhang and Showman (2018a,b) and Komacek et al. (2019). They show that the vertical mixing is tightly tied to the correlation between the abundance on isobars of the species being mixed and the vertical motions. When a chemical species has an equilibrium background abundance varying with pressure (for example, decreasing with decreasing pressure), the atmospheric circulation naturally causes a correlation between vertical motion and abundances. Where the winds are upward they carry parcel of gas with high chemical abundance, while where winds are downward they carry less chemical abundance. The net effect leads to a net chemical flux from regions with high chemical abundance to regions with low abundance. Efficiency of chemical transport by this correlation can be modulated by two different phenomena. The first one is horizontal mixing: the larger the isobaric winds, the harder it is to maintain a horizontal perturbation due to the vertical winds. The second one is relaxation timescale. Given enough time a chemical species abundance will tend towards local chemical equilibrium. If this equilibrium abundance is only pressure dependent, then the fastest the chemical timescale, the smaller the horizontal variations and thus the smaller the mixing. Following these ideas, Komacek et al. (2019) derived that the vertical mixing coefficient should scale as:

$$ K_{\mathrm{zz}}= \frac{w^{2}}{\frac{1}{\tau _{\mathrm{chem}}}+\frac{1}{\tau_{\mathrm{adv}}}}, $$
(14)

where \(w\) is the vertical velocity, \(\tau _{\mathrm{chem}}\) is the chemical timescale and \(\tau _{\mathrm{adv}}\) the horizontal advective timescale. This formula shows that a vertical mixing coefficient, \(K_{\mathrm{zz}}\), cannot be extracted from an atmospheric flow independently from the chemical species. Otherwise said, different chemical species are advected following different paths in the 3D atmosphere which are equivalent to different 1D vertical mixing coefficients. Finally, the formula is valid for species such as chemical species that are non-conservative. Cloud particles cannot be well represented by this model since they do not disappear but fall down. Intriguingly, despite having investigated settling timescales spanning several orders of magnitude, Parmentier et al. (2016), Zhang and Showman (2018b) and Komacek et al. (2019) all found that the vertical mixing coefficient derived from the global circulation models was rather independent of the particle size. More work is needed to understand this behaviour and derive a relationship estimating vertical mixing for settling particles.

Finally, because the vertical mixing timescale inherently depends on the horizontal mixing, any estimate of the vertical mixing rely on a global averaging on isobars (e.g. Parmentier et al. 2013). If one is lucky, the horizontally averaged circulation on isobars might act like a one-dimensional diffusive column. However, on local scales, mixing is not generally diffusive. As such, it is generally not correct to calculate local \(K_{\mathrm{zz}}\) profiles from 3D circulation. Doing so will likely produce spurious variations of mixing with height (e.g., the local vertical mixing coefficient would go to zero when the local vertical velocity changes sign) and the sensitivity of the outputs of the model to the method to calculate the vertical mixing coefficient should be thoroughly tested to understand which conclusions are robust and which ones are not (see e.g. Helling et al. 2019b, for a discussion).

2.7 Other Aspects

2.7.1 Choice of Equations

All the global circulation models described in the studies above do not solve the same set of dynamical equations. The primitive equations (used in e.g. the SPARC/MITgcm) are the standard dynamical equations appropriate for circulations in stratified atmospheres having horizontal length scales greatly exceeding their vertical scales (see Vallis 2006 or Showman et al. 2010 for a discussion of equation sets in atmospheric dynamics). These conditions of stratification and large aspect ratio are generally met for the large-scale flow in planetary atmospheres, including hot Jupiters, where the typical horizontal length scales are \(10^{4}\text{--}10^{5}~\text{km}\) and atmospheric scale heights are \(\sim 200\text{--}500~\text{km}\), leading to aspect ratios of \(\sim20\text{--}500\). Large aspect ratio and stable stratification generally allow the vertical momentum equation to be replaced by local hydrostatic balance, i.e., \(\partial p/\partial z=-\rho g\), where \(p\) is pressure, \(\rho \) is density, \(z\) is height, and \(g\) is gravity; this balance means that spatially and temporally varying meteorological perturbations to the density are generally in hydrostatic balance with meteorological perturbations to the vertical pressure gradient. Showman et al. (2008, 2018b) performed a scaling analysis of the vertical momentum equation for the large-scale flow on hot Jupiters, which suggests that, for large-scale flows, hydrostatic balance is valid locally to typically \(\sim1\%\) or better.Footnote 12 Importantly, the primitive equations make no assumptions that density variations are small, nor do they set any explicit limits on the wind speeds, in contrast to some other reduced equation sets in atmospheric dynamics.Footnote 13 However, the fastest wind speeds in hot-Jupiter circulation models commonly exceed the speed of sound, at least in local regions (the speed of sound in an ideal-gas hydrogen atmosphere is \(2.2~\text{km}\,\text{s}^{-1}\) at \(1000~\text{K}\), rising to nearly \(4~\text{km}\,\text{s}^{-1}\) at \(3000~\text{K}\)). This raises the question of whether acoustic shocks may form in the atmospheres of hot Jupiters. Shocks, by definition, are sharp features with small horizontal scales across the shock, which would necessarily imply a local violation of hydrostatic balance. Thus, an important question is the extent to which shocks are important and whether the primitive equations break down in a significant way. Furthermore, in fluid planets, the smooth connection between the deep interior and the observable atmosphere can lead to flows spanning a large vertical extent that can get close to violate the traditional approximation used in the primitive equations (Mayne et al. 2019).

To address this issue, several groups have performed 3D simulations of hot Jupiters using the fully compressible (Euler or Navier-Stokes) equations. Dobbs-Dixon and collaborators were the first to tackle this problem (e.g., Dobbs-Dixon and Lin 2008; Dobbs-Dixon et al. 2010, 2012; Dobbs-Dixon and Agol 2013). Their simulations produce circulations that are qualitatively very similar to those of primitive-equation GCMs, although differences in radiative forcing, implementation of friction, and other modeling details prevent a direct, head-to-head benchmark comparison.

Mayne et al. (2014) and Mendonça et al. (2016) presented new, non-hydrostatic hot Jupiter GCM codes, the former based on the UK Met Office dynamical core, and the latter a custom code called THOR, which both allow several levels of approximation to be made in the dynamical equations. The most sophisticated system is the fully compressible, non-hydrostatic Euler equations with radially varying gravity and no thin-shell approximation (“full” in the notation of Mayne et al. 2014); these are followed by equation sets that successively introduce various simplifications, including constant gravity (“deep”), thin-shell and traditional approximations (“shallow”), followed by the addition of hydrostatic balance (the primitive equations).12 The inclusion of multiple choices of approximation in a single code is extremely useful, because it allows head-to-head comparisons between the results of different equation sets, where the radiative forcing, domain, and other aspects of the numerics are held fixed. As a benchmark comparison, these studies both evaluated the “shallow hot Jupiter” test case defined in Menou and Rauscher (2009). Mendonça et al. (2016) showed that, when integrated with the “full” equations, this shallow hot-Jupiter benchmark yields very similar behavior to its primitive-equation implementation (Menou and Rauscher 2009; Heng et al. 2011b; Bending et al. 2013), as shown in the top row of Fig. 15. Mayne et al. (2014)’s “full” implementation of this test case did not show such good agreement, although subsequent work (Mayne et al. 2017) suggested that the discrepancies may have resulted from their model’s implementation of numerical diffusion rather than inherent differences between the primitive and fully compressible equations.

Fig. 15
figure 15

Comparisons of solutions from the primitive equations (left) versus full Navier-Stokes equations (right) in the equilibrated “shallow hot Jupiter” test case defined by Menou and Rauscher (2009), in which strong day-night heating (implemented with a Newtonian cooling scheme) occurs in an atmosphere with a surface pressure of 1 bar. Shown are zonal-mean zonal winds versus latitude and pressure. Left panel shows the primitive equation implementation of this test case from Heng et al. (2011b). Right panel shows the fully compressible, Euler-equation implementation of the test case, from Mendonça et al. (2016)

Mayne et al. (2017) presented further comparisons between primitive and fully compressible equations sets, for a hot-Jupiter setup extending deeper into the atmosphere. Day-night forcing was implemented with Newtonian cooling. No frictional drag was included near the base of the model, which implies that the winds spin up over time in these models. Comparisons between the primitive equations and the “full” set after 10,000 days of integration are shown in Fig. 16. Qualitatively, the two implementations produce very similar circulations, with roughly the same day-night temperature differences, qualitatively similar spatial temperature patterns, and maximum zonal-mean zonal winds that differ by only \(\sim10\%\). Still, a number of important quantitative differences occur, most notably the equatorial jet in the primitive equation variant (Fig. 16, left) extends considerably deeper than in the “full” variant (Fig. 16, right). Because these are not equilibriated—the zonal jet is continuing to spin up over time—it is unclear if these differences reflect true differences in the equilibrium state that would be achieved between the two equation sets, or rather simply reflect differences in the rates at which the zonal jet spins up in the two equation sets. Further comparisons and diagnostics of this type are needed to clarify the situation.

Fig. 16
figure 16

Comparisons of primitive equations (left) versus full Navier-Stokes equations (right) in spin-up models of hot Jupiters after 10,000 days of integration, from Mayne et al. (2017). Planetary parameters of HD 209458b are adopted. Day-night forcing is implemented with Newtonian cooling. Top row shows temperature (colorscale, K) and winds (arrows) at \(0.2~\text{bar}\); bottom row shows zonal-mean zonal winds (\(\text{m}\,\text{s}^{-1}\)) versus latitude and pressure

2.7.2 Hot Jupiter Inflation Mechanisms

Many hot Jupiters have observed radius larger than expectations from standard planetary evolution models (e.g., see Fig. 1 in Komacek and Youdin 2017). Strong stellar incident irradiation together with the assumption of efficient global homogenization of the incoming energy significantly reduce the interior cooling and help to sustain the large radius (e.g., Guillot et al. 1996; Fortney et al. 2007), but is still not sufficient to explain the observed radius. Various mechanisms involving external energy injection to the planetary interior or further reducing the interior cooling have been proposed. For thorough reviews of this topic, readers are referred to, for example, Fortney et al. (2010), Baraffe et al. (2014) and Dawson and Johnson (2018). Here we only briefly survey several hydrodynamic mechanisms related to the understanding of atmospheric dynamics of hot Jupiters.

If a small fraction (on the order of 1%) of the stellar irradiated energy is deposited near or below the radiative-convective boundary (RCB), the planetary contraction can be significantly slowed down or even halted (Guillot and Showman 2002; Komacek and Youdin 2017). Showman and Guillot (2002) and Guillot and Showman (2002) proposed that if significant vertical wind shear develops in the deep layers approaching the RCB where the stratification becomes weak, hydrodynamic instability such as the Kelvin-Helmholtz instability could occur and dissipate the kinetic energy associated with the mean flow into thermal energy via turbulent cascade near the RCB. A “Mechanical Greenhouse” mechanism has been proposed by Youdin and Mitchell (2010), in which forced turbulence (which may come from shear-wind instabilities) mix higher entropy in the outer radiative zone towards the interior, pushing the RCB to depth and significantly reducing the cooling. These mechanisms are qualitatively attractive, but the quantitative details involving the development and organization of small-scale turbulence and their interactions with the large-scale fields remain inconclusive. So far only very few studies exist in this direction (e.g., Fromang et al. 2016; Menou 2019), and future investigations unifying our understanding of the global-scale flow and small-scale (over lengthscales much smaller than grid sizes typically used in global models) turbulence are desired. Additional uncertainty of the above mechanisms also arises from our lack of understanding on the deep circulation of hot Jupiters (see Sect. 2.7.4).

Another mechanism driven by the atmospheric dynamics was proposed by Tremblin et al. (2017b). They developed a 2D radiative-advective non-grey model of hot Jupiters equator with a parametrized meridional and vertical mass transport. The simplicity of the model allows one to solve directly for the steady-state solution rather than integrating forward in time, which impedes most global circulation models to solve for a converged steady state. The model shows that for non-zero vertical transport, the deep atmospheric temperature pressure profile was converging towards a hot adiabat rather than an isotherm as in a one dimensional radiative-convective equilibrium (e.g. Fortney et al. 2008; Guillot 2010; Parmentier and Guillot 2014). When large enough vertical transport is assumed, the downward transport of energy is large enough to explain the inflated radii of hot Jupiters. This behavior was further explored with a 3D GCM by Sainsbury-Martinez et al. (2019) who used a Newtonian cooling scheme and integrated their models for several thousands of Earth years. They found that for large enough relaxation timescale for the temperature (e.g. larger then 3000 years at 200 bars) the model converges towards a hot adiabat as in the 2D model.

The large-scale circulation models proposed by Tremblin et al. (2017b) and Sainsbury-Martinez et al. (2019) is attractive in the sense that no (yet unknown) small-scale turbulence is required, but it has not been fully demonstrated. First, the mechanism behind the results is still unclear. To mechanically transport thermal energy from low to high pressure, one requires a net downward heat flux (\(\overline{\omega T'}>0\), where \(\omega \) is the vertical velocity in pressure coordinates, \(T'\) is the isobaric temperature variation relative to the global mean, and the overbar denotes a global mean at isobar surface) through a level near the upper deep layer given that the lower boundary is impermeable. In the deep layer, this is a thermally indirect circulation, meaning that the circulation works against buoyancyFootnote 14 (e.g., Holton and Hakim 2012). It cannot occur spontaneously within the deep layer by its own, and driving forces to the deep layers from the upper photosphere have to be involved. It is not yet clear what is driving the thermally indirect circulation (i.e., where \(\overline{\omega T'}>0\)) throughout the deep layer.

Secondly, both the 2D and the 3D models parametrize a fundamental part of the problem. The former parametrizes the vertical flow whereas the latter one parametrizes the heat transfer through a Newtonian relaxation scheme. Whereas both studies argue that their models are efficient to transport energy downward in most of the parameter space they explored, it is worth noting that the radiation and the circulation are intimately linked. Whether the actual physical parameters of hot Jupiters allow for the aforementioned thermally indirect circulation to develop is not clear. Particularly, a more recent study by Mendonça (2020) using a 3D GCM with semi-grey radiative transfer predicts a much smaller warming of the deep atmosphere then Sainsbury-Martinez et al. (2019), which would not be able to explain the inflated radius of hot Jupiters.

Finally, one often needs to be cautious when applying the simulated results which are in a closed domain to real hot Jupiters. The current simulations are often dependent on the choice of deep boundary conditions (see Sect. 2.7.4). For example, considering a more realistic situation wherein the GCM is connected to a cold planetary interior (equivalent to the assumed cold initial deep layers in Sainsbury-Martinez et al. 2019), the downward heat transport forced by the photospheric circulation could eventually trigger a sharp thermal inversion near the interface between the upper hot adiabat and the lower-entropy interior. Hence a strong and sharp stratification may suppress the circulation further down and stop the heat transport. Additionally, hot Jupiters likely have a convective interior and it is not yet clear how the presence of a radiative-convective boundary would affect the deep atmospheric circulation (see e.g. Rauscher and Showman 2014).

2.7.3 Magnetic Coupling

In the atmospheres of hot Jupiters, sodium, potassium, calcium and aluminium become ionized when dayside temperatures become higher than \(2000~\text{K}\) (Batygin et al. 2013; Helling et al. 2019a,b), rendering the atmosphere significantly conductive. Interactions between the magnetic field and the flow are therefore expected.

Magnetic fields are expected to interact with the atmospheric circulation in two ways. First they dissipate kinetic energy through Ohmic dissipation and hence slow down the winds. Second they affect the atmospheric circulation differentially in different directions, leading to a change in the circulation patterns. The dissipation effect was recognized early on and implemented in several GCMs by simply adopting frictional drag terms in the hydrodynamic equation sets (e.g. Perna et al. 2012). Further studies proposed that observational measurements, such as smaller than expected hot spot offsets, smaller than expected wind speeds and higher than expected day/night contrast could be caused by the presence of additional magnetic drag (Komacek et al. 2017; Kreidberg et al. 2018; Arcangeli et al. 2019; Koll and Komacek 2018).

In all the aforementioned studies, the additional drag terms tend to drive the atmospheric circulation from a jet-dominated towards a day-to-night flow Komacek and Showman (2016). However, as pointed out by Batygin et al. (2013) zonal flows should be much more stable than meridional flows when the magnetic field is aligned with the spin axis of the planet.Footnote 15

More self-consistent magnetohydrodynamic models were performed by Rogers and Komacek (2014). They confirmed that the presence of magnetic coupling leads to a reduced meridional flow. Furthermore, when the conductivity is allow to change due to large day-night temperature difference, Rogers (2017) show that the direction of the zonal jet can start to oscillate, leading to a variation in the hot spot offset of up to 20 degrees. Such a large variation in the hot spot offset may be responsible for the time variability observed in the Kepler phase curves of two hot Jupiters (Armstrong et al. 2016; Jackson et al. 2019).

The coupling between magnetic field and the circulation leads to the downward transport of stellar irradiated energy to planetary interior through the Ohmic dissipation. Stellar energy is deposited at the photosphere which fuels a vigorous atmospheric circulation. This circulation produces currents that can connect to deeper atmospheric layers and dissipate energy (Batygin and Stevenson 2010). If more than \(\sim 1\%\) of the stellar energy can be deposited near the RCB, this could explain that many hot Jupiters have an inflated radius (e.g. Guillot and Showman 2002; Menou 2012a; Thorngren et al. 2019). However, as shown by Rogers and Showman (2014) and Rogers and Komacek (2014), the simulations show that magnetic drag slows the deep winds significantly more than predicted by Menou (2012b), leading to Ohmic heating \(\sim100\) times too small to explain the inflated radii.

Additional magnetic effects in the atmospheres of hot Jupiters that have been explored theoretically but have not been included in GCMs include the presence of thermal instabilities above the photospheric layers triggered by local heat deposition by Ohmic dissipation Menou (2012c) or the presence of an atmospheric dynamo triggered by the large day/night variation in conductivity (Rogers and McElwaine 2017), which could change the magnetic field topography.

Finally, it is worth pointing out that no fully satisfactory model currently exists to compare to the observations of ultra hot Jupiters. The only MHD model that has been applied to hot Jupiter originates from stellar interior modeling from the Rogers group. The anelastic approximation used systematically leads to smaller wind speed than predicted by models solving the primitive or fully compressible equations. In addition, no radiative transfer scheme was used (only Newtonian cooling schemes) and therefore the model outputs cannot be directly compared to observations. On the other hand, models including non-grey radiative transfer are pure hydrodynamic models that only consider the magnetic effect through the inclusion of simple drag terms, which could lead to an incorrect atmospheric circulation pattern in very hot atmospheres. In the coming decades, coupling the two approaches will be a necessary step to interpret the observations of ultra hot Jupiters.

2.7.4 Deep Boundary Conditions and Integration Times

Gaseous planets lack a distinctive boundary between the atmosphere and the interior, and the choices of bottom boundary conditions in GCMs could affect the results both at the photospheric levels and in the deep layers. Most hot Jupiter GCMs adopted a slip-free, non-permeable bottom boundary at a pressure of around 100 or 200 bars (e.g., Showman et al. 2009; Heng et al. 2011a; Rauscher and Menou 2012a; Mendonça et al. 2016; Mayne et al. 2017). Compared to the shallow hot Jupiter simulations with a bottom pressure of about 1 bar (e.g., Heng et al. 2011b; Mayne et al. 2014), flows in the “deep” simulations are typically more time-invariant. Recently, Carone et al. (2020) proposed that for rapidly rotating hot Jupiters like WASP-43b, the simulated domain should extend to a deeper pressure (\(\sim700~\text{bars}\)) to properly capture dynamics in the observable layers.

The radiative timescale increases rapidly with increasing pressure due to the increased opacity and atmospheric mass. Convergence of the whole simulated domain is therefore bottlenecked by the extremely long radiative timescale in the deep layers. Recent GCM experiments by Wang and Wordsworth (2020) showed that several hundred simulated years are needed for the deep flow to converge if assuming an 80-bar bottom boundary pressure. Mendonça (2020) showed that several tens of simulated years are required for the equatorial jet to equilibrate. Such long integration time is yet challenging for GCMs coupled with non-grey radiative transfer, chemistry and cloud microphysics, and is also demanding on the conservation properties of dynamical cores. Proper simulating strategies should be investigated in the near future to cop requirements of both convergence and comprehensive physics.

In GCMs utilizing radiative transfer schemes, a horizontally isotropic net heat flux is typically applied at the bottom boundary. In this case the deep thermal structure is unconstrained and can evolve according to the dynamics driven by the upper atmosphere. In the case of inflated hot Jupiters, the deep model layer may have reached the convective zone given the likely high interior entropy (e.g., Thorngren et al. 2019). Scaling analysis and global convection models of planetary interior with uniform surface cooling predict small (∼ a few Kelvin) isobaric temperature variation in the convective zone (e.g., Stevenson 1991; Kaspi et al. 2009; Showman et al. 2011; Showman and Kaspi 2013). In the absence of strong molecular gradient, the interior is expected to be nearly adiabatic. This means that in the presence of convection, the zeroth-order temperature structure in the upper convective zone could be constrained by the specific interior entropy of the planet, which does not evolve over timescales relevant for GCM integration time. In addition, interactions between the stratified layer and the interior convection may trigger a wealth of turbulence and waves that propagate upward and interact with the mean flow (see Sect. 4). The extent of which the interior convection can affect the atmospheric circulation and deep thermal structure remains unexplored. In the other way around, it remains an open question whether and how interior convection will differ from those forced by uniform surface cooling when the hot-Jupiter-like atmospheric circulation is applied to the surface condition of interior convection models.

The interior of giant planets is expected to be electronically conducting, and interaction of the magnetic field with the flow leads to Ohmic dissipation and retards the flow (e.g., Liu et al. 2008). The large-scale flow in the interior is expected to be significantly slower than that near the photosphere. Convection leads to a nearly barotropic state of the interior, witch together with the slower winds (which implies a small global Rossby number) may result in a Taylor–Proudman effects (Pedlosky 2013) that tends to drag down winds in the deep GCM layers (e.g., Schneider and Liu 2009). To crudely represent this effect, some GCMs adopted a frictional drag near the bottom boundary that relax winds towards zero over characteristic drag timescales (e.g., Liu and Showman 2013; Komacek et al. 2017; Tan and Komacek 2019; Carone et al. 2020). Although easy to implement, the choice of drag timescale and the exact form of the drag are rather loosely chosen. In addition, recent gravity measurements on Jupiter by the Juno spacecraft have revealed that the zonal jets of Jupiter could penetrate down to \(\sim3000~\text{km}\) below the cloud deck (Kaspi et al. 2018; Guillot et al. 2018). Similarly, the zonal jets of Saturn may extend to \(\sim9000~\text{km}\) below the cloud deck which is constrained by gravity measurements from the Cassini spacecraft (Iess et al. 2019). The implication is that, for cooler planets in which the conducting region is far below the GCM domain, whether a deep drag that relaxes winds toward zero is questionable in general conditions.

Understanding the interactions and coupling between the photospheric level dynamics and interior dynamics is in a pressing need given the above issues. A self-consistent coupling between them is numerically challenging partly because traditionally two different sets of equations are used in different parts—primitive or fully compressible equations are used in the upper atmosphere and anelastic approximation is used in the interior. But more critically, in order to achieve a statistically steady state with realistic computational cost, global interior models are overly forced with non-dimensional dynamical parameters that are many orders of magnitude different to realistic values (e.g., Showman et al. 2011), while models of the upper atmosphere do not suffer from this issue. Proper modeling strategies are likely needed to circumvent this challenge. For example, without a fully coupling between the two parts, modeling of the individual part can be used as idealized boundary condition of the other part.

3 Warm Jupiters

A bit farther from their stars than hot Jupiters lie the warm Jupiters, which we define to be approximately Jupiter-mass planets with effective temperatures of \(\sim300\) to \(1000~\text{K}\). They lie just below hot Jupiters in the upper left corner of Fig. 2, with global-mean incident stellar and radiated IR fluxes ranging from a few hundred to a \(\text{few} \times 10^{4}~\text{W}\,\text{m}^{-2}\). As with hot Jupiters, their interior fluxes are expected to be relatively weak, perhaps 10 to 1000 times less than the incident stellar flux they receive (depending on the planet’s effective temperature, history, and other factors), implying that, like hot Jupiters, the warm Jupiters should have atmospheres that in the time-mean are nearly in radiative equilibrium with their parent star.

As yet, there are far fewer observational constraints on the atmospheres of warm Jupiters than hot Jupiters, but this situation may change in the future. Over 40 transiting warm Jupiters have been discovered. Although most of these are Kepler detections that are difficult to follow up due to their great distance from Earth, over a dozen have been discovered by groundbased surveys and the Transiting Exoplanet Survey Satellite (TESS) mission around brighter, closer stars amenable to further observational characterization.Footnote 16

The greater orbital separations of warm Jupiters from their stars will lead to key differences in behavior relative to hot Jupiters. In particular, the tidal effects that drive hot Jupiters into a (presumed) state of synchronization and that circularize their orbits are less dominant at the greater orbital distances of warm Jupiters. Therefore, unlike typical hot Jupiters—where it is generally a reasonable assumption that the rotation rate is equal to the synchronous value, and the obliquity and orbital eccentricity are zero—warm Jupiters should be expected to exhibit a range of rotation rates, obliquities, and orbital eccentricities. This may lead to a wider range of possible behavior.

Let us quantify the distances at which these transitions occur. The tidal synchronization timescale from a primordial rotation rate \(\Omega _{p}\) is (Guillot et al. 1996)

$$ \tau \sim Q\left (\frac{R_{p}^{3}}{G M_{p}}\right )\Omega _{p} \left (\frac{M }{M_{*}}\right )^{2} \left (\frac{a_{\mathrm{orb}}}{R_{p}}\right )^{6} , $$
(15)

where \(G\) is the gravitational constant, \(M_{*}\) is the stellar mass, \(a_{\mathrm{orb}}\) is the orbital semi-major axis, and \(Q\), \(R_{p}\) and \(M_{p}\) are the planet’s tidal dissipation factor, radius, and mass, respectively. Evaluating this expression using values appropriate for a typical hot Jupiter yields

$$ \tau \sim 1\times 10^{6}\left (\frac{Q}{10^{5}}\right ) \left (\frac{a_{\mathrm{orb}}}{0.05~\text{AU}}\right )^{6}~\text{yr} , $$
(16)

where we have used a solar mass for the star, a Jupiter mass and rotation rate for the planet, along with a planetary radius of 1.2 Jupiter radii, which is typical of hot Jupiters. For a canonical hot Jupiter 0.05 AU from its star and adopting \(Q\sim 10^{5}\) (an appropriate time-averaged value for the giant planets in the solar system), the equation yields a spindown time of \(\sim 10^{6}~\text{yr}\)—which underlies the common expectation that hot Jupiters should be close to a state of synchronous rotation. However, the synchronization timescale increases greatly with only modest increases in semi-major axis: for a planet 0.2 AU from its star, the spindown timescale increases to 4 Byr, comparable to system ages. Thus, EGPs beyond \(\sim 0.15~\text{AU}\) around sunlike stars will generally not be synchronized. They may thus exhibit a range of rotation rates and obliquities.

We now discuss, in turn, the effect of non-synchronous rotation, non-zero obliquity, and non-zero eccentricity on the atmospheric circulation of warm Jupiters, as currently understood. All of the research done so far assumes the incident stellar flux greatly exceeds the interior flux, which should be valid in relatively old systems for planets inward of \(\sim1~\text{AU}\). We close this section by offering a few remarks about the way these regimes may be modified by an interior convective flux.

3.1 Non-synchronous Rotation

Several GCM investigations have explored the effect of non-synchronous rotation under conditions of zero obliquity and eccentricity (Showman et al. 2009; Rauscher and Kempton 2014; Showman et al. 2015; Penn and Vallis 2017). Non-synchronous rotation exerts two effects—first the differing rotation rate changes the strength of the Coriolis parameter and its gradient \(\beta \), thereby influencing the Rossby number, the Rossby deformation radius, and other factors that depend on rotation. Second, non-synchronous rotation implies that the planet no longer has permanent daysides and nightsides, but rather that the dayside heating pattern sweeps in longitude across the planet over time. This effect might at first seem trivial, but in fact it can significantly alter the planetary wave modes that arise from the day-night heating pattern, with consequent implications for the superrotation, hot-spot offset, and other aspects of the circulation. Indeed, GCM experiments demonstrate that both effects critically influence the dynamical behavior.

Showman et al. (2009) performed simulations of HD 189733b at rotation rates 0.5, 1, 1.5, and 2 times the synchronous value, and Rauscher and Kempton (2014) performed a similar study for both HD 189733b and HD 209458b, with the goal of understanding how the two effects listed above influence the circulation regime and IR lightcurves. For HD 189733b, the rapidly rotating models in both studies showed the emergence of high-latitude eastward jets in addition to the primary superrotating equatorial jet. The slower rotating model exhibited a robust equatorial jet flanked by strong westward flow. In this case, the equatorial jet was fast and—despite the slower rotation—narrower than the equatorial jet in the synchronously rotating model. These changes do not follow the trends that occur in synchronously rotating GCM experiments when rotation rate alone is varied (e.g., Showman et al. 2008), which suggests that the longitudinal migration of the dayside heating pattern is critical in setting the detailed jet properties.

Strikingly, however, the slowly rotating HD 209458b experiment in Rauscher and Kempton (2014)—corresponding to a rotation period of 6.6 days—developed a qualitatively different circulation pattern comprising a strong westward equatorial jet, a dayside hotspot shifted westward of the substellar point, and day-night flow across most of the terminator. They report that equally slowly rotating models of HD 189733b (i.e., using a rotation period of 6.6 days, three times the synchronous value) also exhibit a similar circulation, with a fast westward equatorial jet rather than superrotation. The dynamical mechanism causing this transition remains unclear, but it appears that the slow rotation and migrating dayside heating patterns are key factors. Interestingly, Mendonça (2020) found in synchronously rotating hot-Jupiter models that rotation periods longer than 5 days allowed equilibrated states containing either eastward or westward equatorial jets, which may be a related phenomenon. As expected, synthetic IR lightcurves for the models with eastward jets reach maximum flux before the secondary eclipse, but the slowly rotating HD 209458b model exhibits an IR flux peak after secondary eclipse. Thus, this phenomenon presents a clear prediction for future lightcurve observations, as well as having a distinct Doppler-shift signature during transit (Rauscher and Kempton 2014).

Showman et al. (2015) suggested that the circulation of non-synchronously rotating hot and warm Jupiters systematically splits into two dynamical regimes, depending on whether the radiative time constant is shorter or longer than the solar day. If the radiative time constant is less than the solar day, the day-night heating pattern (diurnal cycle) is strong, and the circulation resembles the canonical hot-Jupiter regime already discussed: the day-night temperature differences are large, and the global wave modes triggered by the strong day-night forcing drive equatorial superrotation through the mechanism described in Sect. 2.

If instead the radiative time constant is greater than the solar day, the amplitude of day-night heating (the diurnal cycle) is expected to be weak. Longitudinal variations of temperature are therefore small, and the circulation is driven by the equator-to-pole gradient in zonal-mean heating, and primarily plays the role of transporting heat meridionally. At low obliquities, the sunlight predominantly illuminates the low latitudes, and the greatest temperature gradients tend to occur in mid-latitudes, with zonal-mean temperature decreasing poleward. Through thermal-wind balance (Eq. (4)), this implies that the strongest winds occur at midlatitudes rather than the equator, and furthermore that those winds are eastward at photosphere levels. The midlatitudes become baroclinically unstable, leading to baroclinic eddies that transport heat poleward; this triggers the formation of Rossby waves, which propagate meridionally away from their latitude of generation. The meridional propagation of the Rossby waves arising from these baroclinic eddies causes the convergence of angular momentum into the instability latitudes (for a review, see Vallis 2006; Showman et al. 2013b), promoting the generation of midlatitude, eddy-driven jets analogous to those observed on the Earth.

This transition tends to correspond approximately to a particular orbital distance inside of which the diurnal cycle is critical and beyond which it becomes less important or irrelevant. To calculate this distance, note that the effective temperature of a planet in energy balance with its star (and assuming zero albedo for simplicity) is

$$ T_{\mathrm{e}} =\frac{1}{\sqrt{2}}\left (\frac{R_{*} }{a_{\mathrm{orb}}}\right )^{1/2} T_{*}, $$
(17)

where \(R_{*}\) and \(T_{*}\) are the stellar radius and effective temperature, respectively. Inserting this expression into Eq. (3) for the radiative time constant yields

$$ \tau _{\mathrm{rad}} \sim \frac{p c_{p}}{g\sigma T_{*}^{3}}\left (\frac{a_{\mathrm{orb}} }{R_{*}} \right )^{3/2} $$
(18)

which implies that, for a given stellar type, the radiative time constant at the IR photosphere increases approximately as semi-major axis to the \(3/2\) power.Footnote 17 The solar day is \(P_{\mathrm{solar}} = 1/|P_{\mathrm{rot}}^{-1} - P_{\mathrm{orb}}^{-1}|\), where \(P_{\mathrm{rot}}\) is the rotation period and \(P_{\mathrm{orb}}\) is the orbital period. Equating the solar day to the radiative time constant and solving for the semi-major axis, we obtain

$$ a_{\mathrm{orb}} \sim P_{\mathrm{rot}}^{2/3} \left [\frac{1}{k} + \frac{g\sigma T_{*}^{3} }{p c_{p}}R_{*}^{3/2}\right ]^{2/3} $$
(19)

where we have used Kepler’s law, \(P_{\mathrm{orb}}^{2} = k^{2} a_{\mathrm{orb}}^{3}\) to relate the orbital period and semi-major axis, where the constant \(k=2\pi /\sqrt{G(M_{*} + M_{p})}\). For a sunlike star, and adopting typical planetary parameters (\(g\approx 20~\text{m}\,\text{s}^{-2}\), \(p=0.25~\text{bar}\), and \(c_{p}=1.3\times 10^{4}~\text{J}\,\text{kg}^{-1}\,\text{K}^{-1}\)), we find that the radiative time constant becomes longer than the solar day for orbital semimajor axes of 0.045 AU, 0.08 AU, 0.13 AU, and 0.2 AU for rotation periods of 0.5, 1, 2, and 4 days, respectively.

Thus, these arguments suggest that, around a sunlike star, planets inward of \(\sim 0.05\text{--}0.1~\text{AU}\) will exhibit strong diurnal cycles, but outward of \(\sim0.15\text{--}0.2~\text{AU}\), the diurnal cycle will be weak, and the circulation will be driven by the zonal-mean stellar heating.

Figures 17 and 18 provide a test of this prediction from Showman et al. (2015) for planets around an HD189733-like star, which is slightly dimmer than the sun. Simulations were performed on a regular grid in rotation rate and orbital semimajor axis. The condition (19) is just met for the lower left, middle, and upper right panels, respectively, of the two figures. The results indeed show the expected regime transition: models in the lower right exhibit essentially canonical hot-Jupiter-like behavior, while models in the upper left are in a regime where zonal temperature variations are minimal, baroclinic instabilities transport heat poleward in the midlatitudes, and the fastest jets occur in midlatitudes. The transition is broad and, not surprisingly, models near the predicted transition exhibit hybrid behavior with fast midlatitude zonal jets superposed on weak equatorial superrotation.

Fig. 17
figure 17

Temperature pattern over the globe (colors) and winds (arrows) for a set of 9 models straddling the transition between hot and warm Jupiters, illustrating a dynamical transition in the behavior, from Showman et al. (2015). Models in the left, right, and middle columns adopt rotation periods of 0.55 days, 2.2 days, and 8.8 days, respectively; the top, middle, and bottom columns adopt orbital semimajor axes of 0.2, 0.08, and 0.03 AU, respectively, for planets orbiting a K0 star with the properties of HD 189733. Models in the lower right have radiative time constants shorter than their solar days, and exhibit canonical hot-Jupiter circulation patterns with large day-night temperature differences and equatorial superrotation. Models in the upper right have radiative time constants longer than their solar days, and instead exhibit little day-night temperature variation; baroclinic instabilities transport heat toward the poles, and the fastest wind speeds tend to occur in midlatitudes rather than the equator. Model \(\text{H}\Omega _{\mathrm{med}}\) and \(\text{W}\Omega _{\mathrm{slow}}\) are synchronously rotating; the former has parameters identical to HD 189733b. The solid vertical bars in each panel show the substellar longitude at the time of these snapshots

Fig. 18
figure 18

Zonal-mean zonal winds versus latitude and pressure for the same 9 models shown in Fig. 17. There is a striking transition from strong equatorial superrotation for models in the lower right to mid-latitude eastward jet streams for models in the upper left. The fastest rotating, least irradiated model (upper left panel) exhibits multiple zonal jets in each hemisphere, and illuminates the dynamical continuum between the warm Jupiters and Jupiter and Saturn themselves. From Showman et al. (2015)

Interestingly, the fastest rotating, least-irradiated simulation of this set—in the upper left panel of Figs. 17 and 18—exhibits multiple eastward jets in each hemisphere, reminiscent of the pattern of multiple zonal jets on Jupiter and Saturn themselves. The zonal jets in this simulation exhibit meridional spacings that are qualitatively consistent with Rhines scaling, which implies that if the jet speeds were slower (as might be expected under even weaker stellar irradiation), then the jets would have smaller meridional spacings and they would be more numerous, becoming even more similar to Jupiter and Saturn. The transition of the jet structure along with the change in the rotation rate has been studied for conditions appropriate for terrestrial planets and qualitatively similar results to those in Fig. 18 were obtained (e.g., Williams 1988).

In the appropriate temperature range, a similar transition can occur among planets of differing atmospheric composition receiving a given flux of starlight. Lewis et al. (2010) showed that the atmospheric circulation regime of the hot Neptune GJ 436b, which has an effective temperature of 650 K, depends on the metallicity. Menou (2012a), Kataria et al. (2014) and Charnay et al. (2015) found just the same phenomenon for the super Earth GJ 1214b, whose effective temperature is 550 K. For these planets, atmospheres with metallicities \(\gtrsim30\) times solar have high opacities, and therefore low photosphere pressures, leading to short radiative time constants at the photosphere. This leads to strong day-night forcing at the photosphere, causing substantial day-night temperature contrasts and a fast equatorial jet whose maximum eastward winds are at the equator. On the other hand, at metallicities of a few times solar or less, the gas opacities are smaller, so the photosphere levels are deeper, leading to larger radiative time constants at the photosphere. These models have comparatively weaker day-night forcing, leading to temperatures that are nearly constant with longitude, with the fastest zonal winds occurring in midlatitude jets rather than at the equator. These two regimes are directly analogous to those occurring in Figs. 1617, as well as to that occurring in Fig. 9, and they occur for essentially the same reasons. The difference is that the change in radiative time constant that causes the transition is brought about not by a difference in the incident stellar flux but rather by a change in the photosphere pressure, via the metallicity.

3.2 Non-zero Obliquity

We next turn to consider the effects of non-zero obliquity. Only a few papers have so far addressed the influence on obliquity on hot and warm Jupiters (Langton and Laughlin 2007; Rauscher 2017; Ohno and Zhang 2019a,b), using both shallow-water and 3D models.

Let’s first consider some basic aspects. When obliquity is zero, of course, there is no seasonal cycle and the starlight predominantly irradiates low latitudes. When the obliquity is non-zero, the sunlight received by each hemisphere varies throughout the year, leading to seasonal effects, but additionally, the annual-mean irradiation changes, shifting annual-mean sunlight away from low latitudes toward high latitudes. For small obliquities, the strongest instantaneous (daily mean) starlight received remains at low latitudes, even near the time of summer solstice. However, when the obliquity reaches approximately \(18^{\circ }\), the daily mean insolation quickly shifts so that, at summer solstice, the maximum starlight is received at the summer pole—even though the annual-mean insolation still maximizes at low latitudes (Ohno and Zhang 2019a). This implies that, if the radiative time constant is less than the planet’s year, the seasonal effects will exhibit a large increase in amplitudes for obliquities crossing a threshold of \(\sim 18^{\circ }\). A second transition occurs when the obliquity exceeds \(54^{\circ }\), above which even the annual mean irradiation is greater at the poles than the equator. If planets with such large obliquities have short radiative time constants (relative to the planet’s year), they will of course have extremely strong seasons. If their radiative time constants are long compared to the year, they will lack strong seasons, but they will exhibit the unusual situation of exhibiting warm poles and a cold equator, with the atmospheric circulation transporting heat from high latitudes toward low latitudes.

These insolation changes can trigger rich transitions in dynamics, but basic aspects can be anticipated by simple arguments. At small obliquities, low latitudes tend to be hotter than high latitudes, and through thermal-wind balance (Eq. (4)), this poleward-decreasing temperature pattern implies that the midlatitude zonal winds should become more eastward with altitude. This is precisely the way Earth behaves, and is likewise consistent with the warm-Jupiter models shown in the upper left corner of Figs. 1718. On the other hand, if the poles are hotter than the equator, as expected at high obliquities and long radiative time constants, then thermal-wind balance implies that the midlatitude winds become more westward with altitude. If, during a given season, strong seasonal effects imply that one pole is warmer than the equator while the other pole is colder than the equator, then thermal-wind balance implies that the warmer (colder) hemisphere will exhibit midlatitude winds becoming more westward (eastward) with altitude.

Rauscher (2017) presented the first 3D study of how obliquity affects the circulation and observables of EGPs. They considered a rapidly rotating (10-hour period) giant planet on a 10-day period orbiting a sunlike star. Under these conditions, the planet’s effective temperature is \(880~\text{K}\); the rotation rate and irradiation level of their models are thus very similar to the “\(\text{W}\Omega_{\mathrm{fast}}\)” models of Showman et al. (2015) (see Figs. 17 and 18). Given the expectation that the diurnal cycle would likely be unimportant, Rauscher (2017) applied diurnally averaged (i.e., longitudinally axisymmetric) heating and investigated obliquities of \(0^{\circ }\), \(3^{\circ }\), \(10^{\circ }\), \(30^{\circ }\), and \(90^{\circ }\). Their low-obliquity simulations (\(\psi \lesssim 30^{\circ }\)) exhibit a warm equator, cold poles, and eastward zonal jets maximizing in the mid-to-high latitudes (Fig. 19), consistent with expectations from thermal-wind balance, and qualitatively in agreement with the rapidly rotating, weakly irradiated models of Showman et al. (2015).

Fig. 19
figure 19

GCM simulations from Rauscher (2017) showing the effect of obliquity on rapidly rotating warm Jupiters. Top, middle, and bottom rows show simulations with obliquities of \(30^{\circ }\), \(60^{\circ }\), and \(90^{\circ }\), respectively. Left column shows the zonal-mean temperature at the IR photosphere throughout one orbital period. Right column shows the zonal-mean zonal wind at the norther summer solstice

Rauscher’s high-obliquity (\(60^{\circ }\) and \(90^{\circ }\)) cases, however, exhibit a strong seasonal cycle and a very different circulation pattern. Near and after solstice, the summer hemisphere is much warmer than the winter hemisphere, and the hottest location on the planet is the summer pole. As expected from thermal wind balance, this sign reversal in the meridional temperature gradient—with high latitudes warmer than the low latitudes—leads to a sign reversal in the thermal-wind shear, and therefore the zonal flow is to the west. This effect is strongest at \(90^{\circ }\) obliquity, where the global flow is dominated by a single, broad westward zonal jet with speeds of \(\sim 1~\text{km}\,\text{s}^{-1}\). Figure 19 clearly shows that, over the majority of the seasonal cycle—including much of the winter—both poles are warmer than the equator. This behavior is consistent with the fact that, over an annual average, the poles receive more starlight than the equator for obliquities greater than \(54^{\circ }\), an effect which is modest for Rauscher’s \(60^{\circ }\) case but severe for their \(90^{\circ }\) case.

Figure 20 presents IR lightcurves resulting from Rauscher (2017)’s models for the particular situation of transiting planets (where, depending on the orientation of the planetary rotation vector relative to Earth, the sub-Earth latitude can lie anywhere between zero and a latitude equal to the obliquity). Generally, when the obliquity is small, or when the obliquity vector is large but in the sky plane (such that neither pole preferentially tilts toward Earth), the lightcurves are flat. However, if the rotation is oriented such that one pole aims toward Earth, the lightcurve exhibits significant variations, as the Earth-facing pole alternates between colder winters and warmer summers. If the planet’s rotation axis is tilted toward Earth (i.e., if the projection of the rotation vector into the planet’s equatorial plane is aimed directly at Earth), then the flux maxima occur after secondary eclipse, which is caused by the thermal lag whereby the summer hemisphere is hottest \(\sim1/8\) orbit after the solstice (see Fig. 19). If the rotation vector tilts in other directions, the situation is complicated, and the flux maxima of the lightcurve can occur either before or after secondary eclipse. The degeneracies in Fig. 20 imply that it will be hard to disentangle circulation behavior and obliquity for planets with unknown obliquity, but Rauscher showed that the acquisition of dayside flux maps from ingress/egress eclipse mapping can help break this degeneracy. A key takeaway point is that the flux peaks in Fig. 19 result from a completely different physical mechanism than those for close-in hot Jupiters.

Fig. 20
figure 20

IR lightcurves of oblique warm Jupiters from the GCM simulations of Rauscher (2017). These lightcurves are calculated for transiting planets, such that the sub-Earth latitude can lie be-tween zero and the obliquity. Different colors represent GCM simulations performed at different obliquities. For a given obliquity, different line styles represent different possible orientations of the planet’s rotation axis relative to the line-of-sight to Earth. Orbital phase is defined such that transit occurs at 0 and 1, while secondary eclipse occurs at 0.5

Ohno and Zhang (2019a) presented shallow-water models of warm Jupiters that explore the diverse range of regimes as a function of obliquity and radiative time constant, using an idealized model that represents the day-night forcing with a Newtonian cooling scheme in which the radiative time constant is an external parameter (Fig. 21). Unlike Rauscher (2017), they allowed a diurnal cycle, rather than adopting a diurnally averaged, axisymmetric forcing. They found five classes of behavior. As expected, when the radiative time constant is shorter than the stellar day (\(\tau_{\mathrm{rad}}=0.1~\text{day}\)), the diurnal cycle leads to a strong day-night thermal (height) contrast. When the radiative time constant is longer than the stellar day but shorter than the orbital period (\(\tau_{\mathrm{rad}}=5~\text{days}\)), the diurnal cycle is weak and the dynamical structure is more zonally symmetric; as expected, the seasonal cycle is weak when the obliquity is less than \(\sim 18^{\circ }\) but strong when it is greater than \(18^{\circ }\). When the radiative time constant is longer than the orbital period (\(\tau _{\mathrm{rad}}=100~\text{days}\)), then the seasonal cycle is weak at any obliquity; if the obliquity is less than \(54^{\circ }\), the height is thickest at the equator and thinnest at the poles, and the predominant winds are eastward. When the obliquity is greater than \(54^{\circ }\), the height is thicker at the poles than the equator, and the predominant winds are westward.

Fig. 21
figure 21

Regimes of behavior for warm Jupiters as a function of radiative time constant, length of day, and orbital period, from Ohno and Zhang (2019a). When the radiative time constant is less than the stellar day, the day-night contrast is large (regime I). At \(P_{\mathrm{rot}} < \tau _{\mathrm{rad}}< P_{\mathrm{orb}}\), the diurnal cycle is weak; when obliquity is less than \(\sim 18^{\circ }\), the seasonal cycle is weak, whereas when obliquity is greater than \(18^{\circ }\), it is strong (regimes II and III, respectively). When \(\tau _{\mathrm{rad}}>P_{\mathrm{orb}}\), the seasonal cycle is weak and the annual-mean insolation drives the circulation; for obliquities less than \(54^{\circ }\), the height is thicker at the equator than the poles and the winds are predominantly eastward, whereas obliquities greater than \(54^{\circ }\), the reverse is true, and the winds are predominantly westward

3.3 Orbital Eccentricity

EGPs with semimajor axes greater than \(\sim 0.05~\text{AU}\) exhibit a range of orbital eccentricities from zero to nearly one. In contrast to non-synchronous rotation and non-zero obliquities—which alter the distribution of starlight across the globe but not the total starlight hitting the planet—a non-zero eccentricity means that the total stellar power intercepted by the planet varies throughout its orbit. Many eccentric planets are transitional objects: the incident stellar flux received at periapse puts them in the hot-Jupiter regime, whereas the orbit-mean flux is more typical of warm Jupiters. Prominent examples amenable to observational characterization include the transiting planets HAT-P-2b (\(e=0.5\)), HD 17156b (\(e=0.67\)), HD 80606b (\(e=0.93\)), which have stellar fluxes that vary throughout their orbit by factors of 9, 27, and 828, respectively. The time variable forcing places eccentric EGPs in a novel regime that is both inherently interesting and may yield insights into atmospheric processes not easily obtainable from planets in circular orbits.

Generally, eccentric planets are not expected to exhibit rotation periods equal to their orbital periods. At large eccentricities, tidal effects are orders of magnitude stronger near periapse than throughout the rest of the orbit when the planet is farther from the star. Therefore, a common expectation is that for EGPs whose periapses lie very close to the star, tides act to drive eccentric planets toward “pseudo-synchronization,” a state in which the planet’s rotation rate equals the instantaneous rate of orbital motion when the planet is near periapse. Debate exists about the details, but for large eccentricities, this state corresponds to rotation periods much shorter than the orbital period. The rotation period is of course constant throughout the orbit, and given the varying rate of orbital motion, a pseudo-synchronous rotation state implies that, for high eccentricities, the stellar day will be much shorter near apoapse than near periapse.

Eccentric planets can be subject to several possible thermal regimes, due to the fact that their mean temperature—and therefore characteristic radiative time constant at the photosphere—vary throughout the orbit. The first regime corresponds to the situation where the radiative time constant is shorter than the stellar day throughout the orbit, so that the diurnal cycle (day-night radiative forcing) is important at all orbital phases. A second regime occurs when the periapse is small enough that, near periapse, the radiative time constant is shorter than the stellar day, whereas near apoapse it is longer than the stellar day. In this regime, the diurnal cycle is important near periapse, leading to an intense burst of day-night radiative forcing, but the diurnal cycle is not important throughout the rest of the orbit. A third regime occurs when the entire orbit—including the periapse—is sufficiently far from the star that the diurnal cycle is irrelevant throughout the orbit. So far, most eccentric EGPs whose atmospheres are being characterized lie in the first or second regime.

Several investigations of the atmospheric circulation of eccentric EGPs have been performed (Langton and Laughlin 2008; Kataria et al. 2013; Lewis et al. 2010, 2014, 2017). Langton and Laughlin (2008) adopted an idealized 2D model, whereas the other studies adopted a 3D GCM with sophisticated non-grey radiative transfer. For eccentricities that are small—less than perhaps 0.2—the effects of eccentricity on the qualitative circulation regime are modest (Lewis et al. 2010; Kataria et al. 2013). But for larger eccentricities, the investigations performed to date all agree that the intense stellar irradiation near periapse passage plays a dominant role in driving the circulation. This “flash heating” causes transient formation of a Matsuno-Gill-like pattern of standing equatorial waves that can drive equatorial superrotation, which can persist throughout the orbit even though it is not strongly forced throughout the more distant parts of the orbit (Kataria et al. 2013; Lewis et al. 2014, 2017). The flash heating event can also trigger shocks and other wave fronts that propagate from day to night. Still, the 3D studies performed to date indicate that this time-variable forcing appears to produce a time-variable version of the circulation regime already familiar from hot Jupiters on circular orbits, as opposed to a qualitatively different, totally new circulation regime. As mentioned, these studies tend to be in the first or second regimes described above, and one could imagine that eccentric EGPs in the third regime would lack Matuno-Gill patterns and would not have strong equatorial jets (at least not via the mechanism described in Sect. 2). Future work is needed to explore this issue.

For large orbital eccentricities and small periapse distances, the pseudo-synchronous rotation rate is fast, and this leads to an equatorial jet much narrower than occurs in GCM experiments of hot Jupiters on circular orbits. Kataria et al. (2013) performed a systematic investigation across a wide range of semi-major axes and eccentricities, exploring both pseudo-synchronous and synchronous rotation rates, and she showed that the meridional half-widths of the equatorial jet are reasonably well matched by the equatorial deformation radius, as predicted by the analytical theory (see Sect. 2.3).

These studies have important implications for observations of eccentric EGPs. The increase in temperature caused by periapse passage tends to lag the pulse of stellar insolation, such that the maximum temperatures occur after the planet is already past periapse and receding farther from its star. The overall timescale for this thermal pulse to damp out provides an observational measure of the atmosphere’s radiative time constant, a quantity than is difficult to infer for planets on circular orbits. The large day-night temperature difference can persist for times considerably past periapse passage, and if the rotation period is short, this can lead to a “ringing” phenomenon where the flux received at Earth oscillates up and down as the hotter hemisphere rotates in and out of view (Langton and Laughlin 2008; Cowan and Agol 2011a; Kataria et al. 2013; Lewis et al. 2017).

Synthetic lightcurves show that the observational signatures of transiting eccentric EGPs can be subtle and depend on the orientation of the planet’s orbit relative to Earth. Temporal flux variations in IR lightcurves reflect both spatial effects—in which hot and cold regions on the EGP move in and out of view as the planet rotates—and temporal effects in which the entire planet heats and then cools as it passes through periapse. Disentangling these effects can be difficult. Depending on the orbital alignment relative to Earth, the secondary eclipse can occur either before or after periapse passage, and, partly as a result, the peak IR flux can either lead or lag the secondary eclipse (Kataria et al. 2013). In most cases, the peak IR flux lags the periapse passage, although the timing depends on the orbital configuration, since the orbital alignment controls the time at which the newly flash-heated hot hemisphere first rotates into view as seen from Earth (Kataria et al. 2013). As a result of these issues, the range of possible lightcurve signatures is diverse.

Spitzer IR lightcurves have been obtained for HAT-P-2b over its full 5-day orbit (Lewis et al. 2013) and for HD 80606b near its periapse passage (Laughlin et al. 2009; de Wit et al. 2017). These observations show a large spike in the IR flux due to flash heating, which peaks around the time of periapse passage and then slowly decays. In general, this flux spike is reproduced in GCM simulations of these two planets (Lewis et al. 2014, 2017), but the detailed way the flux decays in the simulations exhibits discrepancies with the observations, especially at times \(\gtrsim 1~\text{day}\) after periapse, which suggests possible influences of clouds, chemistry, or interior heat flux that are not yet self-consistently included in the simulations (Lewis et al. 2017).

4 Brown Dwarfs and Directly Imaged Planets

Brown dwarfs are fluid hydrogen objects thought to form like stars but that contain insufficient mass to fuse hydrogen. They range in mass from approximately several Jupiter masses (\(M_{J}\)) up to the stellar mass limit of \(\sim80~M_{J}\).Footnote 18 Despite their large mass range, the hydrogen equation of state implies that these objects have radii close to Jupiter’s radius over nearly two orders of magnitude in mass, from less than \(1~M_{J}\) to \(80~M_{J}\) (e.g., Guillot 1999; Chabrier et al. 2000). Their lack of an internal heat sourceFootnote 19 implies that, like Jupiter, brown dwarfs cool down over time. Nevertheless, massive brown dwarfs contain so much internal energy that they cool down very slowly and, even after billions of years, can exhibit temperatures at the IR photosphere exceeding \(1000~\text{K}\). Brown dwarfs were first definitively discovered in the mid 1990s (Nakajima et al. 1995; Oppenheimer et al. 1995), around the same time as the first EGPs were discovered around main sequence stars (Mayor and Queloz 1995), and over the past 25 years the field of brown dwarf discovery and characterization has developed in parallel with the field of exoplanets. Because brown dwarfs are typically isolated—with no nearby star to drown out their light—they are generally much easier to observationally characterize than exoplanets. The quality and quantity of IR emission spectra for brown dwarfs is therefore exquisite by the standards of exoplanets.

Like stars, brown dwarfs have been subdivided into spectral types, the L, T, and Y types, based on their spectral characteristics (Kirkpatrick 2005; Cushing et al. 2011; Kirkpatrick et al. 2012; Leggett et al. 2017). L dwarfs are the hottest, with effective temperatures \(\sim2100\text{--}1300~\text{K}\), as well as reddish colors in the near-IR and rather shallow spectral features indicating the prevalence of cloud decks of silicates and other refractory condensates. T dwarfs are cooler, with effective temperatures of 1300 to \(\sim600~\text{K}\). Their IR spectra exhibit deep spectral features of water and methane, with blue near-IR colors, indicating comparatively cloud-free atmospheres. Y dwarfs are the coldest class, with effective temperatures less than \(\sim600~\text{K}\) and deep spectral features of water, methane, and other species. The disappearance of silicate clouds in T and Y dwarf spectra makes sense, because at their relatively cool temperatures, refractory materials like silicates condense only very deep in their atmospheres, below the optical-depth unity level, where the cloud layers cannot strongly affect their outgoing IR spectra. Since brown dwarfs cool down over time, this spectral sequence is also an evolutionary sequence, meaning a given object can transition from L to T to Y over billions of years. This process happens more quickly for lower mass objects, such that for objects of a given age, a massive brown dwarf may remain an L dwarf while a lower-mass object may have already transitioned to a T or Y dwarf.

A key point is that, starting in the L dwarfs, and proceeding through the T and Y dwarfs, the atmospheric temperatures of brown dwarfs are sufficiently cool that they should magnetically decouple from any background dipole magnetic field they may possess (e.g., Gelino et al. 2002), making them less like stars, and strengthening the analogy between the atmospheric dynamics of brown dwarfs and those of solar system planets like Jupiter and Saturn.

In addition to brown dwarfs, a variety of young, hot EGPs are being directly imaged and characterized. These include the planets around HR 8799, \(\beta \) Pic, and 2M1207, among other systems. Because of the difficulty of discerning the planet from the much-brighter star, planets that can be directly imaged tend to be distant from their stars (10 AU or more) and have photospheric temperatures of \(\sim500~\text{K}\) or more. This implies that they receive negligible stellar flux compared to the interior heat flux they radiate to space; therefore, from a meteorology standpoint, directly imaged EGPs resemble low-mass, low-gravity versions of brown dwarfs, and they will likely exhibit analogous dynamical processes.

Despite exhibiting similar temperature ranges, the regime of brown dwarf atmospheres differs significantly from that of hot Jupiters. Most known brown dwarfs are isolated objects that receive no starlight; they are hot not due to irradiation, but because they are massive and still losing their formidable heat of formation. This implies that, unlike typical hot Jupiters, they transport an enormous internal heat flux through their interiors and into their atmospheres to be radiated to space. For example, a 1500-K brown dwarf has an IR heat flux of \(290{,}000~\text{W}\,\text{m}^{-2}\), while even a 1000-K object radiates an energy flux of \(60{,}000~\text{W}\,\text{m}^{-2}\), hundreds of times greater than Earth’s IR flux to space of \(240~\text{W}\,\text{m}^{-2}\). In the interiors, this prodigious heat flux is expected to be transported by convection, but the transition to optically thin radiative transfer near \(\sim1~\text{bar}\) pressure decreases the radiative equilibrium temperature gradient, ensuring that the temperature structure transitions from nearly adiabatic in the deep interior to a stably stratified atmosphere at low pressure, with the radiative-convective boundary expected to occur typically at a few bars pressure (Fig. 22, Stevenson 1991; Burrows et al. 2001; Baraffe et al. 2014).

Fig. 22
figure 22

Modeled 1D temperature-pressure profiles of brown dwarfs in radiative-convective equilibrium, from Burrows et al. (2006). Twelve models are shown, with effective temperatures ranging from 900 to 2000 K in steps of 100 K (solid curves). On each curve, red represents the convection zones and black represents stably stratified layers in radiative equilibrium. Black dots show the radiative-convective boundaries. Stars show the approximate IR photosphere level where the bulk of the radiation escapes to space. Note that in most cases, the IR photosphere occurs in the stably stratified atmosphere. Understanding observations of IR spectra, variability, surface patchiness, and chemical disequilibrium therefore requires an understanding of the atmospheric structure, including the ability of atmospheric dynamics to mix clouds and chemical tracers in this stratified region

Several lines of evidence indicate that brown dwarfs exhibit vigorous atmospheric circulations. First, the IR spectra and colors of many brown dwarfs—particularly L dwarfs—indicate the presence of silicate clouds (dust) in the visible atmospheres (e.g., Chabrier et al. 2000; Tsuji 2002; Knapp et al. 2004; Kirkpatrick 2005; Cushing et al. 2006). In the absence of dynamics, cloud particles would settle gravitationally, so the existence of clouds implies the presence of atmospheric motions to loft the particles vertically. Second, the transition between L and T dwarfs is a puzzling phenomenon that seemingly requires a role for atmospheric dynamics, either in generating cloud patchiness or in influencing how cloud microphysics depends on spectral type (Ackerman and Marley 2001; Burgasser et al. 2002; Knapp et al. 2004; Burrows et al. 2006; Marley et al. 2010). Third, many brown dwarfs exhibit chemical disquilibrium of CO, CH4 and NH3, which seems to require vertical mixing of air between deep levels and the atmosphere (Fegley and Lodders 1996; Saumon et al. 2000, 2006; Leggett et al. 2007; Stephens et al. 2009; Miles et al. 2020). Fourth, many brown dwarfs exhibit IR variability over rotational timescales, which probably results from regional patches of differing temperature and cloud opacity—and therefore differing IR fluxes to space—moving in and out of view as the brown dwarf rotates (left panel of Fig. 23, Artigau et al. 2009; Radigan et al. 2012; Apai et al. 2013; Wilson et al. 2014; Buenzli et al. 2014; Buenzli et al. 2015; Metchev et al. 2015; Yang et al. 2016; Lew et al. 2016; Vos et al. 2019; Vos et al. 2020; Zhou et al. 2020; Bowler et al. 2020; for reviews, see Biller 2017 and Artigau 2018). Interestingly, the shape of the lightcurves changes on short timescales, indicating that the spatial pattern of patchiness changes over time—presumably due to time-evolving organization of turbulence, vortices, or other atmospheric structures over the globe. In some situations, the lightcurves can be used to place constraints on the size and spatial distribution of atmospheric spots on the globe (Karalidi et al. 2016; Apai et al. 2017). Fifth, Doppler imaging of the closest brown dwarf to Earth, Luhman 16B, provides direct confirmation of a patchy surface structure (Crossfield et al. 2014, right panel of Fig. 23), a technique that may be extended to other brown dwarfs in the future.

Fig. 23
figure 23

Observations indicating surface patchiness on brown dwarfs. Left panel shows the IR flux versus time measured in J-band (\(\sim 1.2~\upmu \text{m}\)) for brown dwarf SIMP0136 by Artigau et al. (2009). The J-band flux varies by about 5% with a period of 2.4 hours, indicating that regional patchiness in clouds and temperature move in and out of view as the brown dwarf rotates every 2.4 hours, causing significant variability in a disk average. As expected, the period is constant throughout the observing sequence, but the shape of the lightcurve varies substantially on intervals of several days, which implies that the global pattern of cloud patchiness evolves on this timescale. Right panel shows global maps throughout a rotation period of the brown dwarf Luhman 16B obtained by Doppler imaging from Crossfield et al. (2014). The maps indicate regional patchiness on a scale of tens of thousands of km. Note that although the patchiness is real, the preferential north-south elongation of features seen in the maps, as well as some of the lower-amplitude features, are probably analysis artefacts; moreover, at the signal-to-noise ratio of these maps, the method may not be sufficiently reliable to retrieve a zonally banded structure, if any (Crossfield et al. 2014)

Several authors have attempted to constrain atmospheric wind speeds from observations. Allers et al. (2020) recently presented the first true atmospheric speed measurement for a brown dwarf. They independently measured the rotation period of IR variability and the period from radio emission for the brown dwarf 2MASS J1047, showing that the former is about one minute shorter than the latter on this \(\sim1.75\)-hour-period object—a 1% difference. The IR period senses atmospheric cloud patchiness, while the radio period presumably senses the magnetic field rotation rate, which is rooted in the deep interior. Thus these measurements imply that the atmospheric features move eastward relative to the deep interior at a velocity of \(600\pm 300~\text{m}\,\text{s}^{-1}\). This phenomenon is most easily explained by advection of atmospheric features by a fast eastward zonal jet, though a contribution due to zonal wave propagation is not ruled out. In fact, this method is essentially the same approach used to infer the atmospheric wind speeds on Jupiter relative to those in the Jovian interior.

The number of brown dwarfs for which radio and IR periods can both be measured is limited, however, and so several authors have attempted to infer wind speeds solely from IR lightcurves. This requires plausibility arguments about the nature of the atmospheric features causing the variability. Artigau et al. (2009) and Radigan et al. (2012) found that time evolution of their lightcurves could be explained by the assumption that discrete, coherent atmospheric features change in longitudinal separation over time due to advection by an assumed zonal wind; the derived wind speeds are \(\sim300\text{--}500\) and \(30\text{--}50~\text{m}\,\text{s}^{-1}\), respectively. Similarly, Apai et al. (2017) inferred a speed of \(\sim 600~\text{m}\,\text{s}^{-1}\) but suggested this may result from differential propagation of atmospheric waves, rather than advection by zonal jets. Burgasser et al. (2014) and Karalidi et al. (2016) used lightcurves to place constraints on the diameter of atmospheric spots; under the assumption that the spots may coexist with jets of equal width and that the jets obey Rhines scaling, the inferred wind speed is \(\sim600~\text{m}\,\text{s}^{-2}\) or faster. These approaches are worth exploring, but—as the above authors are careful to admit—the assumptions they invoke are non-unique and uncertain, so these estimates should be taken as plausibility inferences rather than tight observational constraints.

Here, we survey the atmospheric dynamics of brown dwarfs, first using basic arguments to sketch the expected dynamical regime, and then summarizing current atmospheric circulation models of brown dwarfs, emphasizing physical insights, basic processes, and observational implications derived from these studies.

4.1 Basic Dynamical Regime

Isolated brown dwarfs rotate rapidly, and this exerts a strong control over their atmospheric dynamics. Doppler broadening of atmospheric spectral emission lines (Reiners and Basri 2008), and rotation periods directly inferred from IR variability (e.g., Artigau et al. 2009; Radigan et al. 2012; Metchev et al. 2015), indicate rotation periods typically ranging from 1.5 to 12 hours. Directly imaged planets likewise rotate rapidly, with periods of order several to 11 hours (Snellen et al. 2014; Zhou et al. 2016).

Because of the fast planetary rotation, the regional-to-global circulation on brown dwarfs and directly imaged planets exhibits Rossby number \(Ro \ll 1\), implying that the large-scale dynamics are rotationally dominated (Showman and Kaspi 2013). Adopting wind speeds of \(50\text{--}1000~\text{m}\,\text{s}^{-1}\) and rotation periods of \(1\text{--}10~\text{hr}\) consistent with the observationally inferred ranges, global-scale flows (\(L\sim 10^{8}~\text{s}\)) yield Rossby numbers of 0.0002 to 0.04.Footnote 20 If the dominant flow length scale is instead a regional scale (\(L\sim 10^{7}~\text{m}\)), similar to the dominant length scales on Jupiter and Saturn, then the Rossby numbers range from 0.002 to 0.4. Over almost the entire range considered, these values are much less than one, implying that the dominant horizontal force balance is between Coriolis and pressure-gradient forces, that is, geostrophic balance (Vallis 2006; Holton and Hakim 2012). The large-scale circulation on Earth, Jupiter, Saturn, Uranus and Neptune are also geostrophically balanced, and therefore a substantial body of literature in the dynamics of rapidly rotating turbulence developed over the past 40 years will likely have important insights for understanding brown dwarfs. Note for a typical brown dwarf with a rotation period of a few hours, the “tropical regime” corresponding to \(Ro\gtrsim 1\) obtains only within a few degrees of the equator.

A second influence of rotation is that it leads to small Rossby deformation radii. Under brown dwarf conditions, the deformation radius in midlatitudes is \(L_{D} = NH/f \sim 500\) to \(5000~\text{km}\)—just a few percent of the planetary radius—for brown dwarfs with rotation periods of 1 to 10 hours. The Rossby deformation radius is a typical length scale that emerges from a variety of processes involving an interplay between buoyancy and rotation. For example, baroclinic instability and convection tend to produce eddies with sizes comparable to the deformation radius. On Jupiter and Saturn, this fact helps to explain the large population of vortices with sizes of \(\sim 1000\text{--}2000~\text{km}\) (e.g. Mac Low and Ingersoll 1986; Li et al. 2004; Vasavada et al. 2006; Choi et al. 2009): they are probably a natural outcome of convection interacting with the overlying, stratified atmosphere. This argument suggests that brown dwarfs will likewise exhibit a large population of eddies with sizes of order \(1000~\text{km}\).

Still, the processes that drive a circulation in brown dwarf atmospheres differ significantly from those on typical irradiated planets. On hot Jupiters, Earth, and other irradiated planets, the large-scale atmospheric circulation results primarily from the gradient in stellar irradiation between the dayside and nightside (or between the equator and pole on planets where the diurnal cycle is weak). But most brown dwarfs lack external source of irradiation, and so this familiar way of generating an atmospheric circulation cannot occur on a brown dwarf. Convective mixing should lead to a nearly constant entropy in their interiors, and if so, the temperature on an isobar near the radiative-convective boundary should vary little with latitude. As a result, one might expect the upward IR flux emitted to space—and therefore the atmospheric temperature-pressure profiles in radiative equilibrium with that upwelling flux—to be nearly independent of longitude or latitude. In such a situation, the stratified atmosphere would exhibit essentially no horizontal temperature contrasts at all.

Still, as we will see, several mechanisms exist for generating a circulation in the stratified atmosphere of brown dwarfs. Because the IR radiation that escapes to space generally originates from within the stratified atmosphere rather than the convection zone (see Fig. 22), understanding the dynamics in this stratified region is critical to explaining observations of lightcurve variability, patchy cloud structure, and chemical disequilibrium.

4.2 Circulation Models of Brown Dwarfs

There exist at least three independent mechanisms for generating an atmospheric circulation in the stably stratified atmosphere of an isolated brown dwarf: (1) interaction of the interior convection with the stratified atmosphere; (2) 1D cloud feedbacks between radiation, vertical mixing, and cloud physics; and (3) multi-D cloud-radiative-dynamical feedbacks that drive an overturning circulation that maintains cloud patchiness. We summarize recent models investigating each of these three mechanisms.

4.2.1 Atmospheric Circulation I: Interaction of Convection with the Stably Stratified Atmosphere

In a pioneering study, Freytag et al. (2010) presented 2D numerical models of convection in a local box, representing a small patch of a brown dwarf, aimed at determining how convection interacts with an overlying, stably stratified atmosphere. Rotation was neglected, and the typical domain size was \(300\text{--}400~\text{km}\), about 0.5% of the radius of a brown dwarf. They coupled the 2D dynamics to a radiative transfer scheme and a dust cycle that computes the transport of a condensable gas species and its condensate—silicate dust—as well as the interconversion between them. In the simulations, convective plumes drip off the top of the convection zone and descend into the interior. The convection interacts with the overlying stratified layer, generating a wide population of small-scale gravity waves (Fig. 24). These waves break, triggering localized regions of vertical mixing in the stratified atmosphere, which can transport the clouds and condensable vapor upward, and generate cloud patchiness on a scale of tens of km. Freytag et al. (2010) analyzed the characteristic velocities in the convection zone and stratified layer, characterized how they depend on effective temperature, and derived effective vertical eddy diffusivities for upward mixing of tracer due to the breaking gravity waves. This helps to explain the cloud decks on typical L dwarfs, and may be relevant for understanding the L/T transition, but its lengthscale is too small to generate the lightcurve variability shown in Fig. 23. The 2D geometry and neglect of rotation also implies that the convection in their models cannot generate Rossby waves and has nothing to say about whether zonal jets, large-scale coherent vortices, or large-scale cloud patchiness are expected on brown dwarfs.

Fig. 24
figure 24

Entropy anomalies (top) and velocity and dust distribution (bottom) versus horizontal distance and height in a local 2D box model of convection from Freytag et al. (2010). This simulation corresponds to a brown dwarf with effective temperature \(1858~\text{K}\) and \(g=1000~\text{m}\,\text{s}^{-2}\). The convective zone roughly corresponds to \(z<0\), and convective plumes can be readily seen there. They interact with the overlying stratified layer at \(z>0\) to generate small-scale gravity waves

Showman and Kaspi (2013) presented the first global-scale models of the atmospheric circulation on brown dwarfs. They demonstrated that rotation plays a key role in the dynamics, leading to behavior on large scales (\(\gtrsim 1000~\text{km}\)) greatly differing from the small-scale convection simulated by Freytag et al. (2010). Showman and Ingersoll (1998) presented 3D spherical-shell, global models of convection in the interior to understand the global-scale convective organization. They solved the Navier-Stokes equations subject to the anelastic approximation, which are valid when the dynamical perturbations to the background density are small and the wind speed is much less than the speed of sound, as expected to be appropriate in the interiors of brown dwarfs. They demonstrated that the rotation exerts a strong control over the structure of the interior circulation, and they characterized how the characteristic convective velocities vary with heat flux, rotation rate, and radius in the interior. Convection occurs more efficiently near the poles than the equator, leading to equator-to-pole temperature differences that may reach a few K in the interior.

Next, Showman and Kaspi (2013) constructed analytic models for the circulation in the stratified atmosphere. As mentioned previously, the absence of external stellar irradiation impinging on brown dwarfs suggests that they lack the day-night or equator-pole radiative forcing that drives the circulation in the tropospheres of most solar system planets. Showman and Kaspi (2013) suggested that the convection will generate a wide population of atmospheric waves, including both gravity waves (Freytag et al. 2010) and Rossby waves, and that these will propagate vertically and drive a mean flow comprising zonal jets, turbulence, and coherent structures. Figure 25 illustrates the mechanism. Gravity and Rossby waves are generated by convection near the radiative-convective boundary and propagate vertically into the atmosphere, where they break or dissipate, causing a zonal acceleration, which in steady state is balanced by an equal and opposite zonal acceleration due to a meridional flow (blue curves). The ascending and descending motion associated with this circulation advects entropy vertically; because entropy increases with height in a stably stratified atmosphere, this implies that, on an isobar, low entropy—hence cold—air will accompany regions of ascent, whereas high entropy—hence warm—air will accompany regions of descent. Such a circulation is thermally indirect; ascending air is denser than descending air, and so the circulation increases rather than decreases the potential energy of the circulation. In steady state, this potential energy is destroyed by radiation, which causes cooling in the warm regions and warming in the cool regions. Of course, such a circulation cannot occur in isolation; it is driven by the absorption of the waves propagating upward from below, and is therefore referred to as a wave-driven circulation. Such circulations are well known in the solar system community; they occur in the stratospheres of Earth, Jupiter, Saturn, Uranus, and Neptune (e.g. Andrews et al. 1987; Conrath et al. 1990; West et al. 1992; Moreno and Sedano 1997).

Fig. 25
figure 25

Illustration of a large-scale, wave-driven atmospheric circulation, as occurs in the stratospheres of solar system planets (including Earth, Jupiter, and Saturn) and as Showman and Kaspi (2013) proposed also occurs in the stratified atmospheres of brown dwarfs. Vertically propagating gravity and Rossby waves are generated by convection near the radiative-convective boundary and propagate upward into the atmosphere (orange wavy arrows), where they are absorbed or break. Their damping generates a coherent eddy acceleration (black ⊙ symbol, representing a vector coming toward the viewer), which in turn drives a mean wind (black and red ⊙ symbols, respectively, representing vectors coming toward the viewer). In steady state, the eddy acceleration is balanced by a Coriolis acceleration (black ⊗ symbol, representing a vector pointing away from the viewer) associated with a slow secondary circulation (blue curve). The ascending and descending branches of this circulation advect entropy vertically, leading to horizontal temperature contrasts

Showman and Kaspi (2013) constructed an analytic theory for such a circulation with the aim of predicting the characteristic horizontal temperature contrasts and wind speeds in brown dwarf atmospheres. The efficiency with which wave absorption drives the mean flow is difficult to predict from first principles, so Showman and Kaspi (2013) considered it a free parameter, described by \(\eta \), which is a dimensionless number giving the ratio of the power per area used to drive the circulation to the total IR heat flux radiated to space. No fundamental theory for \(\eta \) in planetary stratospheres yet exists, but observations of Earth and Jupiter indicates \(\eta \sim 10^{-3}\). This number is generally yet unknown in the context of brown dwarfs, and here we apply a similar number \(\eta \sim 10^{-3}\) in the following analysis. Showman and Kaspi (2013) showed that, to within numerical factors of order unity, the characteristic horizontal temperature differences \(\Delta T_{\mathrm{horiz}}\), horizontal wind speeds \(\Delta u\), and vertical wind speeds \(w\) are given by

$$\begin{aligned} &\frac{\Delta T_{\mathrm{horiz}}}{T} \sim \eta ^{1/2} \frac{NH}{c_{s}} \end{aligned}$$
(20)
$$\begin{aligned} &\frac{\Delta u }{c_{s}} \sim \eta ^{1/2} \, l L_{D} \end{aligned}$$
(21)
$$\begin{aligned} &w \sim \eta ^{1/2} \frac{H }{\tau _{\mathrm{rad}}} \frac{c_{s} }{NH} \end{aligned}$$
(22)

where, as before, \(NH\) is the characteristic horizontal phase speed of long-wavelength gravity waves, \(c_{s}\) is the speed of sound, \(L_{D}\) is the deformation radius, \(l\) is the horizontal wavenumber of the flow, corresponding to a characteristic horizontal length scale \(L\approx \pi /l\), and \(\tau _{\mathrm{rad}}\) is the radiative time constant. For a vertically isothermal, ideal-gas H2 atmosphere, \(NH/c_{s} \approx 0.4\), which is order unity. The length scale of the flow is unknown, but values several times the deformation radius are plausible, which suggests that \(l L_{D}\) likewise is order unity. For small wave-driving efficiencies (\(\eta \ll 1\)), these equations therefore predict that the fractional temperature differences are small, the wind speeds are much less than the speed of sound, and the timescale for air to advect vertically over a scale height, \(H/w\), is much longer than the radiative time constant. For \(\eta \sim 10^{-3}\) and typical brown dwarf conditions, these equations predict that \(\Delta T_{\mathrm{horiz}}\sim 30~\text{K}\), \(\Delta u\sim 10^{2}~\text{m}\,\text{s}^{-1}\), and \(H/w\sim 10\tau _{\mathrm{rad}}\sim 10^{6}~\text{s}\), corresponding to vertical velocities of order \(10^{-2}~\text{m}\,\text{s}^{-1}\), given a typical brown dwarf scale height of 5 to 10 km. Note that these velocities correspond to vertical velocities in the stably stratified atmosphere associated with the large-scale meridional, wave-driven overturning circulation; convective velocities in the underlying convection zone are of course greater.

To further investigate whether convective perturbations can drive an atmospheric circulation and understand how the flow organizes, Zhang and Showman (2014) presented global simulations of the stratified atmosphere, using the 1.5-layer shallow-water equations. Random, small-scale mass sources and sinks were added to represent convective perturbations at the base of the atmosphere. Radiation was represented using a Newtonian cooling scheme that damps horizontal thickness perturbations, consistent with the expectation that the radiative equilibrium temperature structure is nearly independent of longitude and latitude. Zhang & Showman found that when the convective forcing was strong and/or the radiative damping was weak, the flow spontaneously organized into a banded pattern comprising zonal jets with superposed vortices. But when the convective forcing was weak and/or the radiative damping was strong, the atmospheric turbulence damped before it could self-organize into jets, and the flow instead comprised isotropic turbulence and vortices (Fig. 26). Zonal wind speeds ranged from tens to hundreds of \(\text{m}\,\text{s}^{-1}\), consistent with the analytical predictions by Showman and Kaspi (2013). The transition between the jet-dominated and vortex-dominated regimes is relevant to brown dwarfs; for plausible forcing amplitudes, jets and banding tend to occur when \(\tau_{\mathrm{rad}}\gtrsim 10^{5}\text{--}10^{6}~\text{s}\). This suggests that some brown dwarfs may be dominated by spots, whereas others may be dominated by bands.

Fig. 26
figure 26

Shallow-water simulations of the global circulation on isolated brown dwarfs, from Zhang and Showman (2014). (A) Dynamical behavior as a function of the convective forcing rate (abscissa) and strength of damping (ordinate), here represented by the radiative time constant. Note that short radiative time constant implies strong damping of the flow, and vice versa. Orange circles represent models exhibiting a banded flow with prominent zonal jets, which occurs with strong forcing and/or weak damping. Blue circles depict models dominated by horizontally isotropic turbulence, which occurs with weak forcing and/or strong damping. An example of the turbulent regime is shown as the lower left inset and the jet-dominated regime is shown in the upper right inset. Colors in the insets represent geopotential anomaly (in unit of \(10^{5}~\text{m}^{2}\,\text{s}^{-2}\)). (B) A typical synthetic lightcurve from these simulations, calculated by integrating the height (thickness) field over the Earth-facing hemisphere. The solid blue envelope in the main panel shows the time evolution of the overall variability amplitude, which ranges over almost an order of magnitude on timescales of thousands of hours. The insets show the actual lightcurves over short timescales (four rotation periods long each) at three different intervals. These insets illustrate how the lightcurve shapes vary with time, at some times exhibiting complex multi-peaked structures (left and right insets) and at other times exhibiting nearly sinusoidal behavior (middle inset)

These models have implications for explaining surface patchiness and lightcurve variability. Simulated lightcurves calculated from Zhang and Showman (2014)’s models shows that the lightcurves can exhibit a variety of timescales—including not only the rotational modulation but longer, dynamical timescales as well—and allows for changes qualitatively similar to those seen in observed lightcurves (Fig. 26). The overall amplitude of the variability can vary from nearly zero to \(\sim2\%\), and the lightcurves over rotational timescales can alternate between sinusoidal and multi-peaked shapes, similar to those seen in observations (Fig. 23), depending on the evolving global pattern of the turbulence.

Continuing this line of inquiry, Showman et al. (2019) presented high-resolution, 3D numerical simulations of the atmospheric circulation on brown dwarfs forced by small-scale convective perturbations at the base of the stratified atmosphere. Radiative heating/cooling was represented using a Newtonian cooling scheme that, like in Zhang and Showman (2014) damps thermal perturbations and therefore removes energy from the flow. The convective perturbations trigger a wide population of waves, which propagate upward into the atmosphere and induce a wave-driven circulation through a mechanism similar to that described in Fig. 25. Horizontal temperature variations of tens of K on isobars and zonal winds reaching \(200\text{--}400~\text{m}\,\text{s}^{-1}\) occur (Fig. 27), qualitatively consistent with the predictions of Showman and Kaspi (2013). Showman et al. (2019) found that strong zonal jets are confined to low latitudes when the radiative damping is strong (Fig. 27a), but the jets occur at all latitudes when the radiative damping is weaker (Fig. 27c).

Fig. 27
figure 27

Global models of atmospheric circulation on brown dwarfs forced by small-scale convective perturbations at the base of the atmosphere, with no external stellar irradiation, from Showman et al. (2019). Left column shows temperature (K) and right column shows zonal wind (\(\text{m}\,\text{s}^{-1}\)) in three simulations with different radiative time constants—\(10^{5}\), \(10^{6}\), and \(10^{7}~\text{s}\) in the top, middle, and bottom rows, respectively. Rotation period is 5 hours, typical for brown dwarfs

In these 3D models, the equatorial jet often exhibits an oscillation analogous to well-known stratospheric oscillations on solar system planets. In Earth’s equatorial stratosphere, vertically stacked eastward and westward jets—and associated temperature anomalies—migrate downward over time with a period of approximately two years. This phenomenon, dubbed the Quasi-Biennial Oscillation (QBO), is driven by the upward propagation of atmospheric waves—including equatorially trapped Rossby, Kelvin and inertio-gravity waves—that are absorbed in the stratosphere (Baldwin et al. 2001). Similar oscillations have been observed on Jupiter and Saturn. On Jupiter, the oscillation has a period of four years and is called the Quasi-Quadrennial Oscillation (QQO). The oscillation on Saturn has a period of 15 years (half a Saturn year), and is called the Saturn Semi-Annual Oscillation (see Showman et al. 2018a for a review). These are the first full 3D models of a giant planet to capture this type of oscillation. In the simulations, the oscillation typically has periods of several to \(\sim12~\text{years}\), broadly similar to the periods of the QQO and SAO. These models highlight the likelihood that long-term oscillations may be a generic feature on warm-to-cool giant planets and brown dwarfs. A key implication is that brown dwarfs will probably be variable not only on short timescales, but also on multiannual or multidecadal timescales as well.

There is a subtle difference in the models of Zhang and Showman (2014) and Showman et al. (2019) as to whether jets exist when damping is very strong. In Showman et al. (2019), models of a given wind speed exhibit more coherent zonal banding when the forcing and damping are weak than when they are strong, because in the former case, the flow has ample time to reorganize into a strong series of zonal jets, whereas in the latter case, the forcing and damping continually attempt to disrupt the banding pattern as it forms. However, some degree of zonal banding always occurs—none of the models in Showman et al. exhibit entirely isotropic turbulence, no matter how short the radiative time constant. This occurs because, in a 3D model, the radiation damps the horizontal temperature variations, which via thermal wind suppresses the vertical shear of the zonal wind. However, the barotropic mode—that is, the height-independent component of the zonal flow—is not damped by radiative forcing. This differs from Zhang and Showman (2014) because their 1.5-layer model—which assumes a quiescent abyssal layer—lacks a barotropic mode. In the 3D formulation, not only strong radiative damping but also strong frictional drag at the base of the atmosphere is required to suppress the jets entirely. Such frictional damping could be supplied by Lorentz forces that act to brake the deep roots of the zonal jets that penetrate into the deep interior, but this is a subtle issue that requires further investigation.

The vertical motions that cause the temperature contrasts in these models (Fig. 27) will transport condensate species upward, which can lead to cloud formation if the condensation level occurs in the atmosphere, as expected for L and early T dwarfs. Very small cloud condensate particles (radii less than perhaps 0.1 μm) exhibit very slow gravitational settling, so they may become well mixed by the circulation into quasi-ubiquitous layers of fine particles. Larger particles will settle more quickly after condensation, and since the regions of condensation are confined to ascending air—which covers only a fraction of the globe—this will naturally lead to cloud patchiness on large scales. Tan (2018) confirmed this picture with numerical simulations similar to those described here but that include coupling to passive cloud tracers that advect with the circulation and settle gravitationally but do not influence the dynamics (i.e., no cloud radiative feedback). As expected from the qualitative considerations above, these models produce regionally patchy cloud decks, which can modulate the outgoing IR flux and lead to lightcurve variability.

4.2.2 Atmospheric Circulation II: 1D Cloud-Radiative Feedbacks in an Atmospheric Column

Recent work demonstrates that the feedback between radiation, dynamics, and cloud formation can lead to spontaneous emergence of an atmospheric circulation, and therefore represents a totally distinct means of generating an atmospheric circulation on brown dwarfs than the wave-driven circulation mechanism just described. Here, we survey two distinct feedbacks that can lead to atmospheric circulation on brown dwarfs, each of which can separately lead to a vigorous atmospheric circulation and regional-scale cloud patchiness—even in the absence of convective forcing that might trigger any wave-driven circulation.

The importance of cloud radiative feedbacks on brown dwarfs is easy to appreciate. Imagine an atmospheric vertical column with an opaque cloud deck—perhaps several thousand km wide—surrounded by cloud-free air. In the cloudy region, the IR radiation to space escapes from top of the cloud, at relatively low pressure, whereas in the cloud-free region, the radiation to space occurs from much deeper pressures, perhaps a few bars. The two regions will therefore experience extremely different vertical profiles of the radiative heating/cooling, which will lead to horizontal temperature differences on isobars. These horizontal temperature contrasts will drive an overturning circulation that, in turn, can advect cloud condensate vertically, and in principle might be capable of maintaining the cloud patchiness.

To illustrate, imagine an initially quiescent, cloud-free atmosphere in radiative equilibrium, with an effective temperature of \(T_{\mathrm{eff}} = 1000~\text{K}\), corresponding to an IR flux to space \(F=\sigma T_{\mathrm{eff}}^{4} = 6\times 10^{4}~\text{W}\,\text{m}^{-2}\), which is likely radiated to space from a pressure near \(1~\text{bar}\) (Fig. 22). Because the atmosphere is in radiative equilibrium, this same flux is supplied to the atmosphere from below, such that the next heating of the column is zero. Now for sake of argument suppose that a cloud appears at a higher altitude, where the temperature is cooler—let us say \(900~\text{K}\). Radiation to space from the cloudy region will occur from the top of the cloud deck, corresponding to a smaller flux \(4\times 10^{4}~\text{W}\,\text{m}^{-2}\). The perturbation in outgoing IR flux is \(\Delta F\approx 2\times 10^{4}~\text{W}\,\text{m}^{-2}\). Integrated over an atmospheric column, this differential heating induces a rate of temperature change over time \(\Delta F g/c_{p} p \sim 10^{-2}~\text{K}\,\text{s}^{-1}\), where we have adopted \(g\approx 500~\text{m}\,\text{s}^{-2}\), \(c_{p}=1.3\times 10^{4}~\text{J}\,\text{kg}^{-1}\,\text{K}^{-1}\), and \(p\approx 1~\text{bar}\).Footnote 21 Over only a few brown dwarf rotations, such a heating differential would cause temperature changes \(\gtrsim 100~\text{K}\) within the atmospheric column.

Such thermal changes due to the appearance of clouds cause two distinct processes relevant for driving an atmospheric circulation. First, they cause an adjustment of the thermal structure of the entire atmospheric column, which changes the atmospheric stratification profile and potentially alters the altitude ranges where vertical mixing and clouds occur. Second, if such a cloudy atmospheric column is surrounded by columns of relatively cloud-free air, the horizontal differences in the vertical heating/cooling profile causes horizontal temperature contrasts, and the resulting horizontal pressure-gradient forces can potentially drive an overturning circulation. The key question is whether either of these process is self-generating and self-sustaining, that is, whether an initially quiescent atmosphere will experience spontaneous emergence of atmospheric circulation due to these processes. Recent work indicates that the answer is yes.

Tan and Showman (2019) constructed a 1D, time-dependent model of the vertical structure in a brown dwarf atmosphere to investigate the coupling between clouds, radiation, and mixing in a vertical atmospheric column. This model is similar in spirit to a long history 1D brown dwarf atmosphere models (e.g Burrows et al. 1997, 2006; Marley et al. 2002, 2010; Saumon and Marley 2008), except that most of these models assume steady state and solve for the radiative-convective-equilibrium, whereas Tan and Showman (2019) allowed for time dependence associated with time-varying deviations from radiative equilibrium. They adopted a simple cloud scheme in which clouds settle gravitationally and are vertically mixed in regions where convection should occur. Their models show that the equilibrium state is generally not steady; radiative-cloud feedbacks cause the thermal and cloud structure to vary, with growth and collapse of cloud decks occurring episodically on timescales of \(\sim1\) to \(\sim30~\text{hours}\) (Fig. 28). As foreshadowed above, temperature varies on isobars by \(\sim 100~\text{K}\) as the cloud deck grows and decays in a cyclical fashion. Because of the time-varying cloud opacity, the IR photosphere pressure moves vertically by two scale heights, which causes the effective temperature to vary by hundreds of K over a cycle. Tan and Showman (2019) showed that the existence of variability is robust over a wide range of parameters including number density of cloud condensation nuclei (CCN), gravity, interior entropy, and the implementation of the vertical mixing—although the exact oscillation period and other details depend on these parameter values. A key point is that the oscillatory behavior is entirely spontaneous and does not require any externally imposed perturbations. This contrasts with other 1D models where variability emerged only in response to imposed thermal perturbations (Robinson and Marley 2014).

Fig. 28
figure 28

Results of 1D time-dependent models of brown dwarfs showing the development of cyclical oscillations in the cloud structure and thermal flux radiated to space. From Tan and Showman (2019). (A) Time-mean cloud mixing ratio (\(q_{c}\)) as a function of pressure. (B) Time evolution of cloud mixing ratio as a function of time and pressure in logarithmic scale. The dotted red line is the pressure at which the contribution function is at a maximum, which is where the majority of flux escape from the atmosphere. (C) The time-mean temperature as a function of pressure is represented as the thick solid line; the dashed line represents enstatite (MgSiO3) condensation curve; and the grey lines represent the envelope of the variation range of the temperature-pressure profile. (D) The corresponding outgoing thermal flux as a function of time

The growth and collapse of cloud decks shown in Fig. 28 results from several interacting feedbacks. In the growth phase, a positive feedback promotes upward expansion of the cloud top, as follows: an upwardly perturbed cloud top cools by adiabatic expansion, leading to cooler cloud-top temperatures and smaller IR flux radiating to space from the cloud top. Since the temperature profile above the cloud top is determined by the strength of this IR flux, a lessening of this OLR causes the entire overlying \(T/p\) profile (above the cloud) to cool off. For a given temperature on an isobar near the top of the cloud, a decrease in the overlying temperatures acts to destabilize the air near the cloud top, allowing the cloud to expand upward. But this growth process is eventually halted by a negative feedback that kills the cloud: although the cloud itself remains relatively well-mixed, a large upward heat flux against the base of the cloud generates stable stratification below the cloud, which shuts off the source of moisture (i.e. condensable chemical species) into the cloud. Without a source of condensable material to replenish it, gravitational settling of particles then inevitably depletes the cloud of particulate matter, causing the cloud to dissipate. The cycle can then repeat.

To reiterate, this is a solely 1D feedback that does not require large-scale dynamics. On a brown dwarf, distinct atmospheric vertical columns in different location on the globe cloud in principle independently undergo such oscillatory behavior; the interactions between such columns will not be negligible, a topic requiring a multidimensional dynamical model to investigate. It is also important to emphasis that the 1 to 30 hour period of the oscillation has nothing to do with rotation—it purely results from the interactions between radiation, convective mixing, and cloud physics—yet its similarity to the rotation timescale means that brown dwarf lightcurves can vary over short (rotational) timescales. Note that this 1D feedback occurs more easily when the cloud condensation occurs near the top of the convective zone. This is because convection below the cloud base is needed to replenish the condensable vapor after clouds gravitationally settling out—otherwise without replenishment, the clouds will dissipate away and the system will be eventually inactive. 1D model atmospheres where the cloud deck would form in the stratified atmosphere could potentially lack variability from this mechanism (Tan and Showman 2019).

4.2.3 Atmospheric Circulation III: Multi-D feedbacks Between Clouds, Radiation, and Dynamics

Next, we turn to consider the second cloud-radiative feedback. Imagine an initially quiescent, cloud-free atmosphere in radiative equilibrium, and imagine seeding a small cloud into a vertical air column that is surrounded by cloud-free air. Being at higher altitudes where temperatures are cooler, the IR flux to space from the top of the cloud is less than that of the surrounding, cloud-free regions. If the initial air columns were in radiative equilibrium—receiving the same heat flux from below as they radiate to space—then this decrease in outgoing heat flux caused by the cloud implies that, integrated vertically over an atmospheric column, the cloudy air column experiences net heating. This heating promotes ascent, causing the cloud to grow to higher altitudes. Adiabatic expansion decreases the temperature still further, causing greater net heating of the column and further promoting growth of the cloud. Gierasch et al. (1973) demonstrated this system is linearly unstable, in that small perturbations will grow over time, and he suggested its possible applicability to giant planets. The mechanism is thought to promote convective self-aggregration in Earth’s tropics (e.g., Wing et al. 2017). However, the mechanism has not been previously investigated for giant planets in the fully developed, nonlinear regime.

Tan and Showman (2020a) have shown that this mechanism can cause the spontaneous emergence of a self-generating atmospheric circulation on brown dwarfs, even when the atmosphere would be quiescent in the absence of clouds. Figure 29 shows a numerical simulation from Tan and Showman (2020a) demonstrating that this mechanism indeed constitutes a positive feedback that can generate atmospheric turbulence. Tan and Showman (2020a) solved the primitive equations for the stratified atmosphere on an \(f\)-plane (corresponding to a 3D Cartesian domain where the Coriolis parameter \(f\) is assumed constant) over a vertical domain from 10 to 0.01 bars. The initial, cloud-free atmosphere was set to be in radiative equilibrium. Through the positive feedback just described, a small cloud seeded into the domain (top left panel) leads to net heating of the column, ascending motion, cloud growth, and horizontal divergence of the flow. The Coriolis force deflects the horizontally diverging air, leading to the generation of a warm-core, anticyclonic, geostrophically balanced, cloud-covered vortex with a diameter of \(\sim20{,}000~\text{km}\) on a timescale of tens of hours.

Fig. 29
figure 29

Growth of a cloudy vortex on an \(f\)-plane from an initially cloud-free environment through a positive feedback between radiation, dynamics, and cloud development. Plotted is the cloud mixing ratio (colorscale) and winds (vectors) at 0.4 bars. In this model, \(f=3\times 10^{-4}~\text{s}^{-1}\). From Tan and Showman (2017)

Long-time integrations demonstrate that this process is self-generating and leads to an active atmospheric circulation that reaches a statistical, fluctuating turbulent equilibrium. Figure 30, from Tan and Showman (2020a), show snapshots of the equilibrated state for several different rotation rates. Cloudy and cloud-free vortices, filaments, and other time-variable turbulent structures emerge, with characteristic length scales that are smaller the faster the rotation rate. This trend matches expectations, since the geostrophic adjustment process should produce vortices with sizes close to the deformation radius, which decreases with increasing rotation rate. However, vortex mergers can cause an inverse energy cascade that transfers energy to larger and larger scales, so the characteristic flow length scale in the fully turbulent state (Fig. 30) can potentially become greater than the deformation radius, depending on the dissipation rate of the kinetic energy (implemented as a bottom frictional drag in the models). For typical brown-dwarf parameters, horizontally averaged wind speeds are typically \(300\text{--}500~\text{m}\,\text{s}^{-1}\), and the local, outgoing IR flux to space can vary by nearly a factor of two.

Fig. 30
figure 30

3D simulations of cloudy brown dwarf atmospheres on the \(f\)-plane from simulations performed in Tan and Showman (2020a). Shown are cloud abundance at 0.2 bars (left), temperature at 0.2 bars (middle) and outgoing IR flux to space (right) for three different rotation rates in the top, middle, and bottom rows. Cloud-radiative-dynamical feedbacks lead to an active atmospheric circulation

Global, fully spherical models of this type from Tan (2018) and Tan & Showman (submitted) likewise show the development of complex, time evolving turbulence involving cloudy and cloud-free patches. Figure 31 shows the outgoing IR flux in the fully equilibrated state for models with rotation periods of 2.5, 5, 10, and 20 hours. Again, as expected, the dominant turbulent length scales are smaller when the rotation rate is faster (rotation period is shorter). The equatorial dynamics differ from those at higher latitudes because of the larger local Rossby numbers and the presence of equatorially trapped waves, which couple to the cloud-radiative feedbacks in complex ways. Locally, the outgoing IR flux again varies by a factor of two (Fig. 30), and integrated over the disk, the lightcurve variability on rotational timescales can exceed 1%, consistent with those commonly observed on brown dwarfs (e.g., Metchev et al. 2015).

Fig. 31
figure 31

Global circulation models of isolated brown dwarfs with a dual-band radiative transfer scheme, coupled to a cloud cycle, showing the development of highly time-variable turbulence and cloud patchiness. Shown is the total outgoing IR flux, which varies locally by a factor of two due to local variations in cloud opacity and temperature structure. The characteristic size of the turbulent structures decreases with decreasing rotation period (increasing rotation rate). From Tan (2018) and Tan & Showman (submitted)

4.3 Summary of Brown Dwarf Dynamics

Studies performed to date demonstrate that at least three independent mechanisms can lead to a vigorous circulation in the stratified atmospheres of brown dwarfs, despite the absence of any stellar irradiation that could produce day-night or equator-to-pole thermal contrasts. These mechanisms include:

  • Cloud-free dynamics: Convection impinges on the stably stratified atmosphere, triggering waves and driving turbulence, zonal jets, and vortices. This wave-driven circulation can organize clouds into a patchy structure on large scales, even if cloud radiative feedbacks are not considered.

  • 1D cloud feedback: A cloud-radiative-mixing feedback in 1D vertical columns drives cyclical growth and collapse of cloud decks in an atmospheric column. This causes IR flux variability of hundreds of K in effective temperature over periods of \(\sim1\text{--}30~\text{hours}\).

  • Multi-D cloud feedback: A cloud-radiative-dynamical instability in 3D drives a circulation that maintains the cloud patchiness. Wind speeds reach \(300\text{--}500~\text{m}\,\text{s}^{-1}\), horizontal temperature differences on isobars reaches 100–200 K, and local variations in effective temperature reach hundreds of K, corresponding to local variations in the outgoing IR flux of a factor of two.

We note that the global-scale models for these three processes performed to date have essentially been idealized process studies that have aimed to investigate each of the three processes in isolation. The behavior when the three mechanisms are combined has not yet been simulated in global-scale numerical models but will likely involve even richer behavior than the models shown here that have investigated them one by one. Future modeling efforts beyond the idealized framework summarized here include the use of GCMs with non-grey radiative transfer along with equilibrium or transport-induced disequilibrium chemistry, which are essential for straight model-observation comparisons. Other than clouds, possible transport-induced chemical patchiness and its radiative effect may also contribute to affect the global circulation as well as lightcurve variability, similar to those explored in the context of hot Jupiters (see Sect. 2.6). Utilizing convection resolving models and the inclusion of both the convective and stratified zones would help to understand turbulence and wave generation in the stratified layers and their effects on large-scale flow. Convection models may also be used to examine the relatively small-scale cloud formation and its associated dynamics, as well as effects of chemical reactions in modulating the convective instability in brown dwarf atmospheres proposed by Tremblin et al. (2017a) and Tremblin et al. (2019).

5 Irradiated Brown Dwarfs

Although most close-in, highly irradiated EGPs have masses less than a few \(M_{J}\), and most known brown dwarfs are isolated objects receiving negligible stellar flux, there exists a population of massive brown dwarfs on close-in orbits that receive intense stellar irradiation. The irradiation levels they experience are comparable to typical hot Jupiters, and yet because they are massive, they transport interior heat fluxes that rival or even exceed the energy flux absorbed from their parent star. As such, they populate the upper right corner of Fig. 2. While currently a less well-known subclass than that of hot Jupiters or field brown dwarfs, they are equally important, for they provide the link in the chain of understanding between brown dwarfs and hot Jupiters. One can think of them as hot Jupiters with enormous convective heat fluxes, or as brown dwarfs that are modified by introducing external stellar forcing.

These objects fall into two distinct populations. The first are brown dwarfs orbiting main-sequence stars in close-in orbits. The rarity of these objects relative to less massive, close-in EGPs was discovered early, a phenomenon dubbed the “brown dwarf desert” (e.g., Marcy and Butler 2000; Sahlmann et al. 2011). Nevertheless, about two dozen transiting, close-in brown dwarfs of this type have been discovered. Their orbital periods—and hence rotation periods if synchronously rotating—range typically from 3–5 days, implying that rotation plays a dynamical role similar to that on hot Jupiters. The second population are brown dwarf companions to white dwarf stars, which survived a common envelope phase during which the more massive companion swelled during the red giant phase to envelope the less-massive companion. Friction caused by orbital motion of the companion through the red giant’s tenuous outer layers causes the brown dwarf to slowly spiral in. Once the red-giant phase ends, the brown dwarf is left orbiting the remnant white dwarf on an extremely tight orbit with a period of typically 1–4 hours. Because they are likewise expected to be synchronously rotating, their rotation periods should be equally short, implying that rotation plays a more controlling role in their atmospheric dynamics. Moreover, because the white dwarf’s radiation is primarily in the ultraviolet, where typical atmospheric gases are much more absorbing than in the visible, such brown dwarfs may experience very different atmospheric thermal structures and chemistry than their counterparts orbiting sunlike stars.

Lightcurves and other observations of both populations have been acquired, suggesting that these objects, like hot Jupiters, experience large day-night temperature differences. Examples are presented in Fig. 32, which shows IR lightcurves from KELT-1b (Beatty et al. 2019), a 27-\(M_{J}\) brown dwarf orbiting an F star on a 1.2-day orbit, and WD0137-349B, which is a 53-\(M_{J}\) brown dwarf orbiting its primary white dwarf every 1.9 hours. Like hot Jupiters, the IR lightcurve of KELT-1b reaches a peak flux substantially before secondary eclipse (Fig. 32, left panel), suggesting the possibility of a similar dynamical regime involving equatorial superrotation and global-scale waves driven by the day-night thermal forcing. In contrast, the brown dwarf orbiting a white dwarf lacks such a phase-curve offset (Fig. 32, right panel), suggesting that differing processes operate at fast rotation. Several such brown-dwarf-white-dwarf systems have been characterized, including NLTT 5306 (Steele et al. 2013), WD0137-349 (Casewell et al. 2015; Longstaff et al. 2017), EPIC 21223532 (Casewell et al. 2018a), WD1202-024 (Rappaport et al. 2017) and SDSS J141126.20+200911.1 (Littlefair et al. 2014; Casewell et al. 2018b).

Fig. 32
figure 32

Spitzer IR lightcurves of irradiated brown dwarfs. Left panel shows 3.6 μm (grey) and 4.5 μm phase curves of the \(27~M_{J}\) brown dwarf KELT-1b, orbiting a main-sequence (F5) star in a 1.2-day orbit. The phase curves shows two secondary eclipses (at phase 0 and 1) and one transit (at phase 0.5). Modified from Beatty et al. (2019). Dashed line shows the bottom of the secondary eclipses, which represents the radiation from the star alone. Right panel shows 3.6 μm (black), 4.5 μm (blue), 5.8 μm (green) and 8.0 μm phase curves of WD0137-349, which comprises a \(53~M_{J}\) brown dwarf orbiting a white dwarf in a 1.9-hour orbit. Two full orbits are shown. In both phase curves, 0 and 1 represent the times when the brown dwarf’s dayside is aimed toward Earth, whereas phase 0.5 and 1.5 represent times when its nightside is aimed toward Earth. From Casewell et al. (2015)

Note that, from the perspective of observational characterization, brown dwarfs around main-sequence stars face similar observational challenges as typical hot Jupiters—the planetary flux is typically thousands of times dimmer than the parent star, so teasing out the planetary signature from total signal is challenging. In the case of brown-dwarf-white-dwarf systems, however, the IR flux from the brown dwarf often dominates, aiding observational characterization of the brown dwarf.Footnote 22 Many brown dwarfs of this type are cataclysmic variables (e.g., Hellier 2001), in which the brown dwarf sheds mass onto the white dwarf. In this situation, the brown dwarf’s non-spherical shape and the presence of an accretion disc around the white dwarf can complicate an interpretation of the observations (Santisteban et al. 2016). Still, some observed brown-dwarf-white-dwarf systems—including the one shown in Fig. 32—are “detached,” exhibiting sufficiently great separations for these effects to be minimal.

Motivated by these systems, Tan and Showman (2020b) performed idealized GCM experiments of synchronously rotating, strongly irradiated giant planets with rotation periods ranging from 80 hours to 2.5 hours to characterize the role of fast rotation in shaping their circulation. A subset of these models is shown in Fig. 33. Under identical forcing conditions, faster rotation leads to larger day-night temperature differences in the observable atmosphere (pressures less than a few hundred mbar), which results from the inhibition of day-night heat transport caused by rapid planetary rotation. Moreover, because of the decrease in the Rossby deformation radius at fast rotation, the equatorial waveguide, Matsuno-Gill pattern, and the equatorial superrotating jet all become confined within \(10^{\circ }\) of the equator. As a result, in contrast to hot Jupiters, the overall dayside temperature pattern does not exhibit any prominent eastward offset when the rotation rate is very fast (essentially, the eastward offset still occurs but becomes confined preferentially closer and closer to the equator as rotation rate increases). Infrared lightcurves therefore exhibit flux peaks that coincide with the time of secondary eclipse. This helps to explain the fact that many brown dwarfs on close-in orbits around white dwarfs lack prominent offsets of the IR flux peaks, as shown in Fig. 32. Interestingly, Fig. 33 shows that deeper in the troposphere—at pressures of a few hundred mbar—the dayside temperature develops a “fingering” pattern wherein meridionally narrow, zonally elongated tongues of hot and cold air interleave. Outside the equatorial jet, a prominent pattern of alternating eastward and westward zonal jets emerges. Tan and Showman (2020b) showed that, despite the novel forcing, these zonal jets scale well with the Rhines scale. Further work is needed to understand the implications for brown-dwarf-white-dwarf systems and the connection to zonal jets on Jupiter and Saturn.

Fig. 33
figure 33

Atmospheric circulation in synchronously rotating, strongly irradiated atmospheres at modest to very fast rotation rate, motivated by observations of brown-dwarf-white-dwarf systems with orbital periods of several hours. Top, middle, and bottom rows represent idealized GCM simulations with rotation periods of 40, 10, and \(2.5~\text{hours}\), respectively. The forcing setup is identical in all models; only the rotation period differs. Left and middle columns show temperature (colorscale) and winds (arrows) in snapshots at 7 mbar and 330 mbar, respectively; right column shows the zonal-mean zonal wind (colorscale, in \(\text{m}\,\text{s}^{-1}\)) as a function of latitude and pressure. As rotation period decreases, the day-night temperature difference increases, the equatorial jet narrows, the hemisphere-scale dayside hotspot offset decreases, and numerous off-equatorial zonal jets emerge in the mid-to-high latitudes. From Tan and Showman (2020b)

Lee et al. (2020) simulated the atmospheric circulation of the brown dwarf orbiting the white dwarf WD0137-349 with an orbital period of about 2 hours. They adopted the ExoFMS GCM which utilizes a dual-grey radiative transfer scheme assuming constant opacity in both thermal and visible band. Their results are qualitatively similar to the rapidly rotating cases in Tan and Showman (2020b), showing a meridionally narrower equatorial superrotating jet compared to canonical hot Jupiter simulations, nearly geostrophic flows at mid-high latitudes and large day-night temperature difference. Although using different GCMs and radiative forcing setup, the agreement between Lee et al. (2020) and Tan and Showman (2020b) suggests that strong rotation could robustly and profoundly influence the atmospheric circulation of tidally locked objects.

Future endeavors beyond the above idealized work are needed to obtain a completed understanding on the 3D atmospheres of this group of objects. The spectrum of white dwarfs peak at UV, and significant photochemistry are expected to occur in atmospheres of the companion BDs (Longstaff et al. 2017; Casewell et al. 2018b). Effects of photochemical products and radiative heating on the thermal structures and circulation, along with the effects of hydrogen thermal dissociation/recombination similar to those found in ultra-hot Jupiter atmospheres (Bell and Cowan 2018; Komacek and Tan 2018; Tan and Komacek 2019), are worth to explored in a thorough manner. Given the insufficient day-night heat transport, various mineral clouds are expected to condense in the nightsides of these objects (Lee et al. 2020), and may have interesting effects on lightcurve and dynamics. Strong magnetic field (\(\sim \text{kG}\)) has been inferred from observations of field brown dwarfs (e.g., Kao et al. 2016, 2018). If brown dwarfs around white dwarfs host such strong fields, their strongly irradiated atmospheres may experience significant MHD effects.

6 Discussion

Observations reveal a wealth of behavior in the atmospheres of EGPs and brown dwarfs. The intense daynight radiative forcing on hot Jupiters leads to large day-night temperature differences and planetary-scale waves that, according to most GCMs, drive a fast superrotating jet at the equator. Adjustment of the thermal structure by planetary-scale waves provides a framework for understanding the day-night temperature contrasts and their dependence on the strength of stellar irradiation and other parameters. Analytic theories and idealized numerical models have provided significant insights into these processes. Overall, a basic picture of the dynamical processes that maintain the circulation on hot Jupiters is slowly emerging. But many aspects remain poorly understood; for example, the eddies that maintain the superrotating jet are nonlinear and interact with the shear of the jet, making their structure difficult to predict in detail. The ways the circulation interacts with patchy clouds and the detailed dynamics of vertical and horizontal mixing of trace constituents remain less well understood. Predictions for various regime transitions across the hot Jupiter population—and between hot and warm Jupiters—have yet to be confronted by observational data. An open frontier that has been little studied is the extent to which (and how) the atmosphere and interior of hot Jupiters interact, including the depth to which the equatorial jet extends, the type of winds that may exist in the deep convection zone, as well as the degree of mixing, the nature of the thermal structure, and the modes of energy exchange between the atmosphere and deep interior.

Warm Jupiters are a fascinating target for future observations, and help to bridge the gap between hot Jupiters on the one hand and Jupiter and Saturn on the other. At present, few observational data exist to constrain the dynamical behavior of warm Jupiters, but a growing body of models suggest a wide range of behavior that organizes into a variety of regimes depending on the rotation rate, obliquity, orbital eccentricity, and how the photosphere-level radiative time constant compares to the rotational and orbital timescales.

A wealth of data provide information about the thermal and cloud structure of isolated brown dwarfs, indicating that these worlds are covered by patchy clouds modulated by a vigorous atmospheric circulation. This regime differs from hot and warm Jupiters in being fundamentally driven by interior heat loss. In comparison to hot Jupiters, dynamical models of the global circulation on brown dwarfs are in their infancy and represents a ripe frontier for future research. Models point toward rich interactions between the convective interior and stratified atmosphere and of various subtle feedbacks between dynamics, radiation and clouds that can lead to a spontaneous, self-generating circulation comprising turbulence, vortices, and zonal jets, associated with patchy cloud structures like those observed. The continuum between hot brown dwarf on the one hand and Jupiter and Saturn on the other deserves further study, and can face observational tests from the growing quality of data on cooler T and Y dwarfs. Irradiated brown dwarfs represent an even more open frontier, and can help bridge the gap between hot Jupiters and isolated brown dwarfs. Observations are starting to clarify the thermal structure of these objects, but so far almost no modeling and theory has been performed to explore this regime.

We have argued that there is value in considering these distinct populations together as a class. The synergies that can result from comparisons between the diverse populations have yet to be well exploited but can potentially lead to major insights across these populations. The grand challenge of understanding the behavior of giant planets as an integrated population across many orders of magnitude in internal heat flux, external irradiation, rotation rate, and other parameters beckons us.