Atmospheric Circulation of Hot Jupiters: Dayside–Nightside Temperature Differences. II. Comparison with Observations

2.1.1. Theory

In Komacek & Showman (2016), we derived approximate steady-state analytic solutions for dayside–nightside temperature differences relative to those in radiative equilibrium. This theory was derived using the weak temperature gradient limit, with vertical entropy advection globally balancing radiative heating/cooling (parameterized by a Newtonian cooling scheme), as Komacek & Showman (2016) showed that the vertical entropy advection term is larger (in a globally averaged sense) than the horizontal entropy advection term at photospheric pressure levels ($\gtrsim 10\ \mathrm{mbar}$) in hot Jupiter atmospheres. The solutions in Komacek & Showman (2016) were written for all possible balances in the momentum equation, where the day–night pressure-gradient force is balanced near the equator by either advection or drag and at high latitudes by either the Coriolis force, advection, or drag.

Applying the theory of Komacek & Showman (2016) to explain trends in day–night temperature differences from their GCMs with varying mean molecular weight and heat capacity, Zhang & Showman (2016) derived a uniform expression for the dayside–nightside temperature difference that incorporates all dynamical terms (see their Appendix). In this work, we use their complete solution for dayside–nightside temperature differences in order to compare to infrared observations of hot Jupiter phase curves. Their solution for approximate dayside–nightside temperature differences (relative to that in radiative equilibrium) is

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{\Delta }}T}{{\rm{\Delta }}{T}_{\mathrm{eq}}}\sim 1-\displaystyle \frac{2}{\alpha +\sqrt{{\alpha }^{2}+4{\gamma }^{2}}},\end{eqnarray} \tag{ 1 }$$

where

$$\begin{eqnarray}&&\alpha =1+\displaystyle \frac{({\rm{\Omega }}+{\tau }_{\mathrm{drag}}^{-1}){\tau }_{\mathrm{wave}}^{2}}{{\tau }_{\mathrm{rad}}{\rm{\Delta }}\mathrm{ln}p}\end{eqnarray} \tag{ 2 }$$

and

$$\begin{eqnarray}&&\gamma =\displaystyle \frac{{\tau }_{\mathrm{wave}}^{2}}{{\tau }_{\mathrm{rad}}{\tau }_{\mathrm{adv},\mathrm{eq}}{\rm{\Delta }}\mathrm{ln}p}.\end{eqnarray} \tag{ 3 }$$

In Equations (1)–(3) above, all variables retain their meaning from Komacek & Showman (2016) and Zhang & Showman (2016). ¹

For the comparison with observed results below, we calculate dayside–nightside temperatures at a pressure of $100\ \mathrm{mbar}$, similar to the level of the infrared photosphere for a typical hot Jupiter. We take a planetary radius $a=9.43\,\times {10}^{7}\ {\rm{m}}$, specific gas constant $R=3700\,{\rm{J}}\,{\mathrm{kg}}^{-1}\,{{\rm{K}}}^{-1}$, specific heat ${c}_{p}=1.3\times {10}^{4}\,{\rm{J}}\,{\mathrm{kg}}^{-1}\,{{\rm{K}}}^{-1}$, and gravity $g=9.36\ {\rm{m}}\ {{\rm{s}}}^{-2}$. We assume that in radiative equilibrium the nightside temperature is small relative to the dayside temperature, and hence that ${\rm{\Delta }}{T}_{\mathrm{eq}}={T}_{\mathrm{eq}}$. Additionally, we scale ${\tau }_{\mathrm{rad}}$ from Komacek & Showman (2016) using the power-law relationship with temperature from Showman & Guillot (2002) and Ginzburg & Sari (2016), setting

$$\begin{eqnarray}&&{\tau }_{\mathrm{rad}}(p,{T}_{\mathrm{eq}})={10}^{5}\ {\rm{s}}\left(\displaystyle \frac{p}{100\ \mathrm{mbar}}\right){\left(\displaystyle \frac{1800{\rm{K}}}{{T}_{\mathrm{eq}}}\right)}^{3}.\end{eqnarray} \tag{ 4 }$$

2.1.2. Comparison with Observations

Figure 2 shows theoretical predictions for the fractional dayside–nightside temperature difference A with a fixed rotation period of $3.5\ \mathrm{days}$ and ${\tau }_{\mathrm{drag}}$ varying in the range of ${10}^{3}-\infty \ {\rm{s}}$, compared to the set of observed fractional dayside–nightside brightness temperature differences ${A}_{\mathrm{obs}}$. Here we take a long rotation period to effectively bracket the lower limit of expected dayside–nightside temperature differences, as a decreased rotation period increases the resulting dayside–nightside temperature differences (see Figure 3 and related discussion below). The theoretical predictions at low ${\tau }_{\mathrm{drag}}$ match well the trend of steadily increasing ${A}_{\mathrm{obs}}$ with increasing equilibrium temperature. However, note that the predictions for ${\tau }_{\mathrm{drag}}={10}^{5}-\infty \ {\rm{s}}$ regularly underpredict the individual and band-averaged values of ${A}_{\mathrm{obs}}$. ²

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

