4 Atmospheric circulation: Possible regimes and influence on evolution

4.1 Basic parameter regime

Rotation plays a central role in the atmospheric dynamics of Pegasi planets, and HD 209458b in particular. The ratio of nonlinear advective accelerations to Coriolis accelerations in the horizontal momentum equation is u f^-1 L^-1 (called the Rossby number), where u is the mean horizontal wind speed, $f=2\Omega \sin\phi$ is the "Coriolis parameter'' (e.g. Holton 1992, pp. 39-40), $\Omega$is the rotational angular velocity, $\phi$ is latitude, and Lis a characteristic length scale. Rossby numbers of 0.03-0.3 are expected for winds of planetary-scale and speeds ranging from $100{-}1000\rm\, m\rm\, s^{-1}$. In Sect. 3 we showed that modest asynchronous rotation may occur, in which case the Rossby number could differ from this estimate by a factor up to $\sim$2. These estimates suggest that nonlinear advective terms are small compared to the Coriolis accelerations, which must then balance with the pressure-gradient accelerations. (The Rossby number could be $\sim$1 if the winds reach $2000{-}3000\rm\, m\rm\, s^{-1}$, which we show below is probably the maximum allowable wind speed.)

   

Table 1: Rhines' and deformation lengths for giant planets.

	u	R	$\Omega$	$L_{\beta}$	$L_{\rm D}$	Jet width
	$(\rm\, m\rm\, s^{-1})$	( $10^7\rm\, m$)	( $10^{-4}\rm\, s^{-1}$)	( $10^7\rm\, m$)	( $10^7\rm\, m$)	( $10^7\rm\, m$)
Jupiter	50	7.1	1.74	1.0	0.2	$\sim$1
Saturn	200	6.0	1.6	1.9	0.2	$\sim$2
Uranus	300	2.6	1.0	2.0	0.2	$\sim$2
Netpune	300	2.5	1.1	1.9	0.2	$\sim$2
HD 209458b	?	10	$\sim$0.2	15(u/1 km s -1)1/2	$\sim$4	?

Note. $L_{\rm D}$ calculated at the tropopause using $N=0.01\rm\, s^{-1}$for Saturn, Uranus, and Neptune, and $N=0.02\rm\, s^{-1}$ for Jupiter.

The zonality of the flow can be characterized by the Rhines' wavenumber, $k_{\beta} \sim (\beta/u)^{1/2}$, where $\beta$ is the derivative of f with northward distance (Rhines 1975). The half-wavelength implied by this wavenumber, called the Rhines' scale $L_{\beta}$, provides a reasonable estimate for the jet widths on all four outer planets in our solar system (Cho & Polvani 1996; see Table 1). For HD 209458b, the Rhines' scale is ${\sim}1.5\times10^{10} (u/1000\rm\, m\rm\, s^{-1})^{1/2} \rm\,cm$, which exceeds the planetary radius if the wind speed exceeds about $400\rm\, m\rm\, s^{-1}$.

Another measure of horizontal structure is the Rossby deformation radius (Gill 1982, p. 205), $L_{\rm D} \sim N H/f$, where H is the scale height and N is the Brünt-Vaisala frequency (i.e., the oscillation frequency for a vertically displaced air parcel; Holton 1992, p. 54). At pressures of a few bars, the temperature profiles calculated for irradiated extrasolar giant planets by Goukenlouque et al. (2000) suggest $N\sim0.0015\rm\, s^{-1}$. With a scale height of 700 km, the resulting deformation radius is 40000 km. In contrast, the deformation radii near the tropopause of Jupiter, Saturn, Uranus, and Neptune are of order 2000 km (Table 1).

The estimated Rhines' scale (for winds of ${>}500\rm\, m\rm\, s^{-1}$, which we show later are plausible speeds) and deformation radius of Pegasi planets$\,\,$are similar to the planetary radius, and they are a larger fraction of the planetary radius than is the case for Jupiter, Saturn, Uranus, and Neptune (Table 1). This fact suggests that eddies may grow to hemispheric scale in the atmospheres of Pegasi planets$\,\,$and that, compared with the giant planets in our solar system, the general circulation hot Jupiters may be more global in character. Unless the winds are extremely weak, Pegasi planets$\,\,$are unlikely to have >10 jets as do Jupiter and Saturn.

