A GENERAL CIRCULATION MODEL FOR GASEOUS EXOPLANETS WITH DOUBLE-GRAY RADIATIVE TRANSFER

It has been almost one decade since the first atmospheric measurement of a hot Jupiter (Charbonneau et al. 2002), and yet this class of exotic exoplanet still provides us with many mysteries waiting to be solved. These include the culprit responsible for the stratospheric temperature inversions inferred for many hot Jupiters (Hubeny et al. 2003; Fortney et al. 2008; Spiegel et al. 2009; Zahnle et al. 2009; Knutson et al. 2010), the implication of a super-solar carbon-to-oxygen ratio in at least one hot Jupiter (Madhusudhan et al. 2011), the possible influence of magnetic drag on the atmospheric circulation (Perna et al. 2010a), and a peculiar planet where dawn seems to be much hotter than noon (Crossfield et al. 2010). In the midst of a growing collection of observations, some with unexpected or surprising results, is the continuing development of atmospheric models trying to interpret and understand the measurements.

One set of atmospheric models are numerical ones that simulate the global circulation patterns on close-in gas giants. With masses and radii comparable to Jupiter, but subject to incident stellar fluxes 10,000 times stronger than what Jupiter receives from the Sun, and expected to be tidally locked into synchronous orbits, hot Jupiters exist in an atmospheric regime unlike anything in our solar system, and models of their atmospheric circulation are expanding into uncharted territory. In order to try to understand how atmospheres work in this new regime, the models that have been developed so far represent a range of complexities, various approaches, and the use of different sets of assumptions. Some of the simplest models solve the shallow water or equivalent barotropic equations and use various schemes to include the effect of radiative heating (Cho et al. 2003, 2008; Langton & Laughlin 2007). Others solve the primitive equations of meteorology, using either a Newtonian relaxation scheme for the radiative forcing (Showman & Guillot 2002; Cooper & Showman 2005; Showman et al. 2008; Menou & Rauscher 2009; Rauscher & Menou 2010), a dual-band radiative transfer scheme (Heng et al. 2011, similar to the one we present here), or more complex non-gray radiative transfer (Showman et al. 2009). Finally, there are also models that solve the full set of fluid equations, using dual-band flux-limited diffusion for the radiative transfer (Dobbs-Dixon & Lin 2008; Dobbs-Dixon et al. 2010).

Here we present an updated version of our previous general circulation model that now includes a "double-gray" radiative transfer scheme. Fluxes throughout the atmosphere are separated into optical and infrared components, each with its own absorption coefficient, which is constant for the optical band and can scale as a power law with pressure for the infrared band. We use standard two-stream radiative transfer equations to solve for the vertical fluxes and calculate heating rates from those. In addition to producing self-consistent radiative fluxes and heating rates, including the observable infrared flux emerging from the top boundary, this new code has the advantage of maintaining only a moderate level of complexity. This facilitates comparison between simulated results and analytic profiles (Guillot 2010, see also Hansen 2008) and makes it easier to clearly identify the effects due to changes in opacity.

Our new code is described in detail in Section 2. We demonstrate the functionality of our code by presenting a model of a generic hot Jupiter (Section 3) and investigate issues related to the atmosphere's behavior at deep pressures (Section 3.1). We then study the effect of magnetic drag on atmospheric circulation by applying a simplified drag scheme to our model, with drag strengths ranging from weak to strong (Section 4). In Section 5 we summarize our results.

The new version of our code presented here solves the primitive equations of meteorology ⁴ using the dynamical core originally developed by Hoskins & Simmons (1975), converted for hot Jupiter studies in Menou & Rauscher (2009) and used in subsequent papers (Rauscher & Menou 2010; Perna et al. 2010a, 2010b; Miller-Ricci Kempton & Rauscher 2011). The previous version of this code (IGCM1) heated the atmosphere through a Newtonian relaxation scheme. Here we use an updated version of the code (IGCM3, de F. Forster et al. 2000) with several improvements, including a simple radiative transfer scheme, which we have adapted to model gaseous exoplanets. This code also contains a standard dry convection scheme. The static stability throughout each vertical column is calculated, and any layers found to be convectively unstable are adjusted to constant potential temperature (entropy), with the vertical integral of enthalpy conserved.

2.1. Radiative Transfer

The primary driver of atmospheric circulation on close-in extrasolar planets is heating from the stellar irradiation, which is absorbed and re-emitted throughout the atmosphere. We employ the common assumption that the incident radiation is predominantly at optical wavelengths, while the re-emitted radiation is in the infrared. We can then separate our treatment of radiative transfer into two wavelength ranges, with downward optical flux absorbed as it penetrates deeper into the atmosphere and infrared flux absorbed and emitted isotropically throughout the atmosphere. We calculate the upward and downward net fluxes in each vertical column of the atmosphere and then solve for the localized specific heating as dT/dt = (g/c_p )dF/dP, where dT/dt is the temperature tendency, g is the gravitational acceleration (assumed to be constant through the model), c_p is the specific heat, and dF/dP is the vertical derivative of the net vertical flux (pressure, P, is our vertical coordinate). Upward flux, in the direction of decreasing pressure, is defined to be positive. For an atmosphere in local radiative equilibrium (i.e., without dynamics) dF/dP = 0, either because there is a constant flux through the atmosphere (e.g., an outward heat flux from the planet interior) or because the optical heating is exactly balanced by infrared cooling. The heating rate feeds into the energy equation of the dynamics solver (Q in Equation (1d) of Menou & Rauscher 2009).