There are three ways to increase dayside–nightside temperature differences relative to the prediction for the case with ${\tau }_{\mathrm{drag}}=\infty$. First, as seen in Figure 2, atmospheric drag can be strong, with a characteristic timescale shorter than ${10}^{5}\ {\rm{s}}$. In the context of Lorentz "drag," the variations from planet to planet could be caused by changing values of internal magnetic field strength or varying atmospheric composition and hence thermal ionization fraction. Second, we are here considering a relatively long rotation period (relevant for HD 209458b), and considering shorter rotation periods will cause the predicted value of A to increase. This is because if the Coriolis force is stronger (due to a larger Ω), it can support a larger day-to-night pressure gradient (which is related to the day–night temperature contrast) in steady state. This can be seen in Figure 3, which shows how the theoretical prediction of dayside–nightside temperature differences varies with the rotation period and incident stellar flux, assuming no atmospheric drag. Though incident stellar flux is the dominant factor controlling the dayside–nightside temperature contrast, planets with very short rotation periods (where the Coriolis force is stronger than the advective terms) have significantly larger day–night temperature differences. We will include the effects of rotation rate to directly compare with observations in Section 2.2. Lastly, as will be discussed further in Section 2.2, any process that shifts the optical-depth-unity surface to lower pressures can explain an increase in A. One example of this is supersolar metallicities, which have been shown to increase the day–night temperature difference in GCMs (e.g., Showman et al. 2009; Kataria et al. 2015; Wong et al. 2016). Another possibility is high-altitude clouds, which are expected to form largely on the nightside and western limb (Lee et al. 2016; Oreshenko et al. 2016; Parmentier et al. 2016), where it is coldest. As a result, clouds could increase the amplitude of the phase curve and thereby raise our theoretical predictions for A closer to observed values.

Using the theory from Komacek & Showman (2016) and Zhang & Showman (2016), we can make predictions of A for individual planets, assuming that they are tidally locked and hence that their orbital period is equal to their rotation period and choosing a drag strength. In Figure 4, we plot theoretical predictions of A for all planets from Figure 1 along with their corresponding observations, setting ${\tau }_{\mathrm{drag}}=\infty$. We choose ${\tau }_{\mathrm{drag}}=\infty$ for these comparisons because this gives lower limits on the resulting predictions for fractional dayside–nightside temperature differences. The theoretical predictions at $100\ \mathrm{mbar}$ (near the expected infrared photosphere of a typical hot Jupiter) match well with ${A}_{\mathrm{obs}}$ in the high equilibrium temperature regime (${T}_{\mathrm{eq}}\gtrsim 2000\ {\rm{K}}$) and underpredict the values for all planets with ${T}_{\mathrm{eq}}\lesssim 2000\ {\rm{K}}$. However, there is a strong dependence of the theoretical predictions for A on pressure, with fractional day–night temperature differences near unity at pressures $\lesssim 10\ \mathrm{mbar}$ and near zero at pressures $\gtrsim 1\ \mathrm{bar}$. As a result, changing photosphere pressures between 10 and 100 mbar in our theory brackets the range of fractional dayside–nightside temperature differences for all observed planets. Note that the short rotation period of $\approx 0.81\ \mathrm{days}$ for WASP-43b does lead to a prediction of significantly larger day–night temperature contrast than for a slower-rotating planet (e.g., HD 209458b). However, this increase in predicted day–night temperature difference due to the short rotation period is not by itself enough to provide a full quantitative explanation of the observed day–night temperature contrast of WASP-43b.

Zoom In Zoom Out Reset image size

Given that we are using a value of Ω consistent with the expected rotation rate of each planet, there are two possible explanations for the underprediction of observed dayside–nightside temperatures by our theory assuming a photosphere pressure of $100\ \mathrm{mbar}$ at low ${T}_{\mathrm{eq}}$: strong drag affecting the ability of the circulation to erase dayside–nightside temperature differences, or photosphere pressures being $\lesssim 100\ \mathrm{mbar}$. For the case of WASP-43b, the photosphere pressure must be lower than $25\ \mathrm{mbar}$ to explain the phase curve amplitude without invoking strong atmospheric drag. Given that Lorentz forces are expected to be weak at low ${T}_{\mathrm{eq}}\lesssim 1400\ {\rm{K}}$ (Rogers & Komacek 2014), it is unlikely that magnetic effects alone can explain the discrepancy between theoretical predictions and observations at low ${T}_{\mathrm{eq}}$.

A decreased photosphere pressure would be related to the presence of the large column abundance of infrared absorbers, which in this case could be cloud coverage muting the outgoing infrared flux. Nightside clouds would naturally tend to be more prevalent in the atmospheres of cooler hot Jupiters than those in hotter hot Jupiters. This is simply because cooler temperatures allow a greater variety of chemical species to condense. This could help explain why we are underpredicting the value of the dayside–nightside temperature difference primarily at the lowest equilibrium temperatures, while matching the observed dayside–nightside temperature differences at higher equilibrium temperatures. Specifically, clouds are expected to be present with an effective cloud coverage near unity on the nightside regardless of ${T}_{\mathrm{eq}}$, but the aerosol particle size is expected to increase with ${T}_{\mathrm{eq}}$ (Spiegel et al. 2009; Heng & Demory 2013; Parmentier et al. 2016). Using the difference in transit radii between the line center and wing of potassium and sodium lines, Heng (2016) showed that hot Jupiters with ${T}_{\mathrm{eq}}\lesssim 1300\ {\rm{K}}$ are cloudier than hotter planets. Additionally, observations have shown that clouds are nearly ubiquitous for cool planets with ${T}_{\mathrm{eq}}\lesssim 2000\ {\rm{K}}$ (Demory et al. 2013; Esteves et al. 2015; Sing et al. 2016). Another possibility is supersolar metallicities, which increase the opacity of the atmosphere and thereby adjust the photosphere upward. However, it is unclear why this effect would no longer be at work in the high-${T}_{\mathrm{eq}}$ regime.