An upper limit on the atmospheric wind speed can be derived from shear-instability considerations. We assume that no zonal winds are present (u(P₀)=0) in the convective core, a consequence of synchronization by tidal friction. The build-up of winds at higher altitudes in the radiative envelope is suppressed by Kelvin-Helmholtz instabilities if the shear becomes too large. This occurs when the Richardson number becomes smaller than 1/4 (cf. Chandrasekhar 1961), i.e. when

$$Ri={N^2\over ({\rm d}u/{\rm d}z)^2} < {1\over 4}\cdot$$ (7)

In the perfect gas approximation, assuming a uniform composition,

$$N^2={g\over H}(\nabla_{\rm ad}-\nabla_T),$$ (8)

where ${\cal R}$ is the universal gas constant divided by the mean molecular mass, $c_{\rm p}$ is the specific heat, $H={\cal R}T/g$ is the pressure scale height, $\nabla_{\rm ad}={\cal R}/c_{\rm p}$ is the adiabatic gradient, $\nabla_T={\rm d}\ln T/{\rm d}\ln p$, T is temperature, and p is pressure. The hypothesis of uniform composition is adequate despite hydrogen dissociation because the timescales involved are generally much longer than the dissociation timescales. It can be noted that a shear instability could appear at larger Richardson numbers in the presence of efficient radiative diffusion (Zahn 1992; Maeder 1995). This possibility will not be examined here.

The maximal wind speed at which Kelvin-Helmholtz instabilities occur can then be derived by integration of Eq. (7):

$$u_{\rm max}(P)\sim {1\over 2}\int_{P_0}^{P} \left[{\cal R}T(\nabla_{\rm ad}-\nabla_T)\right]^{1/2} {\rm d}\ln p.$$ (9)

The value of $u_{\rm max}$ thus derived at $P\sim 1$ bar is of the order of $3000\rm\,m\,s^{-1}$, to be compared to the winds of Jupiter, Saturn, Uranus and Neptune which reach $100{-}500\rm\,m\,s^{-1}$ (Ingersoll et al. 1995). Our estimate for Pegasi planets, assuming the winds are measured in the synchronously-rotating frame, may be uncertain by a factor of $\sim$2 due to the possibility of a non-synchronously-rotating interior. In comparison, the expected speed of sound at the tropopause of HD 209458b is $\sim$ $2400\rm\, m\rm\, s^{-1}$.

A characteristic timescale for zonal winds to redistribute temperature variations over scales similar to the planetary radius Rthen stems from $\tau_{\rm zonal}\mathrel{\hbox{\hbox to 0pt{ \lower.6ex\hbox{$\sim$ }\hss}\raise.4ex\hbox{$>$ }}}R/u_{\rm max}$.

The radiative heating timescale can be estimated by a ratio between the thermal energy within a given layer and the layer's net radiated flux. In the absence of dynamics, absorbed solar fluxes balance the radiated flux, but dynamics perturbs the temperature profile away from radiative equilibrium. Suppose the radiative equilibrium temperature at a particular location is $T_{\rm rad}$ and the actual temperature is $T_{\rm rad} + \Delta T$. At levels close to optical depth unity, the net flux radiated towards outer space is then $4\sigma T_{\rm rad}^3 \Delta T$ and the radiative timescale is

$$\tau_{\rm rad}\sim \frac{P}{g}\frac{c_{\rm p}}{4\sigma T^3},$$ (10)

where $\sigma$ is the Stefan-Boltzmann constant.

Figure 4 shows estimates of $\tau _{\rm zonal}$ and $\tau _{\rm rad}$ for HD 209458b calculated using the temperature profiles from the "hot'' (thin line) and "cold'' (thick grey line) models from Sect. 2. The zonal timescale is estimated by calculating the maximum wind speed that can exist as a function of pressure given the static stability associated with each model, while the radiative time is calculated using the temperatures and heat capacities shown in Fig. 1. At pressures exceeding 0.1bar, radiation is slower than the maximal advection by zonal winds, but by less than one order of magnitude. The consequent day/night temperature difference $\Delta T_{\rm day-night}$ to be expected is:

$${\Delta T_{\rm day-night}\over \Delta T_{\rm rad}}\sim 1 - {\rm e}^{-\tau_{\rm zonal}/\tau_{\rm rad}},$$ (11)