The incident optical flux at the top of each column of the atmosphere, with latitude ϕ and longitude θ, is F(ϕ, θ) = (1 − A)F_inccos (ϕ)cos (θ) for the day side ⁵ and F = 0 on the night side, where F_inc is the stellar flux incident on the top of the atmosphere at the substellar point and A is the Bond albedo, which here we assume to be zero. We assume a well-mixed optical absorber throughout the atmosphere and negligible scattering, so that on the day side the downward flux in a column cos ⁻¹μ_⋆ degrees away from the substellar point (μ_⋆ = cos (ϕ)cos (θ)) is given by the well-known equation (e.g., Stephens 1984)

$$\begin{eqnarray} F_{\downarrow \mathrm{vis}} (P) &=& (1-A) \mu _\star F_{\mathrm{inc}} \exp \left(-\frac{1}{\mu _\star } \int _z^\infty \kappa _{\mathrm{vis}}du \right) \nonumber \\ &=& (1-A) \mu _\star F_{\mathrm{inc}} \exp \left(-\frac{1}{\mu _\star } \frac{\kappa _{\mathrm{vis}}}{g} P\right), \end{eqnarray} \tag{ 1 }$$

where we have converted from path length (du = ρdz) to pressure using the equation for vertical hydrostatic equilibrium (assumed in the primitive equations) and by integrating from the top of the atmosphere (z = ∞, P = 0) down to the current pressure. Note that the absorption of the incident flux, as well as the resulting heating of the atmosphere, is entirely determined by the wavelength-integrated absorption coefficient κ_vis. Lacking any reflection from a solid surface or clouds, the upward optical flux is always zero (F_↑vis = 0).

Our assumption that all of the incoming stellar irradiation is at visible wavelengths leads to the upper boundary condition for our infrared flux: F_↓IR(P = 0) = 0. The bottom boundary condition is an upward heat flux from the interior: F_↑IR(P = P₀) = σT⁴ _int. Throughout the atmosphere upward and downward infrared fluxes are absorbed and emitted as for a gray atmosphere, with a wavelength-integrated absorption coefficient that scales with pressure as a power law: ⁶

$$\begin{equation} \kappa _{\mathrm{IR}}= \kappa _{\mathrm{IR},0} (P/P_{\mathrm{ref}})^{\alpha }. \end{equation} \tag{ 2 }$$

The infrared fluxes are

$$\begin{equation} F_{\uparrow,\downarrow \mathrm{IR}}(P) = \int \left(1-\exp \left[-\frac{1.66}{g} \int \kappa _{\mathrm{IR}}{\it dP}\right]\right) \frac{d \sigma T^4}{{\it dP}} {\it dP}, \end{equation} \tag{ 3 }$$

where the integrations are from the bottom boundary (P = P₀) to the given pressure for the upward flux and from the top boundary (P = 0) for the downward flux, and we have used the standard diffusivity factor 1.66 to account for integration of the isotropic radiation over all zenith angles (see, e.g., Stephens 1984).

If both the infrared and optical absorption coefficients are constant with pressure (α = 0), our radiative scheme (with the dynamics turned off) should reproduce the analytic profiles from Guillot (2010, see also Hansen 2008). Note that for this type of "double-gray" atmosphere with constant absorption coefficients, the temperature–pressure profiles can never be convective; as shown in Appendix B, for a convective region to exist the power-law index, α, must be greater than a threshold value based on the ratio of the specific gas constant to the specific heat (R/c_p ).

We initialize our model runs with no winds and with temperatures everywhere set to the analytic nightside profile (Equation (27) of Guillot 2010). For the first (planet) day of the run, the fluxes from the radiative transfer scheme are not used and instead a Newtonian relaxation scheme heats each point in the atmosphere toward the value it would have in local radiative equilibrium, with a relaxation timescale set to 0.1 planet days. During this time, the atmosphere is also allowed to respond by developing winds. After the first day, the Newtonian scheme is no longer used and the heating is entirely determined by fluxes from the radiative transfer scheme.

To test the implementation of our radiative transfer in the code, we set the infrared opacity to be constant with pressure (α = 0 in Equation (2)) and compared our temperature–pressure profiles against the analytic solutions from Guillot (2010). ⁷ We ran several tests in which the downward optical flux at the top boundary was identical at all points around the globe, so that there are no horizontal gradients and no atmospheric motion. This effectively turns the code into a one-dimensional radiative transfer solver. For a given set of absorption coefficients and boundary fluxes, we began with initial profiles set to the analytic solutions and then let the code run until each column was in local radiative equilibrium: dF/dP = 0, with heating equal to cooling, or with a constant flux set by the bottom boundary condition: F_{↑, IR} = σT⁴ _int. We considered a profile to be equilibrated once it was within less than 1% of dF/dP = 0.