That the theory matches best with observations at high ${T}_{\mathrm{eq}}$ is expected from Komacek & Showman (2016), due largely to the very short radiative timescales at these high values of incident stellar flux. As a result, cloud formation is no longer an effective way to increase the phase curve amplitude in this regime. Additionally, although Lorentz forces are very strong in the high-temperature regime, due to the steep dependence of thermal ionization on temperature (Perna et al. 2010; Rauscher & Menou 2013; Rogers & Komacek 2014), the radiative timescale becomes so short that the observable atmosphere of the planet (at and above the photosphere) is close to radiative equilibrium. This forces large day–night temperature differences at the photosphere (Komacek & Showman 2016), as seen for our high-${T}_{\mathrm{eq}}$ predictions in Figure 4.

To understand the relevant dynamics governing the trend of increasing dayside–nightside temperature differences with increasing incident stellar flux found in Section 2 in greater detail, we turn to numerical (GCM) simulations. As in Komacek & Showman (2016), we use the MITgcm (Adcroft et al. 2004) to solve the hydrostatic primitive equations (Equations (1)–(5) in Komacek & Showman 2016). All aspects of the model dynamics and the drag parameterization, Shapiro filter, and initial-temperature–pressure profile are the same as in Komacek & Showman (2016). However, instead of prescribing heating/cooling rates using a Newtonian cooling formalism as in the theory of Komacek & Showman (2016), we couple a double-gray radiative transfer scheme, similar to that of Heng et al. (2011a) and Rauscher & Menou (2012), to the MITgcm dynamical core. This scheme enables us to both test the qualitative predictions of Komacek & Showman (2016) using a more realistic model and produce simulated phase curves for comparison with observations.

The plane-parallel, two-stream approximation of radiative transfer equations for the diffuse, azimuthally averaged intensity I is as follows (e.g., Liou 2002, Chapter 4.6), considering absorption only and omitting scattering:

$$\begin{eqnarray}&&{\mu }_{1}\displaystyle \frac{{{dI}}^{+}}{d\tau }={I}^{+}-B(\tau ),\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&-{\mu }_{1}\displaystyle \frac{{{dI}}^{-}}{d\tau }={I}^{-}-B(\tau ),\end{eqnarray} \tag{ 6 }$$

where ${I}^{+}$ is the upward intensity and ${I}^{-}$ is the downward intensity of diffuse radiation in the infrared, τ is infrared optical depth, ${\mu }_{1}$ is the mean infrared zenith angle associated with a semi-isotropic hemisphere of radiation, and $B(\tau )$ is the Planck function at the local temperature $T(\tau )$. The double-gray assumption utilizes a separation between short-wavelength (visible) insolation and long-wavelength (thermal) emission. Equations (5) and (6) are solved to give the diffusive thermal flux ${F}_{\mathrm{th}}^{\pm }=\pi {I}_{\mathrm{th}}^{\pm }$, using boundary conditions of zero downward thermal flux at the top of the atmosphere and a prescribed net upward thermal flux at the bottom of our domain. We use the hemi-isotropic (or hemispheric) closure for the thermal band, which assumes that the thermal radiation is isotropic in each hemisphere (upwelling and downwelling). This is a reasonable choice among the other possible closures used for the two-stream approximation, as it ensures energy conservation (Pierrehumbert 2010).

We consider the stellar insolation as only a direct beam source without diffuse components, consistent with our assumption of no scattering. In this limit, the visible flux ${F}_{{\rm{v}}}^{\pm }$ is simply

$$\begin{eqnarray}&&{F}_{{\rm{v}}}^{-}={\mu }_{{\rm{v}}}{F}_{0}\exp \left(-\displaystyle \frac{{\tau }_{{\rm{v}}}}{{\mu }_{{\rm{v}}}}\right),\quad \quad {F}_{{\rm{v}}}^{+}=0,\end{eqnarray} \tag{ 7 }$$

where F₀ is the visible flux incident on the top of the atmosphere at a given longitude and latitude, ${\mu }_{{\rm{v}}}$ is the local zenith angle of the insolation, and ${\tau }_{{\rm{v}}}$ is its optical depth at a given vertical location in the atmosphere.

We utilize the reliable and efficient numerical package TWOSTR (Kylling et al. 1995) to solve the radiative transfer equations for a plane-parallel atmosphere in the two-stream approximation. TWOSTR is based on the general-purpose multistream discrete ordinate algorithm DISORT (Stamnes et al. 1988) and incorporates all the advanced features of that well-tested and stable algorithm. One notable feature of TWOSTR is an efficient treatment of internal thermal sources that may be allowed to vary either slowly or rapidly with depth (Kylling & Stamnes 1992). The applicability of TWOSTR to atmospheric calculations has been tested thoroughly in Kylling et al. (1995). We repeated those tests with our version of the algorithm and performed our own tests to ensure the correct implementation in the MITgcm (see the Appendix for an example).

To couple TWOSTR to the MITgcm, the temperature–pressure structure of each vertical atmospheric column is fed from MITgcm to calculate fluxes with TWOSTR. Then, we calculate the thermodynamic heating rate per unit mass, q, by taking the divergence in pressure coordinates of the net vertical flux F as in Showman et al. (2009),

$$\begin{eqnarray}&&q=g\displaystyle \frac{\partial F}{\partial p},\end{eqnarray} \tag{ 8 }$$

where $F={F}_{{\rm{v}}}^{+}-{F}_{{\rm{v}}}^{-}+{F}_{\mathrm{th}}^{+}-{F}_{\mathrm{th}}^{-}$. This heating rate is then used for the next iteration of the dynamical equations.