where $\Delta T_{\rm rad}$ is the day-night difference in radiative equilibrium temperatures. Rough estimates from Fig. 4 suggest that $\tau_{\rm zonal}/\tau_{\rm rad} \sim0.3$ at 1 bar, implying that $\Delta T_{\rm day-night}/\Delta T_{\rm rad}\sim$0.3. If $\Delta T_{\rm rad}=1000\rm\,K$, this would imply day-night temperature differences of 300K at 1 bar. Values of $\Delta T_{\rm day-night}$ even closer to $\Delta T_{\rm rad}$ are likely given the fact that slower winds will lead to an even more effective cooling on the night side and heating on the day side.

The small radiative time scale implies that, for the day-night temperature difference to be negligible near the planet's photosphere, atmospheric winds would have to be larger than the maximum winds for the onset of shear instabilities.

Figure 4: Left: characteristic time scales as a function of pressure level. $\tau _{\rm zonal}$ is the minimal horizontal advection time (dashed). $\tau _{\rm rad}$ is the timescale necessary to cool/heat a layer of pressure P and temperature T by radiation alone assuming the gas is optically thick (solid). (The plotted radiative timescale is an underestimate away from optical depth unity, i.e., at pressures less than about 0.3 bars or greater than about 3 bars.) For each case, the thin black and thick grey lines correspond to the hot and cold models from Sect. 2 and Fig. 1. Right: approximate cooling/heating rate as a function of pressure (see Eq. (18)).

4.2 Nature of the circulation

An understanding of the horizontal temperature difference and mean wind speed is desirable. It is furthermore of interest to clarify the possible geometries the flow may take. For Pegasi planets, we envision that the dominant forcing is the large-scale day-night heating contrast, with minimal role for moist convection. This situation differs from that of Jupiter, where differential escape of the intrinsic heat flux offsets the solar heating contrast and moist convection plays a key role. For Pegasi planets, such an offset between intrinsic and stellar fluxes cannot occur, because the intrinsic flux is 10⁴times less than the total flux. The fact that the day-night heating contrast occurs at hemispheric scale - and that the relevant dynamical length scales for Pegasi planets$\,\,$are also hemispheric (Sect. 4.1) - increases our confidence that simple analysis, focusing on the hemispheric-scale circulation, can provide insight.

Because the Rossby number is small, the dominant balance in the horizontal momentum equation at mid-latitudes is between the Coriolis force and the pressure-gradient force (geostrophic balance). We adopt the primitive equations, which are the standard set of large-scale dynamical equations in a stably-stratified planetary atmosphere. When differentiated with pressure (which is used here as a vertical coordinate), this balance leads to the well-known thermal wind equation (e.g. Holton 1992, p. 75):

$$f{\partial {\vec v}\over\partial\ln p}= -{\cal R}{\vec k}\times \nabla_{\rm H} {\vec T}$$ (12)

where ${\vec v}$ is the horizontal wind, p is pressure, $\nabla_{\rm H}$ is the horizontal gradient, T is temperature, and ${\vec k}$ is the unit upward vector. Assuming the interior winds are small compared to the atmospheric winds, this implies that, to order-of-magnitude,

$$\vert{\vec v}\vert \sim {{\cal R}\over f R}\Delta T_{\rm horiz}\Delta\ln p$$ (13)

where $\vert{\vec v}\vert$ and $\Delta T_{\rm horiz}$ are the characteristic wind speed and horizontal temperature difference in the atmosphere and $\Delta\ln p$ is the difference in log-pressures from bottom to top of the atmosphere. Non-synchronous rotation of the interior would alter the relation by changing the value of f, but the uncertainty from this source is probably a factor of two or less (see Eq. (6) and subsequent discussion).

The thermodynamic energy equation is, using pressure as a vertical coordinate (e.g. Holton 1992, p. 60),

$${\partial T\over \partial t} + {\vec v}\cdot\nabla_{\rm H} T - \omega {H^2 N^2\over {\cal R}p} = {q\over c_{\rm p}}$$ (14)

where t is time, $\omega = {\rm d}p/{\rm d}t$ is vertical velocity, q is the specific heating rate (erg ${\rm\,g}^{-1}{\rm\,s}^{-1}$), and we have assumed an ideal gas.

A priori, it is unclear whether the radiative heating and cooling is dominantly balanced by horizontal advection (second term on left of Eq. (14)) or vertical advection (third term on the left). To illustrate the possibilities, we consider two endpoint scenarios corresponding to the two limits. In the first scenario, the radiation is balanced purely by horizontal advection: zonal winds transport heat from dayside to nightside, and meridional winds transport heat from equator to pole. In the second scenario, the radiation is balanced by vertical advection (ascent on dayside, descent on nightside).