This testing showed that, although our profiles matched the analytic solution at low pressures, they tended toward isothermal in the deepest levels, rather than approaching an adiabat (as shown in Figure 1). The reason for this is that our numerical domain spans many orders of magnitude in pressure and vertical levels are evenly distributed in log P. Deep in the atmosphere the pressure difference between adjacent levels becomes very large (dP ∼ P), the optical depth greatly exceeds unity (τ = κ_IR P/g), and the exponential term in Equation (3) goes to zero. When the deepest levels are in radiative equilibrium, the temperature profile is entirely determined by the vertical resolution, with the temperature differences between adjacent pressure levels set by Δ(σT⁴) = σT⁴ _int (a profile that exactly matches those shown in Figure 1).

Zoom In Zoom Out Reset image size

2.2. Transition to Flux-limited Diffusion at High Optical Depth

In order to achieve the correct temperature–pressure profiles at depth, we have implemented a scheme in which we transition to calculating the radiative flux using the diffusion approximation at optical depths greater than one. The flux through the optically thick, deep atmosphere is calculated as

$$\begin{equation} F_{\mathrm{IR,diff}} = \frac{16 g \sigma T^3}{3 \kappa _{\mathrm{IR}}} \frac{dT}{{\it dP}}. \end{equation} \tag{ 4 }$$

A temperature profile that increases with pressure results in an upward flux (F_IR > 0).

The transition from fluxes calculated using the main radiative transfer equations (Equations (1) and (3)) and the diffusion approximation (4) occurs based on the local optical depth (τ = κ_IR P/g, taking into account any dependence of κ_IR with pressure). The original radiative scheme begins to fail when the term E = exp (− [1.66/g]∫κ_IR dP) becomes too small, and so we use this as a criterion for when to begin the implementation of the flux diffusion scheme. Through trial and comparison to the analytic solutions, we determined that the fix should begin at the point where E ⩽ 0.01, which in our code corresponds to an optical depth τ > τ_limit = 5.55/(10^A − 10^−A ), where A is the vertical resolution of the model, A = log₁₀(P_max/P_min)/NL. For a model whose bottom boundary is at P_max = 100 bar, whose upper boundary is P_min = 1 mbar, and that has NL = 30 levels, τ_limit = 7 (=0.56 bar in Figure 1). A model using only 10 vertical levels would have τ_limit = 2 (=0.16 bar in Figure 1). This resolution-dependent τ_limit appropriately fixes the resolution dependence of the problem with the original radiative transfer scheme.

At optical depths greater than τ_limit, the infrared flux is calculated as F_IR(P) = (E^0.023)F_{IR, rad} + (1 − E^0.023)F_{IR, diff}, where F_{IR, rad} is the infrared flux calculated from Equation (3). The transition to the diffusion fluxes has the same exponential dependence as the problem with the original radiative scheme, and the factor of 0.023 was chosen so that at the pressure where E = 0.01 (τ = τ_limit), F_IR = 0.9 × F_{IR, rad} + 0.1 × F_{IR, diff}. We found that this correction produced a smooth transition (in both fluxes and the temperature–pressure profiles) from the optically thin to optically thick regimes.

2.3. Advantages of This Scheme

We present a fiducial model for a hot Jupiter with global parameters similar to HD 209458b, as listed in Table 1. We chose HD 209458b parameters so that we could use absorption coefficients and internal and external heat fluxes based on values from Guillot (2010) and, in Section 4, magnetic drag timescales from Perna et al. (2010a). Although there is evidence for a temperature inversion in the atmosphere of HD 209458b (Knutson et al. 2008), we do not include stratospheric absorbers in our model, in the interest of generality.

Table 1. Parameters for our Fiducial Model

Parameter	Value	Units
Radius of the planet, Rp	1 × 108	m
Rotation rate, Ω	2.1 × 10−5	s−1
Gravitational acceleration, g	8	m s−2
Specific gas constant, $\mathcal {R}$	3523	J kg−1 K−1
Ratio of gas constant to heat capacity, $\mathcal {R}/c_P$	0.286
Optical absorption coefficient, κvis	4 × 10−3	cm2 g−1
Infrared absorption coefficient, κIR, 0	1 × 10−2	cm2 g−1
Infrared absorption power-law index, α	0
Infrared absorption reference pressure, Pref	1	bar
Internal heat flux, F↑IR, int	3500	W m−2
Corresponding temperature, Tint	500	K
Incident flux at substellar point, F↓vis, irr	1.06 × 106	W m−2
Corresponding temperature, Tirr	2078	K

Download table as:  ASCIITypeset image

We used a horizontal spectral resolution of T31, which corresponds to ∼4° in latitude and longitude, and used 30 vertical levels, evenly spaced in pressure from 1 mbar to 100 bar. We ran the model for 2000 planet days (=2000 orbital periods), by which point the kinetic energy and average temperature of the atmosphere have stabilized, after the initial spin-up period. All of the results shown in this paper, both for this fiducial model and those in later sections, are for a snapshot of the atmosphere at day 2000. In terms of variability, there is a very minor level of quantitative change from day to day, and on large scales the atmospheric structure is very steady.