In all of our numerical simulations varying incident stellar flux, drag strength, and rotation rate, we use visible and infrared opacities that do not change with these parameters. We do so to better understand the effects that incident stellar flux, drag strength, and rotation rate themselves have on the atmospheric circulation of hot Jupiters. This numerical exploration improves on previous work exploring how the atmospheric circulation varies with drag strength and day–night forcing amplitude using one-layer models (Perez-Becker & Showman 2013) and 3D models with simplified Newtonian cooling (Komacek & Showman 2016). However, this model is not as complex as explorations in this parameter space with the full SPARC/MITgcm model (Showman et al. 2009; Kataria et al. 2015, 2016), which use state-of-the-art radiative transfer with accurate opacities. As a result, this model is a middle rung in a "modeling hierarchy" exploring how the atmospheric circulation, especially the day–night temperature contrast, of hot Jupiters varies with incident stellar flux, drag strength, composition, and rotation rate.

We fix the visible absorption coefficient ${\kappa }_{{\rm{v}}}=4\,\times {10}^{-4}\,{{\rm{m}}}^{2}\,{\mathrm{kg}}^{-1}$ as in Rauscher & Menou (2012). This is the same visible absorption coefficient used in the analytic solutions of Guillot (2010), which is chosen to obtain a good match to the models of Fortney et al. (2008) and Showman et al. (2008). Note that an increased value of the visible opacity (relevant to absorption by possible molecular TiO or VO) would cause a strong thermal inversion in the upper atmosphere, strongly affecting the temperature–pressure profile. We do not consider such an enhanced visible opacity in this work. We set the infrared absorption coefficient to be a power law in pressure

$$\begin{eqnarray}&&{\kappa }_{\mathrm{th}}=c{\left(\displaystyle \frac{p}{1\mathrm{Pa}}\right)}^{b},\end{eqnarray} \tag{ 9 }$$

where $c=2.28\times {10}^{-6}\,{{\rm{m}}}^{2}\,{\mathrm{kg}}^{-1}$ and b = 0.53. This power-law opacity, which is commonly used to account for the effects of collision-induced or pressure-broadened absorption (e.g., Arras & Bildsten 2006; Youdin & Mitchell 2010; Heng et al. 2012; Rauscher & Menou 2012), gives the best fit to the pressure–temperature profile from the analytic solutions of Parmentier & Guillot (2014) and Parmentier et al. (2015) for HD 209458b using the full opacity–pressure–temperature relationship from Freedman et al. (2008).

We take the same planetary parameters for our numerical simulations as Komacek & Showman (2016); see also Section 2.1.1. We systematically vary the incident stellar flux in a sequence of models to determine its effect on the thermal structure and atmospheric circulation. The incident stellar flux can be described using an "irradiation temperature," defined as ${T}_{\mathrm{irr}}={({F}_{\star }/\sigma )}^{1/4}$, where F_⋆ is the incident stellar flux at the top of the atmosphere and σ is the Stefan–Boltzmann constant. Thus, ${T}_{\mathrm{irr}}$ can be thought of as the equilibrium temperature of a blackbody at the substellar point. Equivalently, we can represent the stellar flux through the global-mean equilibrium temperature ${T}_{\mathrm{eq}}$. This is the temperature of a spherical blackbody planet where the heat is efficiently redistributed over all $4\pi$ sr, such that ${T}_{\mathrm{eq}}={T}_{\mathrm{irr}}/{4}^{1/4}$. Note that in these simulations we assume zero planetary albedo such that ${T}_{\mathrm{eq}}$ is an effective control parameter. This is generally consistent with the observed low bond albedos for hot Jupiters (Schwartz & Cowan 2015). We set an internal temperature of ${T}_{\mathrm{int}}=100\ {\rm{K}}$, which is equivalent to an upward thermal flux at the bottom of the domain of ${F}_{\mathrm{th}}^{+}=\sigma {T}_{\mathrm{int}}^{4}$. This prescribed internal flux is small enough to not greatly affect our numerical solutions above the photosphere but large enough to aid in the equilibration of deep pressure levels of our GCM.

In our numerical simulations, we use a horizontal resolution of C32 (approximately equal to a global resolution of 128 × 64 in longitude and latitude). There are 40 vertical levels, with the bottom 39 levels spaced evenly in log pressure between 0.2 and 200 mbar, and a top layer extending from $0.2\ \mathrm{mbar}$ to space. We use a standard fourth-order Shapiro filter, which smooths grid-scale variations and enables the model to maintain numerical stability. We integrate our models until the circulation at photospheric pressures ($\lesssim 1\ \mathrm{bar}$) reaches statistical equilibrium in thermal and kinetic energy. As a result, our simulations are each integrated to $2000\ \mathrm{Earth}\ \ \mathrm{days}$ of model time, and our results are time-averaged over the last $100\ \mathrm{days}$ of model time.

Our main suite of numerical simulations examine how the atmospheric circulation changes with varying incident stellar flux and drag strengths (parameterized by a drag timescale ${\tau }_{\mathrm{drag}}$). This is comparable to our nonlinear suite of simulations varying radiative and drag timescales in Komacek & Showman (2016), except here we solve the actual radiative transfer equations (albeit in a double-gray two-stream approximation) rather than using Newtonian cooling. We use the same drag scheme as in Komacek & Showman (2016), which consists of two parts. First, we include a weak frictional drag at the bottom of the model, extending from $200\ \mathrm{bars}$ to $10\ \mathrm{bars}$. The drag coefficient for this component is largest (while still being relatively weak) at the bottom and decreases with decreasing pressure, becoming zero (i.e., with no applied drag) at and above $10\ \mathrm{bars}$. Second, we optionally add a spatially constant drag term with a drag timescale of ${\tau }_{\mathrm{drag}}$. The total drag coefficient at each pressure level is taken to be the greater of the individual drag coefficients from these two components (see Komacek & Showman 2016, Equation (12)).