Consider the first scenario, where horizontal advection dominates. Generally, we expect ${\vec v}$ and $\nabla_{\rm H} T$to point in different directions, and to order of magnitude their dot product equals the product of their magnitudes. (If the direction of $\nabla_{\rm H} T$ is independent of height and no deep barotropic flow exists, then at mid-latitudes one could argue that winds and horizontal pressure gradients are perpendicular to order Ro. However, the existence of either asynchronous rotation or variations in the orientations of $\nabla_{\rm H} T$ with height will imply that winds and $\nabla_{\rm H} T$ are not perpendicular even if $Ro\ll 1$. Because asynchronous rotation and variation in the directions of $\nabla_{\rm H} T$ with height are likely, and because, in any case, ${\vec v}$ and $\nabla_{\rm H} T$ will not be perpendicular near the equator, we write the dot product as the product of the magnitudes.) An order-of-magnitude form of the energy equation is then

$${\vert{\vec v}\vert \Delta T_{\rm horiz}\over R}\sim {q\over c_{\rm p}}\cdot$$ (15)

The solutions are

$$\vert{\vec v}\vert\sim\left({q\over c_{\rm p}} {{\cal R}\Delta\ln p \over f}\right)^{1/2}$$ (16)

$$\Delta T_{\rm horiz}\sim R\left({q\over c_{\rm p}} {f\over {\cal R}\Delta\ln p}\right)^{1/2}\cdot$$ (17)

A rough estimate of the heating rate, $q/c_{\rm p}$, results from the analysis in Sect. 4.1:

$${q\over c_{\rm p}} = {4\sigma T^3 \Delta T g\over p c_{\rm p}}$$ (18)

where $\Delta T$ is the characteristic magnitude of the difference between the actual and radiative equilibrium temperatures in the atmosphere. The heating rate depends on the dynamics through $\Delta T$. We simply evaluate the heating rate using Eq. (18) with $\Delta T \approx T/2$. The results are shown in Fig. 4.

The key difficulty in applying the equations to Pegasi planets$\,\,$is the fact that $q/c_{\rm p}$ depends on pressure, and $\Delta T_{\rm horiz}$ probably should too, but Eqs. (16) and (17) were derived assuming that $\Delta T_{\rm horiz}$ is constant. We can still obtain rough estimates by inserting values of $q/c_{\rm p}$at several pressures. At $\sim$50-100 bars, where $q/c_{\rm p}\sim 10^{-4}{-}10^{-3}\rm\,K\rm\, s^{-1}$, we obtain temperature differences and wind speeds of 50-$150\rm\,K$ and 200- $600\rm\, m\rm\, s^{-1}$. At 1 bar, where $q/c_{\rm p}$ reaches $10^{-2}\rm\,K\rm\, s^{-1}$, the estimated temperature contrast and wind speed is $\sim$$500\rm\,K$ and $\sim$ $2000\rm\, m\rm\, s^{-1}$. The estimates all assume $f\approx3\times10^{-5}\rm\, s^{-1}$, ${\cal R}=3500\rm\, J\rm\,kg\rm\,K^{-1}$, and $\Delta\ln p\approx 3$.

Later we show that the equations successfully predict the mean wind speeds and temperature differences obtained in numerical simulations of the circulation of HD 209458b. This gives us confidence in the results.

Now consider the second scenario, where vertical advection (third term on the left of Eq. (14)) balances the radiative heating and cooling. The magnitude of $\omega$ can be estimated from the continuity equation. Purely geostrophic flow has zero horizontal divergence, so $\omega$ is roughly

$$\omega \sim Ro {\vert{\vec v}\vert \over R}p.$$ (19)

Strictly speaking, this is an upper limit, because the nonlinear terms comprising the numerator of the Rossby number contain some terms (e.g., the centripetal acceleration) that are divergence-free. Substituting this expression into the energy equation, using Eq. (13), and setting $Ro = \vert{\vec v}\vert/f R$ implies that

$$\vert{\vec v}\vert\sim {R\over N H} \left( {q\over c_{\rm p}} f {\cal R}\right)^{1/2}$$ (20)

$$\Delta T_{\rm horiz}\sim {f R^2\over N H \Delta\ln p}\left({q\over c_{\rm p}} {f\over {\cal R}}\right)^{1/2}\cdot$$ (21)