We performed partial runs at T21 and T42 and found consistency with the T31 results shown here, but the T42 resolution was too computationally expensive for this work. We also tested the strength and order of the hyperdissipation we used to prevent the buildup of noise on small scales, which takes the form dX/dt = ν_diss(− 1)^{p + 1}∇^2p X, where X represents the vorticity, divergence, and temperature fields. For horizontal resolutions of T21, T31, and T42, we tested orders of p = 2, 3, and 4 and values of ν_diss corresponding to e-folding times (for damping the smallest scales) at fractions of a planet day: 0.5, 0.05, 0.005, and 0.0005. The flow quickly became numerically unstable for a damping time of 0.5 and exhibited an obvious buildup of kinetic energy at the smallest scales for 0.05. A comparison of the kinetic energy spectra revealed that the p = 2 order was overdamping kinetic energy on all scales. We chose to use an order of p = 4 and a damping timescale of 0.005 planet days for our T31 fiducial model. Since the radiative timescales throughout our model vary by orders of magnitude, the upper levels may be underdamped and the deepest levels overdamped (Thrastarson & Cho 2011), but we use this form for hyperdissipation until a more robust and physically motivated scheme is developed.

Figure 2 shows temperature–pressure profiles throughout the atmosphere for our fiducial model. Also plotted are the analytic, local radiative equilibrium profiles for the substellar point and the night side (previously seen in Figure 1). As is typical of previous hot Jupiter models, there is a strong day–night temperature difference at high altitude, which decreases with increasing pressure. At moderate pressures (∼1 bar) the primary temperature difference is between the equator and poles. There is also a dynamically created temperature inversion at these pressures. At low pressures the day side remains at temperatures close to radiative equilibrium, but at deeper pressures the winds are better able to alter the temperature structure and the day side is cooled from more efficient homogenization with the night side. Figure 2 also shows that our use of the flux-limited diffusion scheme produces the correct behavior in the pressure–temperature profiles at depth.

Zoom In Zoom Out Reset image size

Figure 3 shows the zonal average of the zonal wind for our fiducial model, with the equatorial super-rotating jet that is standard in hot Jupiter models (Showman & Polvani 2011). This jet extends from the uppermost levels to almost the bottom of the atmosphere, with only the two deepest levels showing westward equatorial flow. Whereas at high altitudes the equatorial jet is accompanied by strong flow (zonal and meridional) from day to night across most regions of the terminator, the equatorial jet dominates the flow for levels below the optical photosphere. Also in agreement with most hot Jupiter models, we find strongly supersonic winds high in the atmosphere (at several times the sound speed), a transonic flow at moderate pressures, and subsonic (only) flow at pressures greater than ∼10 bar.

Zoom In Zoom Out Reset image size

In Figure 4 we plot vertical velocities throughout the model. Although the primitive equations used in this code replace the vertical momentum equation with hydrostatic balance, vertical motions are still present through the continuity equation: $d\omega /{\it dP} = -\nabla \cdot \vec{v}$, where ω = dP/dt and $\vec{v}$ is the horizontal wind. We integrate the continuity equation down from the top of the model, where the boundary condition imposes ω = 0 at P = 0. We then convert to a vertical velocity (dz/dt) through use of the hydrostatic equation: dP/dz = −ρg. We note that snapshots of the vertical wind profiles taken at other points in the run match the main features shown here.

Zoom In Zoom Out Reset image size

The strength of vertical motion in the atmosphere has important implications for the mixing of chemical species throughout the atmosphere, which can introduce chemical disequilibrium and influence the observable properties of the planet (Moses et al. 2011). For example, titanium oxide (TiO) and vanadium oxide (VO) have been identified as possible candidates for the absorbing species responsible for stratospheric temperature inversions, but Spiegel et al. (2009) demonstrated that very strong vertical mixing would be required to keep these absorbers aloft in the atmosphere. From Figure 4 we can see that vertical speeds in our model are generally very slow, at tens of meters per second (compared to horizontal wind speeds of ∼1 km s⁻¹), although in localized regions they can exceed 100 m s⁻¹. Throughout most of the atmosphere the vertical winds are systematically stronger on the night side than the day side, although at pressures less than ∼10 mbar the dayside winds are slightly stronger. There is a trend toward higher vertical speeds at lower pressures, and the polar regions have the strongest vertical motion, upward on the day side and downward on the night side. The direction of vertical motion is consistent with the substellar-to-antistellar flow at low pressures; the divergent flow on the day side causes upward vertical motion, while convergent flow on the night side results in downward motion. The profiles in Figure 4 demonstrate that the strength of vertical mixing depends on the local atmospheric flow structure. However, the horizontal wind speeds are much stronger than vertical wind speeds, so that horizontal mixing can also bring species out of equilibrium (Cooper & Showman 2006) and can add to the difficulty of keeping TiO/VO aloft (Showman et al. 2009). The amount of chemical mixing throughout the atmosphere is inherently a three-dimensional problem.

Figure 5 shows the temperature and flow pattern at the infrared photosphere (τ = 2/3, P = 53 mbar), as well as a map of the infrared flux emitted from the top boundary of the model. (Note that under our model assumptions, including that of zero albedo, the infrared flux encompasses all of the light that would be observed from the planet and so can also be thought of as the bolometric flux.) We again see typical hot Jupiter atmospheric structure, where there is a significant day–night temperature difference, but the equatorial jet has advected the hottest region of the atmosphere eastward of the substellar point. Our radiative scheme self-consistently produces the infrared map that would be observed for this planet, from the upward infrared flux at the top boundary of the model.