As in Komacek & Showman (2016), we examine simulations with varying ${\tau }_{\mathrm{drag}}$ in the range of ${10}^{3}-\infty \ {\rm{s}}$. Simulations with ${\tau }_{\mathrm{drag}}=\infty$ have no large-scale drag at pressures less than $10\ \mathrm{bars}$, but they (and all other simulations presented in this paper) still have the basal drag scheme at the bottom of the domain. We also vary the incident stellar flux at the substellar point (with values of $1.418\times {10}^{4},2.268\times {10}^{5},\,$ $1.148\times {10}^{6},3.628\times {10}^{6},\,$ $8.859\times {10}^{6},\,$ $1.837\times {10}^{7}\,{\mathrm{Wm}}^{-2}$), corresponding to varying global-mean equilibrium temperatures ${T}_{\mathrm{eq}}$ from 500 to 3000 K in steps of $500\ {\rm{K}}$. This results in a grid of 36 simulations, with resulting temperature and wind maps at $80\ \mathrm{mbar}$ shown in Figure 5.

Zoom In Zoom Out Reset image size

The results of this suite of simulations agree well with both theoretical expectations and the simulations using Newtonian cooling in Komacek & Showman (2016); compare Figure 5 here to their Figure 4. Dayside–nightside temperature differences are small at low ${T}_{\mathrm{eq}}$ and become larger with increasing ${T}_{\mathrm{eq}}$, as expected from the theory in Section 2. Additionally, when drag is strong, with a drag timescale ${\tau }_{\mathrm{drag}}\lesssim {10}^{4}\ {\rm{s}}$, characteristic dayside–nightside temperature differences are larger than in simulations with weaker drag. This is expected from the theory in Komacek & Showman (2016), as drag should increase the day–night temperature contrast if the condition ${\tau }_{\mathrm{drag}}\ll {{\rm{\Omega }}}^{-1}$ is satisfied. Given that in these simulations we use a rotation period of $3.5\ \mathrm{days}$ and hence ${{\rm{\Omega }}}^{-1}=4.8\times {10}^{4}\ {\rm{s}}$, we expect that a transition between low and high day–night temperature differences occurs between ${\tau }_{\mathrm{drag}}={10}^{5}\ {\rm{s}}$ and ${\tau }_{\mathrm{drag}}={10}^{4}\ {\rm{s}}$, and this prediction agrees with our current grid of simulations.

As in Komacek & Showman (2016), our simulations show a change between long and short ${\tau }_{\mathrm{drag}}$ in the general character of the atmospheric circulation near the infrared photosphere. At long ${\tau }_{\mathrm{drag}}\gt {10}^{5}\ {\rm{s}}$, there is a dominant eastward jet at the equator. At short ${\tau }_{\mathrm{drag}}$, however, the flow is predominantly from day to night, rather than east–west. This transition may be observationally detectable (Parmentier et al. 2013; Showman et al. 2013a), and indeed the Doppler analysis of HARPS spectra by Louden & Wheatley (2015) found that the atmospheric flow of HD 189733b is likely dominated by an east–west jet, rather than day–night flow.

A subset of our simulations with ${\tau }_{\mathrm{drag}}={10}^{6}\ {\rm{s}}$ and high ${T}_{\mathrm{eq}}$ do not show mirror symmetry of the wind pattern about the equator. These simulations are time-variable, with the eastward jet notably fluctuating in latitude at the location west of the substellar point, where a hydraulic jump has been seen in previous simulations (e.g., Showman et al. 2009). Recent work by Fromang et al. (2016) has shown that shear instabilities can cause hot Jupiter atmospheres to have (potentially observable) time variability, as predicted by Showman & Guillot (2002), Goodman (2009), and Li & Goodman (2010). The GCMs of Showman et al. (2009) predicted on the order of $1 \%$ variability in the secondary eclipse depths of HD 189733b, in agreement with the upper limit on variability found from the observations of Agol et al. (2010). We plan to quantify the observability of time variability in future work, as it has the potential to affect observations of exoplanet atmospheres with JWST.

We post-process the resulting temperature–pressure profiles from our GCMs using TWOSTR in order to calculate light curves as a function of orbital phase. Specifically, we calculate the mean flux that the Earth-facing hemisphere radiates toward Earth, which is a function of orbital phase. Here we write this hemispheric-mean flux radiated toward Earth as ${S}_{\mathrm{av}}(\delta )$, where δ is the orbital phase (taken to be equivalent to the sub-Earth longitude). We assume a circular orbit and zero obliquity, causing the orbital phase to vary linearly with time through the orbit. Here we use the same definition for orbital phase as in Zhang & Showman (2016), such that δ is equal to 0^◦ at primary eclipse and 180^◦ at secondary eclipse. Similarly, $\delta =90^\circ$ and $\delta =270^\circ$ correspond to times in the orbit where the eastern and western terminators of the planet lie at the sub-Earth longitude, respectively. Additionally, note that, as in Zhang & Showman (2016), we define the planetary longitude (here written as λ) such that $\lambda =0^\circ$ at the substellar point and $\lambda =180^\circ$ at the antistellar point.