Numerical estimates using H=700 km, $N\sim0.0015\rm\, s^{-1}$, and $q/c_{\rm p}=10^{-4}{-}10^{-3}\rm\,K\rm\, s^{-1}$ (appropriate to 50-100 bars) yield temperature differences and speeds of $80{-}250\rm\,K$ and $300{-}900\rm\, m\rm\, s^{-1}$, respectively. Here we have used the same values for R and f as before. A heating rate of $10^{-2}\rm\,K\rm\, s^{-1}$, appropriate at the 1 bar level, yields $\vert{\vec v}\vert\sim3000\rm\, m\rm\, s^{-1}$ and $\Delta T_{\rm horiz}\sim 800\rm\,K$. These values are similar to those obtained when we balanced heating solely against horizontal advection.

We can estimate the vertical velocity under the assumption that vertical advection balances the heating. Expressing the vertical velocity w as the time derivative of an air parcel's altitude, and using $w=-\omega/\rho g$, where $\rho$ is density, yields

$$w\sim{q\over c_{\rm p}}{{\cal R}\over H N^2}\cdot$$ (22)

Using $q/c_{\rm p}=10^{-2}\rm\,K\rm\, s^{-1}$ and H=700 km implies that $w\sim 20\rm\, m\rm\, s^{-1}$ near 1 bar. If this motion comprises the vertical branch of an overturning circulation that extends vertically over a scale height and horizontally over $\sim$10¹⁰ cm (a planetary radius), the implied horizontal speed required to satisfy continuity is $\sim$ $3000\rm\, m\rm\, s^{-1}$, consistent with earlier estimates.

The numerical estimates, while rough, suggest that winds could approach the upper limit of $\sim$ $3000\rm\, m\rm\, s^{-1}$implied by the shear instability criterion. This comparison suggests that Kelvin-Helmholtz shear instabilities may play an important role in the dynamics.

The likelihood of strong temperature contrasts can also be seen from energetic considerations. The differential stellar heating produces available potential energy (i.e., potential energy that can be converted to kinetic energy through rearrangement of the fluid; Peixoto & Oort 1992, pp. 365-370), and in steady-state, this potential energy must be converted to kinetic energy at the rate it is produced. This requires pressure gradients, which only exist in the presence of lateral thermal gradients. A crude estimate suggests that the rate of change of the difference in gravitational potential $\Phi$ between the dayside and nightside on isobars caused by the heating is $\sim$ ${\cal R}(q/c_{\rm p}) \Delta\ln p$. The rate per mass at which potential energy is converted to kinetic energy by pressure-gradient work is ${\vec v}\cdot\nabla_{\rm H}\Phi$, which is approximately $\vert{\vec v}\vert {\cal R}\Delta T_{\rm horiz}\Delta\ln p/R$. Equating the two expressions suggests, crudely, that $\Delta T_{\rm horiz}\sim (q/c_{\rm p})(R/\vert{\vec v}\vert)$, which implies $\Delta T_{\rm horiz}\sim 500\rm\,K$ using the values of $q/c_{\rm p}$ and $\vert{\vec v}\vert$near 1 bar discussed earlier.

One could wonder whether a possible flow geometry consists of dayside upwelling, nightside downwelling, and simple acceleration of the flow from dayside to nightside elsewhere (leading to a flow symmetrical about the subsolar point). The low-Rossby-number considerations described above argue against this scenario. A major component of the flow must be perpendicular to horizontal pressure gradients. Direct acceleration of wind from dayside to nightside is still possible near the equator, however. This would essentially be an equatorially-confined Walker-type circulation that, as described earlier, would lie equatorward of the geostrophic flows that exist at more poleward latitudes.

4.3 Numerical simulations of the circulation

To better constrain the nature of the circulation, we performed preliminary three-dimensional, fully-nonlinear numerical simulations of the atmospheric circulation of HD 209458b. For the calculations, we used the Explicit Planetary Isentropic Coordinate, or EPIC, model (Dowling et al. 1998). The model solves the primitive equations in spherical geometry using finite-difference methods and isentropic vertical coordinates. The equations are valid in stably-stratified atmospheres, and we solved the equations within the radiative layer from 0.01 to 100 bars assuming the planet's interior is in synchronous rotation with the 3.5-day orbital period. The radius, surface gravity, and rotation rate of HD 209458b were used ( $10^8\rm\, m$, $10\rm\, m\rm\, s^{-2}$, and $\Omega = 2.1\times10^{-5}\rm\, s^{-1}$, respectively).