Zoom In Zoom Out Reset image size

One of the advantages of this new code is that the radiative scheme self-consistently solves for the infrared and optical fluxes throughout the atmosphere and calculates the resulting local diabatic heating throughout the three-dimensional atmosphere. In Figure 6 we show the fluxes and heating for various regions of the atmosphere. The optical flux is only incident on the day side, is always downward, and heats the atmosphere, while the infrared flux is in both directions throughout the planet but is primarily upward and cools the atmosphere.

Zoom In Zoom Out Reset image size

For an atmosphere in local radiative equilibrium, at high pressures (below the optical photosphere) there would be a non-zero cooling flux, set by the boundary condition of upward flux from the planet interior (see Figure 1). The incident stellar flux sets the top boundary condition, and high in the atmosphere the optical heating flux is mostly balanced by the infrared cooling flux, with a slight offset due to the cooling flux coming from the interior. The heating at each level in the atmosphere is the difference between the net fluxes into and out of the layer and, by definition, would be near zero everywhere for local radiative equilibrium (modulo the small interior heat flux leaking out through the atmosphere).

The left panel of Figure 6 shows that the interplay between radiation and atmospheric dynamics results in flux profiles altered from what we would expect for local radiative equilibrium. At the bottom of the atmosphere, the cooling flux is close to the upward flux of 3500 W m⁻² from the planet interior (set as a boundary condition). In the upper atmosphere the infrared and optical fluxes are close to being globally balanced. The top two plots show that the optical flux only passes through the day side of the planet (μ_⋆ = cos (ϕ)cos (θ) > 0) and that the regions closest to the substellar point (purple line) experience a stronger optical flux that penetrates deeper than areas close to the terminator (green line). The flux at the top of the model is set by the incident stellar flux boundary condition, here chosen to be 1.06 × 10⁶ W m⁻² at the substellar point (and decreasing as μ_⋆, according to Equation (1)).

There is a net heating flux throughout most of the day side (aside from the region near the terminator), with the infrared cooling fluxes unable to completely balance the optical heating. The total dayside heating flux extends deeper than the optical photosphere because of the contribution from infrared flux for a few layers around 1 bar. This is due to the dynamically induced thermal inversion at these pressures (see Figure 2), which leads to the net downward transport of infrared radiation.

The right panel of Figure 6 shows horizontally integrated diabatic heating rates throughout the atmosphere, which are calculated from the fluxes shown in the left panel; locally dT/dt = (g/c_p )dF/dP, where upward flux (in the direction of lower pressure) is defined to be positive (see Section 2.1). When flux increases with increasing pressure, there is heating; this is the case for the downward (negative) optical flux at the substellar point, which goes to zero as pressure increases. Likewise, in the upper atmosphere the upward (positive) infrared flux near the antistellar point decreases with decreasing pressure, which heats this region of the night side. For the globally averaged rates, the optical flux always leads to heating while the infrared generally leads to cooling. The net globally averaged profile has cooling throughout most of the atmosphere, aside from a region near the optical photosphere where there is net heating.

The global net heating rate is 3.9 × 10²¹ W. A planet in global radiative equilibrium should have a net radiative heating rate of zero. However, for a steady-state atmosphere, the dissipation of wind kinetic energy (numerical or physical) must be balanced by a non-zero conversion from thermodynamic potential energy, which in turn must be generated by a non-zero heating rate (e.g., Pearce 1978). Physically, the dissipation of kinetic energy leads to frictional heating and the energetics are balanced. However, in our fiducial model the numerical dissipation of kinetic energy (through hyperdissipation or other numerical loss) is not returned to the atmosphere as heating and so the net radiative heating must be non-zero to compensate. We can take the ratio of the net heating rate and the rate of energy input into the model (the total of the incident stellar flux and the cooling flux from the interior) to determine that in this model the rate of numerical loss of kinetic energy is 12%, a significant value and similar to the numerical loss found in our previous version of this code (Rauscher & Menou 2012).

3.1. The Effects of Internal Heating and Flux Diffusion on the Deep Atmosphere

The temperature structure of the upper atmosphere, above the optical and infrared photospheres, is dominated by the external stellar heating. At deeper pressures, however, the atmospheric structure is set by the strength of the heat flux from the interior or, equivalently, the entropy of the interior adiabat. Here we present two additional models to examine how the behavior of the deepest levels of our model depends on our choice of T_int and our use of the flux diffusion scheme (Section 2.2).

The precise value for a planet's interior heat flux is difficult to disentangle from other uncertainties in system and observed parameters. We chose T_int = 500 K for our fiducial model in order to match the analytic setup of Figure 3 from Guillot (2010), but values of T_int down to ∼100 K may be reasonable. As a test of the influence of the choice of T_int on the atmospheric circulation, we ran a model with T_int = 0 K (and all other parameters identical to our fiducial run), with the resulting temperature structure shown in the left panel of Figure 7.

Zoom In Zoom Out Reset image size