To calculate the local outgoing (top-of-atmosphere) flux in the direction of Earth at each grid point and at each orbital phase (here written $S(\lambda ,\phi ,\delta )$), we use the local temperature–pressure profile as a function of latitude (here written as ϕ) and longitude from the same GCM results as shown in Figure 5 that have been time-averaged over the last 100 days of integration. We then use TWOSTR to calculate the outgoing flux using the particular value of μ relevant for the latitude and longitude of each grid point, along with taking into account the orbital phase. Here $\mu (\lambda ,\phi ,\delta )$ (similar to the visibility in Cowan et al. 2009) is the cosine of the angle between the local normal to each grid point and the line of sight toward Earth, taking into account both the angle between the longitude of each grid point and the sub-Earth longitude and the angle between the latitude of each grid point and the equator. As a result, this post-processing approach correctly accounts for the directionality of the emitted radiation from the planet.

We then integrate up the outgoing flux radiated toward Earth at each grid point to calculate a hemispheric-mean flux radiated toward Earth. This approach is the same as that used in previous work (e.g., Fortney et al. 2006; Showman et al. 2008, 2009; Cowan et al. 2009; Lewis et al. 2010; Kataria et al. 2015; Showman et al. 2015), computing the integral

$$\begin{eqnarray}&&{S}_{\mathrm{av}}(\delta )=\displaystyle \frac{1}{\pi {a}^{2}}{\int }_{-\pi /2}^{+\pi /2}{\int }_{\pi /2-\delta }^{3\pi /2-\delta }S(\lambda ,\phi ,\delta )\mu (\lambda ,\phi ,\delta ){dA},\end{eqnarray} \tag{ 10 }$$

where ${dA}={a}^{2}\cos (\phi )d\phi d\lambda$, with $d\phi$ and $d\lambda$ the subtended latitudinal and longitudinal angles of each grid point, respectively.

The resulting phase curves calculated from each of our GCMs using Equation (10) are shown in Figure 6. Here we have normalized the flux for each phase curve by the average flux radiated toward Earth over an orbit. There are two key parameters we can pull out of these simulated phase curves: the phase curve amplitude and the phase offset of the peak emitted flux from secondary eclipse. In general, as expected from our theory, the phase curve amplitude increases with increasing ${T}_{\mathrm{eq}}$ and decreasing ${\tau }_{\mathrm{drag}}$. The phase offset shows the opposite trend, generally decreasing with increasing ${T}_{\mathrm{eq}}$ and decreasing ${\tau }_{\mathrm{drag}}$. Next, we will directly compare our simulated phase curve amplitudes to observations in Section 3.4, and we will explore trends in phase offsets further in Section 3.6.

Zoom In Zoom Out Reset image size

Picking off the maximum and minimum of the mean flux radiated toward Earth and calculating their effective temperatures ${T}_{\mathrm{eff}}={({S}_{\mathrm{av}}/\sigma )}^{1/4}$, we can calculate the phase curve amplitude A from our numerical simulations. Note that for our double-gray simulations, which only have one infrared band, the effective temperature is equivalent to the brightness temperature. These numerical predictions of the phase curve amplitude for varying equilibrium temperature and drag strengths are shown (along with the observations from Figure 1) in Figure 7.

Zoom In Zoom Out Reset image size

As expected from the theoretical predictions in Section 2 and as seen in observations, the numerical phase curve amplitudes increase with increasing incident stellar flux. Additionally, the presence of strong drag can greatly increase the phase curve amplitude. As expected from the theoretical predictions, the day–night temperature difference is significantly larger for the models with ${\tau }_{\mathrm{drag}}\lt {10}^{5}\ {\rm{s}}$ than for those with ${\tau }_{\mathrm{drag}}\geqslant {10}^{5}\ {\rm{s}}$. However, the exact dependence of the phase curve amplitude on ${\tau }_{\mathrm{drag}}$ for the cases with ${\tau }_{\mathrm{drag}}\geqslant {10}^{5}\ {\rm{s}}$ is subtle, as the ${\tau }_{\mathrm{drag}}={10}^{5}\ {\rm{s}}$ case actually shows smaller phase curve amplitudes than the ${\tau }_{\mathrm{drag}}=\infty$ model. This is because the ${\tau }_{\mathrm{drag}}={10}^{5}\ {\rm{s}}$ is a transition case between an equatorial jet and dayside–nightside flow and lacks the cold midlatitude vortices seen on the nightside of the cases with ${\tau }_{\mathrm{drag}}\geqslant {10}^{6}\ {\rm{s}}$ (see Figure 5). Lastly, note that the trend of increasing A with increasing ${T}_{\mathrm{eq}}$ is weaker for our numerical simulations than for our theoretical predictions (compare with Figure 2). This is likely due to the lack of an opacity temperature dependence in our numerical simulations, which, if included, would shift the photosphere to higher pressures at higher equilibrium temperatures and increase the simulated phase curve amplitudes.

We ran a secondary suite of GCMs varying ${T}_{\mathrm{eq}}$ from 500 to 3000 K and fixing ${\tau }_{\mathrm{drag}}=\infty$, while keeping the rotation period at the value it would have if it were equal to the orbital period of a planet orbiting a star with the same bolometric luminosity as HD 189733. We do so in order to analyze the effects of rotation rate on our results from Section 3.2, as in our main grid of simulations we kept the rotation rate fixed. Note that this set of models is similar to the double-gray GCMs performed in Perna et al. (2012), but here we span a slightly larger range of incident stellar flux. The relationship between the orbital period ${P}_{\mathrm{orb}}$ and ${T}_{\mathrm{eq}}$ for a blackbody with zero albedo is