The intense insolation was parameterized with a simple Newtonian heating scheme, which relaxes the temperature toward an assumed radiative-equilibrium temperature profile. The chosen radiative-equilibrium temperature profile was hottest at the substellar point ($0^{\circ }$ latitude, $0^{\circ }$ longitude) and decreased toward the nightside. At the substellar point, the profile's height-dependence was isothermal at 550 K at pressures less than 0.03 bars and had constant Brunt-Vaisala frequency of $0.003\rm\, s^{-1}$ at pressures exceeding 0.03 bars, implying that the temperature increased with depth. The nightside radiative-equilibrium temperature profile was equal to the substellar profile minus $100\rm\,K$; it is this day-night difference that drives all the dynamics in the simulation. The nightside profile was constant across the nightside, and the dayside profile varied as

$$T_{\rm dayside} = T_{\rm nightside} + \Delta T_{\rm rad} \cos\alpha$$ (23)

where $T_{\rm dayside}$ is the dayside radiative-equilibrium temperature at a given latitude and longitude, $T_{\rm nightside}$ is the nightside radiative-equilibrium temperature, $\Delta T_{\rm rad}=100\rm\,K$, and $\alpha$ is the angle between local vertical and the line-of-sight to the star. For simplicity, the timescale over which the temperature relaxes to the radiative-equilibrium temperature was assumed constant with depth with a value of $3\times10^5\rm\, s$. This is equal to the expected radiative timescale at a pressure of about 5 bars (Fig. 4, left).

The Newtonian heating scheme described above is, of course, a simplification. The day-night difference in radiative-equilibrium temperature, 100 K, is smaller than the expected value. Furthermore, the radiative timescale is too long in the upper troposphere ($\sim$0.1-1bar) and too short at deeper pressures of $\sim10{-}100$ bars. Nevertheless, the net column-integrated heating per area produced by the scheme ( $\rm\,W\rm\, m^{-2}$) is similar to that expected to occur in the deep troposphere of HD 209458b and other Pegasi planets, and the simulation provides insight into the circulation patterns that can be expected.

The temperature used as the initial condition was isothermal at 500 K at pressures less than 0.03 bars and had a constant Brunt-Vaisala frequency of $0.003\rm\, s^{-1}$ at pressures exceeding 0.03 bars. There were no initial winds. The simulations were performed with a horizontal resolution of $64\times32$ with 10 layers evenly spaced in log-pressure.

Figures 5-8 show the results of such a simulation. Despite the motionless initial condition, winds rapidly develop in response to the day-night heating contrast, reaching an approximate steady state after $\sim$400 Earth days. Snapshots at 42 and 466 days are shown in Figs. 5 and 6, respectively. In these figures, each panel shows pressure on an isentrope (greyscale) and winds (vectors) for three of the model layers corresponding to mean pressures of roughly 0.4, 6, and 100 bars (top to bottom, respectively). The greyscale is such that, on an isobar, light regions are hot and dark regions are cold.

Figure 5: EPIC simulation of atmospheric circulation on HD 209458b after 42 days. Panels depict pressure on an isentrope (greyscale) and winds (vectors) for three model levels with mean pressures near 0.4, 6, and 100 bars (top to bottom, respectively). From top to bottom, the maximum wind speeds are 937, 688, and $224\rm\, m\rm\, s^{-1}$, respectively, and the greyscales span (from dark to white) 0.31-0.49 bars, 5.6-7.8 bars, and 92-113 bars, respectively. Substellar point is at $0^{\circ }$ latitude, $0^{\circ }$ longitude.

Figure 6: Simulation results for HD 209458b at 466 days (after a steady state has been reached). As in Fig. 5, panels depict pressure on isentropes (greyscale) and winds (vectors) for three model levels with mean pressures near 0.4, 6, and 100 bars (top to bottom, respectively). From top to bottom, the maximum wind speeds are 1541, 1223, and $598\rm\, m\rm\, s^{-1}$, respectively, and the greyscales span (from dark to white) 0.32-0.51 bars, 5.6-8.1 bars, and 90-127 bars, respectively. Substellar point is at $0^{\circ }$ latitude, $0^{\circ }$ longitude.