With no heating from the interior, a night side in radiative equilibrium would have zero temperature. In our T_int = 0 K model, all of the nightside heating is the result of advection from the day side. The temperatures in the upper atmosphere are still very similar to our fiducial model, a result that should not be surprising, given that the stellar heating dominates at low pressures. The radiative fluxes for this model match those from our fiducial model (see Figure 6) for pressures less than ∼1 bar (both globally and for different regions of the atmosphere). Below the optical photosphere our fiducial model has upward fluxes of ∼3500 W m⁻² (the flux corresponding to our choice of T_int = 500 K), while in this T_int = 0 K model the fluxes quickly decrease toward zero. The heating rates for this model match our fiducial model almost identically. The temperature structure of the deepest levels of this model remains fairly isothermal (at least along the equator).

Similarly, the model run without the flux diffusion scheme (shown in the right panel of Figure 7) also has an upper atmosphere whose temperature structure closely matches our fiducial model, but with an isothermal deep atmospheric structure. Both of these tests show that the upper atmosphere (at the pressures directly probed by observations) is relatively insensitive to the conditions in the deeper atmosphere.

However, any coupling between the atmosphere and planet interior will be sensitive to the temperature structure at deep pressures. The isothermal behavior at deep pressures in these models means that the atmosphere has a high static stability and there is less efficient coupling between layers than in our fiducial model, where the atmosphere trends toward an adiabat. ⁸ While it is not yet clear what is the best way to set the bottom boundary condition for hot Jupiter models, the temperature profiles at deep pressure will control how efficiently the atmosphere can transport energy and momentum vertically.

We note that recent work by Heng et al. (2011) uses a code very similar to the one we have developed to study a range of classes of atmospheres, including hot Jupiters. They set their internal heat flux to zero but choose as a bottom boundary condition to have their "surface" layer be an idealized slab with uniform temperature and finite heat capacity, in order to simulate interaction with the planet interior. The fluxes into and out of the bottom of their model interact with this slab by raising and lowering its temperature. Thus, while similar, their model setup is not identical to the one we show here.

Recent work suggests that hot Jupiter atmospheres should be weakly thermally ionized and—depending on the strength of the planetary magnetic field—the influence of magnetic drag on the atmospheric winds could be strong enough to significantly alter the atmospheric circulation (Perna et al. 2010a). In the near future, observations may be able to constrain the speed of atmospheric winds by precisely measuring the associated Doppler shift present in the transmission spectra of hot Jupiters, giving us an idea of the strength of drag (magnetic or otherwise) in these atmospheres (Miller-Ricci Kempton & Rauscher 2011). This effect is also directly related to the structure of the planet as a whole, as ohmic dissipation of the induced currents could provide an additional heating source and possibly explain the unexpectedly large radii of some hot Jupiters (Batygin & Stevenson 2010; Batygin et al. 2011; Perna et al. 2010b). Simple scaling laws suggest an anti-correlation between two observable properties: the amount by which the hottest region of the atmosphere is offset from the substellar point and the degree of radius inflation (Menou 2012), an issue that can be studied in more detail with numerical models.

The original models that we used to estimate the effect of magnetic drag on hot Jupiter circulation were run with the previous version of our code, in which the radiative heating was determined by a Newtonian relaxation scheme (Perna et al. 2010a). In that work we estimated an average value for the magnetic drag timescale at each pressure level (τ_drag), for planetary magnetic field strengths of 3, 10, and 30 G. In deriving the strength of the magnetic drag, we assumed a dipole planetary magnetic field aligned with the rotation axis, that the thermal ionization was primarily due to potassium, and we compared terms in the induction equation to conclude that the atmosphere was in the purely resistive MHD regime (further details on the assumptions we use can be found in Liu et al. 2008; Perna et al. 2010a; Menou 2012). We then simulated the effect of magnetic drag on the winds by applying a Rayleigh drag to our model, with the form dv/dt = −v/τ_drag. Note that we have applied drag to both the meridional (north–south) and zonal (east–west) winds, although according to our assumptions the drag should only act on the zonal winds. Throughout most of the atmosphere zonal winds dominate the flow, but this treatment of the drag could lead to inaccurate results high in the atmosphere (P ≲ 10 mbar), where the meridional wind speeds become comparable to the zonal wind speeds. The weakest drag (for B = 3 G) had timescales ranging from ∼6 × 10⁶ s at 1 mbar to ∼8 × 10⁸ s at 220 bar, with the medium (B = 10 G) and strong (B = 30 G) drag cases having timescales 10 and 100 times shorter, respectively. Here we employ the same form for the magnetic drag, but now with our upgraded model that includes radiative transfer. This form for the magnetic drag oversimplifies the complex physics but gives us an initial estimate of how it may affect the circulation. Future work will include an improved, more physically accurate scheme for calculating the magnetic drag from local conditions, in particular allowing for horizontally non-uniform drag.

In Figure 8 we show the zonal wind profiles for models with weak, medium, and strong drag. By comparison with our fiducial, drag-free model (Figure 3) we can see that increased drag has three effects on the flow: (1) slower wind speeds, (2) less of a disparity in strength between eastward and westward winds, and (3) the confinement of the super-rotating equatorial jet to higher pressures. This is similar to our results in Figure 3 of Perna et al. (2010a). Note that the model with the strongest amount of drag still has supersonic winds in the upper atmosphere (∼7 km s⁻¹), even though the zonally averaged wind speeds in Figure 8 are subsonic (c_s ≈ 2–3 km s⁻¹); this is because there is strong substellar-to-antistellar flow and the eastward and westward winds cancel each other in the zonal average.