$$\begin{eqnarray}&&{P}_{\mathrm{orb}}=\displaystyle \frac{{\pi }^{1/4}}{4}{\left(\displaystyle \frac{{L}_{\star }}{\sigma }\right)}^{3/4}{({{GM}}_{\star })}^{-1/2}{T}_{\mathrm{eq}}^{-3},\end{eqnarray} \tag{ 11 }$$

where L_⋆ and ${M}_{\star }$ are the stellar luminosity and mass, respectively, and G is the gravitational constant. We calculate ${L}_{\star }=4\pi {R}_{\star }^{2}\sigma {T}_{\star }^{4}$, where ${R}_{\star }$ and T_⋆ are the stellar radius and effective temperature, respectively, and we have assumed that the stellar spectrum can be approximated as a blackbody. Here we take values of stellar radius, mass, and effective temperature relevant to HD 189733: ${R}_{\star }=0.805{R}_{\odot }$, ${M}_{\star }\approx 0.8{M}_{\odot }$, and ${T}_{\star }=4875\ {\rm{K}}$. This results in orbital and rotation periods decreasing from 30.67 to 0.14 days as ${T}_{\mathrm{eq}}$ increases from 500 to 3000 K. Note that the short end of these orbital periods corresponds to a distance from the star of ∼2 stellar radii. As a result, our simulations with very large ${T}_{\mathrm{eq}}$ are included to understand the expected behavior of the circulation in the high-Ω limit. However, note that hot Jupiters with such short orbital periods would likely be destroyed as a result of tidal interactions with their host star. Additionally, planets with long orbital periods are not expected to be tidally locked, and hence this suite of simulations was run solely to understand the change in dynamics with varying ${T}_{\mathrm{eq}}$ and Ω, not to compare with observations.

The results of the suite of simulations consistently varying ${T}_{\mathrm{eq}}$ and Ω are shown in Figure 8. First, note that the atmospheric circulation for ${T}_{\mathrm{eq}}\geqslant 1000\ {\rm{K}}$ is in general similar to the models in Section 3.2, with an eastward equatorial jet dominating the flow pattern. For comparison with the GCMs with fixed rotation rate, we show the dayside–nightside temperature contrast from this grid of simulations along with the results from the larger grid in Figure 7. As in the simulations with consistent rotation rate from Perna et al. (2012), the fractional dayside–nightside temperature differences increase with increasing ${T}_{\mathrm{eq}}$. This is the same trend expected from our theory and found from our grid of simulations with a constant rotation rate shown in Figure 5.

Zoom In Zoom Out Reset image size

The fractional dayside–nightside temperature differences at high ${T}_{\mathrm{eq}}\geqslant 2500\ {\rm{K}}$ match well with those from the simulations with fixed rotation rate. However, the fractional dayside–nightside temperature differences at intermediate $1000\ {\rm{K}}\,\leqslant {T}_{\mathrm{eq}}\leqslant 2000\ {\rm{K}}$ are significantly lower than those from the models with a constant rotation rate. This is likely because of the decreased effect of the midlatitude nightside vortices on increasing the day–night temperature contrast in our runs with consistent rotation rate. This is because the vortices occur at higher latitudes in these simulations with a consistently varying rotation rate than in the simulations with a constant rotation rate.

Note that in simulations with high ${T}_{\mathrm{eq}}\geqslant 1500\ {\rm{K}}$, the temperature and wind patterns are very similar. At such high ${T}_{\mathrm{eq}}$, the maximum zonal-mean zonal wind speed begins to saturate at $\sim 5\ \mathrm{km}\ {{\rm{s}}}^{-1}$. Additionally, the simulations with short rotation periods have a larger equator-to-pole temperature contrast than those in Figure 5. This is because, as discussed by Showman et al. (2015), more rapidly rotating flows can support larger horizontal temperature differences. This is because the increasing Coriolis force with increasing rotation rate can support larger pressure gradients, which are related to the latitudinal temperature difference (see Equation (8) of Showman et al. 2015).

3.6.1. Numerical Results

Post-processing our numerical results from Section 3.2 using the method in Section 3.3, we obtain the offset of the peak of the infrared phase curve along with the amplitude discussed previously. Here, we calculate the offset from secondary eclipse for each broadband infrared phase curve from Figure 6 in order to compare with observational data. Figure 9 shows a comparison between our numerically calculated infrared phase curve offsets as a function of ${T}_{\mathrm{eq}}$ for varying ${\tau }_{\mathrm{drag}}$, along with the observed offsets of the peak of the infrared phase curve for the same low-eccentricity transiting hot Jupiters as in Figure 1. The general trend resulting from our numerical simulations is decreasing eastward phase offsets with increasing ${T}_{\mathrm{eq}}$, as seen from the observations and expected from the theory of Cowan & Agol (2011a) and also from Zhang & Showman (2016). We will compare directly our numerical results with the analytic theory of Zhang & Showman (2016) in Section 3.6.2.

Zoom In Zoom Out Reset image size

Our numerical simulations capture the general trend of observed decreasing infrared phase offsets with increasing incident stellar flux and can explain the entire set of observed infrared phase offsets with varying values of ${\tau }_{\mathrm{drag}}$. However, there is considerable scatter in the observed infrared phase offsets between different wavelength bands, so in the context of our numerical simulations a wide range (${10}^{4}\mbox{--}{10}^{7}\ {\rm{s}}$) of drag timescales must be invoked to explain the observed phase shifts. Considering Lorentz forces as the dominant drag mechanism, this would only be physical if hot Jupiters had a wide range of magnetic field strengths or atmospheric compositions and hence ionization fractions. The former is unlikely if the magnetic field is driven by an internal dynamo (Christensen et al. 2009; Christensen 2010). The latter possibility may have an effect, but note that sodium (the element with the lowest ionization potential) is found in a variety of hot Jupiters regardless of equilibrium temperature (Sing et al. 2016). Additionally, note that the Rayleigh drag formulation that we use is at best a very rough approximation to the true effects of Lorentz forces (Perna et al. 2010; Rauscher & Menou 2013; Heng & Workman 2014; Rogers & Komacek 2014). We will discuss this further in Section 4.