Figure 7: Additional simulation results at 466 days. Arrows are identical to those in Fig. 6, but here greyscale is vertical velocity ${\rm d}\theta /{\rm d}t$, where $\theta$ is potential temperature. Light regions are heating (i.e., $\theta$ is increasing) and dark regions are cooling ($\theta$ is decreasing). From top to bottom, the greyscales span (from dark to white) -0.003 to $0.003\rm\,K\rm\, s^{-1}$, -0.002 to $0.002\rm\,K\rm\, s^{-1}$, and -0.001 to $0.001\rm\,K\rm\, s^{-1}$, respectively. The lowermost panel separates the convective interior from the radiative region, and the implication is that mass exchange can happen across this interface.

The simulation exhibits several interesting features. First, peak winds exceed 1 km s^-1, but despite these winds, a horizontal temperature contrast is maintained. The dayside (longitudes $-90^{\circ}$ to $90^\circ$) is on average hotter than the nightside, but dynamics distorts the temperature pattern in a complicated manner. Second, an equatorial jet develops that contains most of the kinetic energy. The jet initially exhibits both eastward and westward branches (Fig. 5) but eventually becomes only eastward (Fig. 6), and extends from $-30^{\circ}$ to $30^{\circ}$in latitude. Third, away from the equator, winds develop that tend to skirt parallel to the temperature contours. This is an indication that geostrophic balance holds. (Because the thermal structure is independent of height to zeroeth order and we have assumed no deep barotropic flow, horizontal pressure and temperature gradients are parallel.) Nevertheless, winds are able to cross isotherms near the equator, and this is important in setting the day-night temperature contrast. As expected from the order-of-magnitude arguments in Sect. 4.1, the jets and gyres that exist are broad in scale, with a characteristic width of the planetary radius.

In the simulation, the day-night temperature difference (measured on isobars) is about $50\rm\,K$. This value depends on the adopted heating rate, which was chosen to be appropriate to the region where the pressure is tens of bars. Simulations that accurately predict the day-night temperature difference at 1 bar will require a more realistic heating-rate scheme; we will present such simulations in a future paper.

The temperature patterns in Figs. 5 and 6 show that Earth-based infrared measurements can shed light on the circulation of Pegasi planets. In Fig. 6, the superrotating equatorial jet blows the high-temperature region downwind. The highest-temperature region is thus not at the substellar point but lies eastward by about $60^{\circ}$ in longitude. The maximum and minimum temperatures would thus face Earth before the transit of the planet behind and in front of the star, respectively. On the other hand, if a broad westward jet existed instead, the maximum and minimum temperatures would face Earth after the transits. Therefore, an infrared lightcurve of the planet throughout its orbital cycle would help determine the direction and strength of the atmospheric winds.

In the simulation, the intense heating and cooling causes air to change entropy and leads to vertical motion, as shown in Fig. 7. Light regions (Fig. 7) indicate heating, which causes ascent, and dark regions indicate cooling, which causes descent. Air is thus exchanged between model layers. This exchange also occurs across the model's lowermost isentrope, which separates the radiative layer from the convective interior (Fig. 7, bottom). Because upgoing and downgoing air generally have different kinetic energies and momenta, energy and momentum can thus be exchanged between the atmosphere and interior.

Figure 8 indicates how the steady state is achieved. The build-up of the mass-weighted mean speed (top panel) involves two timescales. Over the first 10 days, the mean winds accelerate to $110\rm\, m\rm\, s^{-1}$, and the peak speeds approach 1 km s^-1. This is the timescale for pressure gradients to accelerate the winds and force balances to be established. The flow then undergoes an additional, slower ($\sim$300 day) increase in mean speed to $300\rm\, m\rm\, s^{-1}$, with peak speeds of $\sim$1.5 km s^-1. This timescale is that for exchange of momentum with the interior, across the model's lowermost isentrope, to reach steady state, as shown in Fig. 8, middle.