Zoom In Zoom Out Reset image size

The structure of the horizontal flow changes with increasing drag strength, as does the structure of the vertical flow. In Figure 9 we show profiles of the vertical winds throughout the model with strong drag. Comparing with the fiducial model (Figure 4), we see that two of the main trends remain: vertical motion is generally stronger at lower pressure and is strongest at the poles. Some of the details of the flow do change; in particular, at deep pressures we see slower winds and less of a difference between average speeds on the day and night sides. Finally, although the horizontal wind speeds have been reduced by the strong drag, the drag does not act on vertical motion and the vertical winds at low pressure are slightly stronger in this model than in the fiducial one. This could lead to a slight change in the ability of the atmosphere to vertically mix chemical species.

Zoom In Zoom Out Reset image size

In the model with strong drag, the circulation pattern is transitioning away from the behavior seen in the low- and no-drag models. At high altitude the substellar-to-antistellar flow dominates, with the eastward equatorial jet confined to the night side. In fact, even though the zonal average of equatorial winds is eastward for much of the atmosphere, there are no pressure levels at which eastward equatorial flow is able to completely circle the globe. In Figure 10 we plot the temperature and wind structure at the photosphere, for the strong-drag model. The slower wind speeds and the transition to a primarily substellar-to-antistellar flow structure result in a hot day/cold night temperature structure. In the strong-drag model the hottest region of the atmosphere remains at the substellar point, rather than being advected eastward, as is the case for our fiducial model (compare to Figure 5).

Zoom In Zoom Out Reset image size

As a result of the change in the circulation pattern, we find, as expected, a correlation between increased drag strength and a decreased shift of the hotspot away from the substellar point; in Figure 11 we show maps of the infrared flux emitted by the planet for the weak-, medium-, and strong-drag models. By integrating over the visible portion of the planet, for viewing angles along the equator, we calculated the infrared luminosity that would be observed throughout the planet's orbit, for our fiducial and the three drag models (Figure 12). We used a single model snapshot, instead of snapshots of the planet throughout an orbit, but the level of atmospheric variability is low enough in all models that the shape of the light curve should be unaffected.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We find very little difference between the phase curves for the fiducial, weak-, and medium-drag models; however, the peak in flux from the strong-drag model is well aligned with an orbital phase of 0.5 (when the substellar point is directed toward the observer). The peak in luminosity occurs at orbital phases ranging from 0.467 for the drag-free model to 0.494 for the strong-drag model. It may be challenging to observationally distinguish between these models with current instruments. The best hot Jupiter thermal phase curves so far can only constrain the phase of peak flux to ±0.02 (Knutson et al. 2007, 2009).

The phase curve for the strong-drag model is also more sharply peaked over the day side and has a dimmer night side than the other three models, leading to a larger flux contrast between the day and night sides. The fiducial model has a minimum flux that is 17.4% of the maximum, while this ratio is 12.9% for the strong-drag model. The most sensitive phase curve observations can measure this ratio to within ±3%–7% (Knutson et al. 2007, 2009), meaning that this measurement may also be challenging, if one is to discriminate between these models.

In order to compare the observable properties of radius inflation and hotspot offset (Menou 2012), we must estimate the ohmic heating that would be produced in these magnetic drag models. Although we do not explicitly calculate the heating from ohmic dissipation in these models, it should be equal to the amount of kinetic energy dissipated by magnetic drag, which we can calculate from the kinetic energy fields and prescribed drag timescales. We report this value for each of our models in Table 2. As a percentage of the stellar heating input to our models, we find that the ohmic heating in our weak-, medium-, and strong-drag models is at efficiencies of 0.6%, 3%, and 60%. The amount of extra heating necessary to significantly inflate a hot Jupiter's radius depends, in part, on the pressures at which it is deposited (e.g., Guillot & Showman 2002). Ohmic heating is a non-local process, meaning that the dissipation of induced currents does not necessarily happen in the same location as where magnetic drag brakes the winds, and in this paper we do not estimate where the heating may occur; however, our previous work indicated that enough of the heating may take place deep enough in the atmosphere to contribute to radius inflation (Perna et al. 2010b). Although it is not clear if the weak-drag model would produce enough ohmic heating to be significant, the 3% efficiency found in our medium-drag model should be able to inflate the planetary radius, while the 60% efficiency of our strong-drag model is so high that this heating may be able to completely evaporate the planet (according to the models of Batygin et al. 2011).

It is important to note that our strong-drag model is not self-consistent, in that the ohmic heating is not explicitly included. Regardless of whether the heating was deposited locally to where the drag occurs (and thus would drastically change the temperature structure, amount of ionization, and wind profiles high in the atmosphere), or whether the heating was deposited deeper inside the planet (and so would change the overall static stability of the atmosphere, which could in turn affect the flow pattern of the general circulation), the self-consistent inclusion of heating from ohmic dissipation would alter the atmospheric structure, the amount of kinetic energy dissipated by drag, and the final calculation of total ohmic heating in the planet. A better analysis of the correlation between the eastward shift of the hottest region of the atmosphere and the amount of radius inflation must wait for an improved modeling scheme that accounts for magnetic drag and its associated ohmic heating self-consistently.