More generally, it is theoretically expected that a competition between the speed of the equatorial jet and the radiative cooling of the circulation determines the phase offset (Cowan & Agol 2011a; Zhang & Showman 2016). Given that drag on the winds acts to slow the equatorial jet (both directly through decelerating the winds and indirectly by affecting the wave action that drives the jets), invoking a variety of ${\tau }_{\mathrm{drag}}$ is sufficient to explain observations in the context of our numerical simulations. However, as discussed above, this is likely unphysical for the case of magnetic drag. In general, it is unclear exactly how the effective drag strength should vary with planetary parameters, largely because the most important drag mechanism in hot Jupiter atmospheres has yet to be identified. Future observations could help to determine whether the observed scatter in infrared phase offsets is real, or whether there is instead a clear trend in phase offset with incident stellar flux.

3.6.2. Comparison with Analytic Theory

Zhang & Showman (2016) developed a simple theoretical model, similar to that of Cowan & Agol (2011a), to determine the phase offset of the planetary hot spot for tidally locked gas giants. To do so, they used a kinematic model in the Newtonian cooling approximation, assuming a horizontally constant radiative timescale and a zonally symmetric jet with a prescribed zonal wind speed. Doing so, they derived the following transcendental relation that determines the hot spot offset for a given equatorial jet speed and radiative timescale (see their Appendix B):

$$\begin{eqnarray}&&\sin ({\lambda }_{\max }-{\lambda }_{{\rm{s}}}){e}^{{\lambda }_{\max }/\xi }=-\displaystyle \frac{\eta }{\xi \cos ({\lambda }_{{\rm{s}}})},\end{eqnarray} \tag{ 12 }$$

where ${\lambda }_{\max }$ is the hot spot offset, $\xi ={\tau }_{\mathrm{rad}}/{\tau }_{\mathrm{adv}}$ is the ratio of radiative to advective timescales (where ${\tau }_{\mathrm{adv}}=a/{U}_{\mathrm{jet}}$ is the timescale for the equatorial jet [with speed ${U}_{\mathrm{jet}}$] to advect across a planetary hemisphere), ${\lambda }_{{\rm{s}}}={\tan }^{-1}(\xi )$, and

$$\begin{eqnarray}&&\eta =\displaystyle \frac{\xi }{1+{\xi }^{2}}\displaystyle \frac{{e}^{\pi /2\xi }+{e}^{3\pi /2\xi }}{{e}^{2\pi /\xi }-1}.\end{eqnarray} \tag{ 13 }$$

Zhang & Showman (2016) compared the prediction from Equation (12) to their numerical simulations using Newtonian cooling, finding good agreement throughout parameter space (see their Figure 9).

Here we examine the predictions from Equation (12) for how phase offset varies with equilibrium temperature in the context of our numerical simulations examined in Section 3.6.1. Figure 10 shows a comparison between the solutions to Equation (12) for varying speeds of the equatorial jet and our numerical calculations of infrared phase offsets. For this comparison, we use the relationship between ${\tau }_{\mathrm{rad}}$ and ${T}_{\mathrm{eq}}$ from Equation (4) to calculate analytic phase offsets as a function of ${T}_{\mathrm{eq}}$.

Zoom In Zoom Out Reset image size

The analytic theory of Zhang & Showman (2016) reproduces the general trend of decreasing phase offset with increasing equilibrium temperature from our GCMs with double-gray radiative transfer. However, note that given that there is no theory for the equatorial jet speeds in hot Jupiter atmospheres (which will be discussed further in Section 4), the analytic jet speeds do not vary self-consistently with ${T}_{\mathrm{eq}}$. The phase offsets from our GCMs do not have a monotonic dependence with temperature, due to nonlinear interactions between the equatorial waves driving the superrotating jet. ³ Additionally, the analytic theory shows a slightly steeper relationship between ${\lambda }_{\max }$ and ${T}_{\mathrm{eq}}$ than the GCM calculations with varying ${\tau }_{\mathrm{drag}}$, and a shallower dependence than our GCM calculations with the rotation rate varying consistently with ${T}_{\mathrm{eq}}$. This is likely due to the lack of dependence of ${U}_{\mathrm{jet}}$ on ${T}_{\mathrm{eq}}$ in the analytic theory. Including this would increase the analytic infrared phase offsets at high ${T}_{\mathrm{eq}}$, as the jet speed should increase with incident stellar flux, decreasing ${\tau }_{\mathrm{adv}}$. Similarly, the analytic theory appears to predict a much larger phase offset than found in our main grid of GCM calculations (without a consistently varying rotation rate) at low ${T}_{\mathrm{eq}}=500\ {\rm{K}}$. This is mainly because the jet speeds at such a low value of incident stellar flux are very small when drag is applied, with a value of ${U}_{\mathrm{jet}}\sim 10\,{\rm{m}}\,{{\rm{s}}}^{-1}$ for ${\tau }_{\mathrm{drag}}={10}^{6}\,{\rm{s}}$ and ${U}_{\mathrm{jet}}\lesssim 1\,{\rm{m}}\,{{\rm{s}}}^{-1}$ for ${\tau }_{\mathrm{drag}}\leqslant {10}^{5}\ {\rm{s}}$.