Figure 8: Top: mass-weighted mean wind speed within the model domain ($\sim$0.01-100 bar) during the simulation of HD 209458b shown in Figs. 5-7. Middle: globally-averaged zonal angular momentum per area minus that of the synchronously-rotating state. (That is, ${RA}^{-1} \int u\cos\phi\,{\rm d}m$, where A is the planet's area, u is eastward speed, $\phi$ is latitude, and dm is an atmospheric mass element.) Bottom: globally-averaged flux of kinetic energy across the model's bottom isentrope (at $\sim$100 bars), which is the interface between the radiative layer and the convective interior in this model.

In the simulation, kinetic energy is transported from the atmosphere into the interior at a rate that reaches $2500\rm\,W\rm\, m^{-2}$(Fig. 8, bottom). This downward transport of kinetic energy is $\sim$1% of the absorbed stellar flux and is great enough to affect the radius of HD 209458b, as shown in Paper I. Although the simulation described here does not determine the kinetic energy's fate once it reaches the convective interior, we expect that tidal friction, Kelvin-Helmholtz instabilities, or other processes could convert these winds to thermal energy. The exact kinetic energy flux will depend on whether a deep barotropic flow exists. Furthermore, potential and thermal energy are also transported through the boundary, and pressure work is done across it. Our aim here is not to present detailed diagnostics of the energetics, but simply to point out that energy fluxes that are large enough to be important can occur. We are currently conducting more detailed simulations to determine the sensitivity of the simulations to a deep barotropic flow, and we will present the detailed energetics of these simulations in a future paper.

The evolution of angular momentum (Fig. 8, middle) helps explain why the circulation changes from Figs. 5 to 6. After 42 days (Fig. 5), eastward midlatitude winds have already developed that balance the negative equator-to-pole temperature gradient. But at this time, the atmosphere's angular momentum is still nearly zero relative to the synchronously-rotating state, so angular momentum balance requires westward winds along some parts of the equator. (Interestingly, in these regions, the temperature increases with latitude, as expected from geostrophy.) The circulation involves upwelling on the dayside and equatorial flow - both east and west - to the nightside, where downwelling occurs. This circulation is similar in some respects to the Walker circulation in Earth's atmosphere. After 466 days, however, the atmosphere has gained enough angular momentum that the equatorial jet is fully superrotating. Some westward winds exist at high latitudes, but they are weak enough that the net angular momentum is still eastward.

Application of the arguments from Sect. 4.2 provides a consistency check. The mean heating rate in the simulations is about $1.7\times10^{-4}\rm\,K\rm\, s^{-1}$. Inserting this value into Eqs. (16)-(17), we obtain $\Delta T_{\rm horiz}= 70\rm\,K$and $\vert{\vec v}\vert = 240\rm\, m\rm\, s^{-1}$. These compare well with the day-night temperature difference and mass-weighted mean speed of $50\rm\,K$ and $300\rm\, m\rm\, s^{-1}$ that are obtained in the simulation. (For the estimate, we used $f=3\times10^{-5}\rm\, s^{-1}$, ${\cal R}= 3500\rm\, J\rm\,kg^{-1}\rm\,K^{-1}$, $R=10^8\rm\, m$, and $\Delta p=3$.) If Eqs. (20)-(21) are adopted instead, using values of $N=0.003\rm\, s^{-1}$ and H=300 km, we obtain temperature differences and speeds of $130\rm\,K$ and $470\rm\, m\rm\, s^{-1}$. Compared to the simulation, these estimates are too high by a factor of $\sim$2, but are still of the correct order of magnitude. An important point is that the maximum speed in the simulation substantially exceeds that from the simple estimates; the estimates are most relevant for the mean speed.

4.4 Vertical motion and clouds

As shown in the appendix, hot and cold regions on the planet may have distinct chemical equilibrium compositions. Advection of air between these regions plays a key role for the (disequilibrium) chemistry in the atmosphere. Cloud formation critically depends on whether vertical motions dominate over horizontal motions. This in turns affects the albedo and depth to which stellar light is absorbed.

Clouds form when the temperature of moving air parcels decreases enough for the vapor pressure of condensible vapors to become supersaturated. If dayside heating and nightside cooling are balanced only by horizontal advection, air flowing from dayside to nightside remains at constant pressure but decreases in temperature (hence entropy) because of the radiative cooling. Nightside condensation may occur, and particle fallout limits the total mass of the trace constituent (vapor plus condensates) in this air. As the air flows onto the dayside, its temperature increases and any existing clouds will sublimate. This scenario therefore implies that the dayside will be cloud-free.

In an alternate scenario, dayside heating and nightside cooling are balanced by vertical advection, at least near the sub- and antisolar points. Because entropy increases with altitude in a stably-stratified atmosphere, this scenario implies that ascent occurs on the dayside and descent on the nightside. Clouds may therefore form on the dayside, increasing the albedo and decreasing the pressure at which stellar light is absorbed.

In reality, both horizontal and vertical advection are important, and the real issue is to determine the relative contribution. We are pursuing more detailed numerical simulations to address this issue.