Finally, we calculate the rate of numerical kinetic energy loss in each of our drag models. For our fiducial model the rate of kinetic energy loss was found by comparing the non-zero global net radiative heating to the rate of energy input to the atmosphere, assuming that the net heating was balancing the net numerical loss. The magnetic drag is an explicit sink of kinetic energy, and we subtract this value from the net radiative heating, assuming that the remainder is balancing numerical dissipation. As we report in Table 2, the amount of numerical loss is similar between our fiducial, weak-, and medium-drag models but is zero (to within 1%) in our strong-drag model. It may be that the similarity between the phase curves of the no-, weak-, and medium-drag models is a result of the strong role of numerical dissipation in these models, overshadowing any differences due to the applied magnetic drag. In addition, given the high rate of numerical kinetic energy loss in these models, the ohmic heating rates that we have estimated may be inaccurate; if the loss did not occur and the wind speeds were higher, we would expect for there to be more heating. These issues will need to be addressed in future work.

An exact identification of the source of numerical loss in our models (e.g., whether it is localized or distributed uniformly) is beyond the scope of this paper; however, the lack of any significant loss in our strong-drag model seems to indicate that the numerical loss is related to the much stronger wind speeds in our other models. This is also corroborated by the results from our initial runs where we tested various settings for the strength and order of the hyperdissipation. In models run at lower resolution (T21) we found a similar rate of loss to the models presented here (∼10%–20%), with the spread due to differences in the hyperdissipation used. We found that the weaker the hyperdissipation, the more kinetic energy remained in the atmosphere, and the higher the rate of numerical kinetic energy loss. It may be that we can substantially reduce numerical loss by correctly including all sources of drag in our models, be it magnetic or related to subgrid processes (see, e.g., Li & Goodman 2010), but for now this remains an issue of concern.

We have presented a new radiative transfer scheme for our atmospheric circulation code. It divides all radiation into optical and infrared wavelengths. The optical flux is incident on the top boundary, and its attenuation is controlled by a constant optical absorption coefficient. The infrared flux is absorbed and emitted at each level, as governed by an infrared absorption coefficient that goes as a power law with pressure (and can be constant, as we have chosen for all of the results presented in this paper). Below the infrared photosphere the infrared fluxes are calculated using the flux-limited diffusion approximation, which we found was necessary in order to reproduce analytic temperature profiles correctly.

Using this new code, we present a fiducial model for a generic hot Jupiter and find that the temperature and wind structures agree well with previously published models. We also show vertical velocities, radiative fluxes, and heating/cooling rates as a function of pressure and region of the atmosphere. Our radiative scheme allows us to self-consistently predict the infrared flux emitted by the planet, and we map this as a function of location on the globe, reproducing the standard shift of the brightest region eastward of the substellar point.

We show temperature profiles for models that (1) did not use the flux diffusion scheme or (2) assumed zero flux from the interior. In both cases the deepest pressure levels are mostly isothermal, instead of increasing in temperature toward the inner convective zone, significantly changing the ability of the atmosphere to exchange energy and momentum with the interior. However, the upper atmosphere (including and above the infrared photosphere) is relatively insensitive to the behavior in the deepest levels and has similar observable properties to our fiducial model.

We present models that use our new radiative transfer scheme in combination with a simplified form of magnetic drag, as an improvement on the models presented in Perna et al. (2010a). As in those earlier models, we find that as the drag strength is increased, the structure of the atmosphere changes; the strongest level of drag is able to significantly alter the atmosphere from its drag-free state. As a new result, we are able to measure the amount by which the brightest region of the atmosphere is offset from the substellar point, as a function of the strength of the magnetic drag used in each model. The weak- and medium-drag models have orbital thermal phase curves similar to the drag-free model. The winds are slowed enough in the strong-drag model that the brightest region of the atmosphere is closely aligned with the substellar point and, related to this effect, the day–night flux amplitude in the strong-drag model is slightly larger than for the drag-free model. We estimate the amount of ohmic heating that would be produced for each of our models and find that the medium-drag model should have a significantly inflated planetary radius, while the amount of heating in the strong-drag model is so high that it could perhaps lead to planet evaporation. A better analysis of magnetic effects on hot Jupiter atmospheres must wait for improvements in our code, however, using a more realistic form for the drag and explicitly including the effect of ohmic heating.

Finally, we estimate the rate of numerical loss of kinetic energy in each of our models and find it to be at the level of ∼15%, except for our model with strong magnetic drag, which we calculate to have losses consistent with zero at the percent level. Our strong-drag model has much slower winds than the other models presented here, and we speculate that strong numerical loss is directly related to those high-speed winds.

We thank the anonymous referee for comments that improved the quality of this paper. This work was performed in part under contract with the California Institute of Technology (Caltech) funded by NASA through the Sagan Fellowship Program. K.M. was supported by NASA grant PATM NNX11AD65G.