THE EFFECTS OF IRRADIATION ON HOT JOVIAN ATMOSPHERES: HEAT REDISTRIBUTION AND ENERGY DISSIPATION

A concise list of the parameters and their adopted values are described in Table 1. Most of the parameters are chosen to have typical values, e.g., the surface gravity of the hot Jupiter is g_p = 10 m s⁻². Since our primary interest is to explore the dependence of selected exoplanetary properties on the strength of irradiation and the shortwave atmospheric opacity, we keep the rest of the parameter values to be the same for all of the models (Table 1).

Table 1. Table of Common Parameters and Their Values

Symbol	Description	Units	Value
			(s)
$\tau _{\rm S_0}$	Normalization for optical depth of shortwave absorbers	...	5 × 103, 5 × 104
$\tau _{\rm L_0}$	Normalization for optical depth of longwave absorbers	...	104
	Enhancement factor due to collision-induced absorption	...	1
cP	Specific heat capacity at constant pressure (of the atmosphere)	J K−1 kg−1	14550.4
${\cal R}$	Specific gas constant (of the atmosphere)	J K−1 kg−1	4157.25
$\kappa \equiv {\cal R}/c_P$	Adiabatic coefficienta	J K−1 kg−1	2/7
P0	Reference pressure at bottom of simulation domain	bar	1000
gp	Acceleration due to gravity	m s−2	10
Rp	Radius of hot Jupiter	km	7.1 × 104
tν	Hyperviscous time	dayb	10−6
Nv	Vertical resolution	...	36

Notes. ^aWe term κ the "adiabatic coefficient," following Pierrehumbert (2010), and avoid calling it the "adiabatic index" in order to not confuse it with the ratio of specific heat capacities. ^bExpressed in terms of an exoplanetary day (i.e., rotational period).

Download table as:  ASCIITypeset image

Assuming a Bond albedo of zero, the irradiation⁵ temperature is defined as (Heng et al. 2012)

$$\begin{equation} T_{\rm irr} = T_\star \left(\frac{R_\star }{a} \right)^{1/2}, \end{equation} \tag{ 1 }$$

such that the irradiating flux is ${\cal F}_0 = \sigma _{\rm SB} T^4_{\rm irr}$. (For parameter values appropriate to Earth, one recovers ${\cal F}_0 \approx 1370$ W m⁻², the solar constant.) The irradiating flux ${\cal F}_0$, stellar radii R_⋆, effective stellar temperature T_⋆ and exoplanet–star spatial separation a are related by the expression

$$\begin{eqnarray} a &=& R_\star T^2_\star \left(\frac{\sigma _{\rm SB}}{{\cal F}_0} \right)^{1/2} \approx 0.09 \mbox{ AU} \left(\frac{R_\star }{R_\odot } \right) \nonumber \\ &&\times \left(\frac{T_\star }{6000 \mbox{ K}} \right)^2 \left(\frac{{\cal F}_0}{2 \times 10^8 \mbox{ erg cm}^{-2} \mbox{ s}^{-1}} \right)^{-1/2}, \end{eqnarray} \tag{ 2 }$$

where σ_SB is the Stefan–Boltzmann constant and we have assumed that the exoplanet orbits a Sun-like star. The use of the irradiation, rather than the equilibrium temperature, circumvents the issue of dealing with (confusing) factors of order unity related to assumptions about the efficiency of heat redistribution.

The tidal locking time is (Goldreich & Soter 1966; Bodenheimer et al. 2001)

$$\begin{equation} t_{\rm lock} \approx \frac{8Q}{45 \Omega _p} \left(\frac{\omega _p}{\Omega _p} \right) \left(\frac{M_p}{M_\star } \right) \left(\frac{a}{R_p} \right)^3, \end{equation} \tag{ 3 }$$

where ω_p and Ω_p are the rotational and orbital frequencies of the exoplanet, respectively. For our models, we get t_lock ∼ 10³–10⁸ yr (Q/10⁶)(ω_p/Ω_p). Thus, for a Gyr old solar-like star (which is what we assume), our assumption of tidal locking is plausible and self-consistent if ω_p/Ω_p ∼ 1. If ω_p/Ω_p ∼ 10, then this assumption becomes suspect only for Model C. With the rotational and orbital frequencies being equal, we obtain

$$\begin{equation} \Omega _p \approx \left(\frac{G M_\star }{a^{3}}\right)^{1/2}, \end{equation} \tag{ 4 }$$

where G is the gravitational constant and M_⋆ is the stellar mass. Thus, ${\cal F}_0$, Ω_p, and a need to be varied self-consistently. We note that setting the rotational and orbital frequencies to be equal results in a more constrained problem for our simulations. While the existence of close-in, non-spin-synchronized gas giants remains possible, there is little empirical evidence currently available to support such a modeling effort.

We present 16 models in this paper, divided into subgroups that we term "cold" (C, C1, and C2), "warm" (W, W1, and W2), and "hot" (H and H1). For Model W, the irradiating flux ${\cal F}_0$ is chosen to match the threshold value found by Demory & Seager (2011), below which transiting Kepler giant exoplanet candidates appear to have non-inflated radii. Each of these models is studied both with and without the presence of a temperature inversion. The parameters specific to each model are reported in Table 2.

Table 2. Table of Parameters, and their Values, Specific to Each Model

Model	${\cal F}_0$	Ωp	a	Tinit	Δt	Simulation Time	Cint	t−1ν
	(W m−2)	(s−1)	(AU)	(K)	(s)	(Earth days)	(J K−1 m−1)	(s−1)
C	2 × 104	1.3 × 10−6	0.28	589, 501	600	4500	107	2.1 × 10−1
C1	5 × 104	2.6 × 10−6	0.17	741, 630	240	3000	107	4.1 × 10−1
C2	7 × 104	3.4 × 10−6	0.15	806, 686	240	3000	107	5.4 × 10−1
W	2 × 105	7.5 × 10−6	0.09	1048, 892	180	2000	106	1.2
W1	5 × 105	1.5 × 10−5	0.06	1317, 1121	120	1500	106	2.4
W2	7 × 105	1.9 × 10−5	0.05	1433, 1219	120	1500	106	3.0
H	2 × 106	4.2 × 10−5	0.03	1863, 1585	120	1500	106	6.7
H1	5 × 106	8.4 × 10−5	0.02	2343, 1994	80	1000	106	13.3

Note. Initial temperatures quoted are for γ₀ = 0.5 and 2.0, respectively.

Download table as:  ASCIITypeset image

The adiabatic coefficient κ is related to the thermodynamical properties of the atmospheric gas (Pierrehumbert 2010),

$$\begin{equation} \kappa \equiv \frac{{\cal R}}{c_P} = \frac{2}{2+n_{\rm dof}}. \end{equation} \tag{ 5 }$$

By assuming the atmosphere to be dominated by molecular hydrogen (i.e., the number of degrees of freedom of the gas is n_dof = 5), we get κ = 2/7. Furthermore, the specific gas constant is related to the mean molecular weight μ and the universal gas constant ${\cal R}^\ast$ by ${\cal R} = {\cal R}^\ast /\mu = 4157.25$ J K⁻¹ kg⁻¹. We use μ = 2 to emphasize the detachment from any specific case study, and remark that our results are rather insensitive to small changes in the mean molecular weight (e.g., μ = 2.35 as is typical for a solar mix). It follows that the specific heat at constant pressure is c_P = 14550.4 J K⁻¹ kg⁻¹. We recall that the adiabatic index is γ = 1 + 2/n_dof = 7/5 for an atmosphere dominated by molecular hydrogen and that the sound speed is $c_s = \sqrt{\gamma k_{\rm B} T/ 2 m_{\rm H}}$ with k_B denoting the Boltzmann constant, T the temperature, and m_H the mass of a hydrogen atom.

In this work, we adopt the dual-band approximation for radiative transfer, namely, that the blackbody peaks of the stellar and exoplanetary emission are well separated in wavelength or frequency (see Heng et al. 2011a and references therein). One then requires two broadband opacities to completely describe the radiative transfer: a shortwave/optical opacity κ_S to quantify the absorption of starlight with depth/pressure and a longwave/infrared opacity κ_L to quantify the absorption of thermal emission from the exoplanet. We neglect the effects of scattering (see Heng 2012). Varying κ_S while holding κ_L fixed mimics the effect of extra, shortwave absorbers in the atmosphere, of unspecified chemistry, which cause a temperature inversion when κ_S > κ_L (Hubeny et al. 2003; Burrows et al. 2007b; Burrows & Orton 2010; Knutson et al. 2008).

In the absence of empirical constraints, we assume the shortwave opacity to be constant, which means that it is related to the shortwave optical depth in a simple manner

$$\begin{equation} \tau _{\rm S} = \frac{\kappa _{\rm S} P}{g_p} = \frac{\tau _{\rm S_0} P}{P_0}, \end{equation} \tag{ 6 }$$

where P₀ = 1 kbar is the reference pressure at the bottom of the simulation domain. The longwave optical depth is described by two terms (Heng et al. 2011a),

$$\begin{equation} \tau _{\rm L} = \tau _{\rm L_0} \left[ \frac{P}{P_0} + \left(\epsilon -1 \right) \left(\frac{P}{P_0} \right)^2 \right], \end{equation} \tag{ 7 }$$

where the second term is an approximate way to account for the effects of collision-induced absorption, which introduces an enhancement factor of at the bottom of the simulation domain. The longwave optical depth normalization is $\tau _{\rm L_0} = \kappa _0 P_0/g_p$, where κ₀ is the longwave opacity normalization. We pick = 1. Ignoring collision-induced absorption in the infrared is a good approximation as long as this effect is only significant at atmospheric layers residing below the longwave photosphere (Heng et al. 2011a).

In the absence of clouds/hazes, the ratio γ₀ ≡ κ_S/κ₀ determines if a temperature inversion exists in the atmosphere (Hubeny et al. 2003; see also the discussion in Heng et al. 2012). If γ₀ > 1, the shortwave photosphere sits above the longwave photosphere, and an inversion exists; if γ₀ < 1, the relative locations of the two photospheres are reversed, and there is no inversion. To study the effect of a temperature inversion on the day–night heat redistribution, we examine models with γ₀ = 0.5 and 2. We choose a typical value for the longwave opacity normalization: κ₀ = 0.01 cm² g⁻¹ such that the longwave/infrared photosphere resides at

$$\begin{equation} P_{\rm IR} \sim \frac{g_p}{\kappa _0} = 0.1 \mbox{ bar}. \end{equation} \tag{ 8 }$$

The shortwave opacity follows from the chosen value of γ₀.

The global-mean temperature–pressure profiles are shown in Figure 1 for the two chosen values of γ₀ = 0.5 and 2, and the three representative flux strengths indicated with Models C, W, and H in Table 2. The effect of the temperature inversion is clearly evident: models with γ₀ = 0.5 are hotter in the interior and cooler in the outer parts, while models with γ₀ = 2 display the opposite behavior. We have set the internal heat flux of our hot Jupiter models to be zero. For completeness, we include in the Appendix the zonal-mean zonal wind, temperature, potential temperature, and stream function profiles corresponding to these six models.

Zoom In Zoom Out Reset image size

The Princeton-GFDL FMS dynamical core solves the primitive equations of meteorology (see Heng et al. 2011b and references therein). In the vertical, we choose the simulation domain to span six orders of magnitude in pressure, from P = 1 mbar to 1 kbar; the latter value is chosen to roughly match the transition from the radiative to the convective zone of a hot Jupiter. The vertical resolution is N_v = 36, such that each pressure scale height is covered by about three grid points. The horizontal resolution is 192 × 96 (longitude versus latitude). In the absence of a physically motivated scheme for more sophisticated initial conditions, the simulations are started from a state of rest with a constant temperature T_init.⁶ For = 1, the initial temperature is (Heng et al. 2011a)

$$\begin{equation} T_{\rm init} \approx \left[ \frac{{\cal F}_0}{16 \sigma _{\rm SB}} \left(2 + \frac{\sqrt{3}}{\gamma _0} \right) \right]^{1/4}. \end{equation} \tag{ 9 }$$

The specific value of T_init used to initiate each of the models is reported in Table 2.

The areal heat capacity of the simulation bottom (C_int) is chosen to have a value such that the thermal inertia of gaseous hydrogen matches the radiative time constant at P = 1 kbar (Heng et al. 2011a). For Models C, W, and H, it suffices to adopt C_int = 10⁷, 10⁶, and 10⁶ J K⁻¹ m⁻², respectively, as it has been shown by Heng et al. (2011a) that order-of-magnitude variations in C_int produce small (<1 K) variations in the temperature–pressure profile.

As already noted by Heng et al. (2011b), numerical noise accumulates at the grid scale and has to be dissipated via a "hyperviscosity" in our spectral code. Since the hyperviscous term is neither an original nor a physical part of the primitive equations, its value cannot be specified from first principles. As in Heng et al. (2011a, 2011b), we pick the order-of-magnitude value of the hyperviscous timescale to be as large as possible in order for the simulations to collectively run to completion: t_ν = 2π × 10⁻⁶/Ω_p. Decreasing the hyperviscous timescale—and hence increasing the rate of dissipation (t⁻¹_ν)—by one to two orders of magnitude leaves the hotspot offset largely unchanged (see Figure 19 of Heng et al. 2011a), but decreases the zonal wind speeds by ∼10% (Heng et al. 2011b). Our estimates for the zonal wind speeds are thus upper limits with regards to this numerical issue.

Finally, we note that convective adjustment is applied when the dry lapse rate becomes super-adiabatic, i.e., the Schwarzschild criterion for convective stability is violated. Our convective adjustment scheme is described in Heng et al. (2011a).

Our three-dimensional circulation models, while not directly simulating the onset and development of shocks and shear instabilities, allow us to establish whether the flow is subject to such instabilities and to shocks and to determine the size of the affected region. In the current study, we are primarily interested in exploring the effects of irradiation, while also examining the influence of (shortwave) opacity.

4.1.1. Dissipation from Kelvin–Helmholtz-unstable Gas

The first form of hydrodynamic dissipation that we consider is due to the KH instability in the thermally stratified shear flow. Whether such a flow becomes unstable can be estimated by measurements of the Richardson number

$$\begin{equation} {\cal R}i = N^2 \left(\frac{\partial v_\Theta }{\partial z} \right)^{-2} = -{\cal R} \sigma _0^{\kappa -1} \frac{\partial \theta _T}{\partial \sigma _0} \left(\frac{\partial v_\Theta }{\partial \sigma _0} \right)^{-2}, \end{equation} \tag{ 14 }$$

where

$$\begin{equation} N = \left(\frac{g_p}{\theta _T} \frac{\partial \theta _T}{\partial z} \right)^{1/2} \end{equation} \tag{ 15 }$$

is the Brunt–Väisälä frequency, θ_T ≡ T(P/P₀)^−κ is the potential temperature, and σ₀ ≡ P/P₀. The criterion for the flow to become KH unstable is expressed by the condition ${\cal R}i<1/4$ (Kundu & Cohen 2004; Li & Goodman 2010).

We computed the Richardson number at each point of the flow. For Model C, we find ${\cal R}i>1/4$ everywhere, and hence shear instabilities are not expected at low levels of irradiation. In the warm and hot models, we find ${\cal R}i<1/4$ only in the outermost pressure levels, down to a few millibars. Hence, we conclude from this analysis that the onset and consequent dissipation of shear instabilities is unlikely to play a major role in inflating hot Jupiters.

4.1.2. Dissipation from Shocked Gas

In order to estimate the presence and amount of dissipation due to shocked gas, we begin by keeping track of the Mach number at each location within the flow; since zonal velocities are the most significant ones, we define the Mach number as

$$\begin{equation} {\cal M}(\Theta,\Phi,r)\;=\; \frac{v_\Theta }{c_s}\,. \end{equation} \tag{ 16 }$$

Shocks are likely to occur for ${\cal M}\gtrsim 1$, barring special conditions that allow the pre- and post-shock regions to remain in causal contact (Zel'dovich & Raizer 1966). One may regard our simple estimates of the shock dissipation as upper limits since ${\cal M}\gtrsim 1$ is a necessary but insufficient condition for the development of shocks.

At each pressure level within the atmosphere, flow velocities change with both latitude and longitude. To investigate how deep shocks can penetrate, in Figure 4 we show the maximum Mach number as a function of atmospheric depth for the three representative Models C, W, and H, each with and without temperature inversions. In the coldest model, the flow is always subsonic. Wind speeds increase, as expected, with the strength of the irradiating flux and hence with the day–night temperature gradients. In the W model, the flow can become supersonic (with Mach numbers of ∼1.5–2) at pressure levels above a few tens of bars. In the H model, shock penetration reaches a few bars, with Mach numbers ∼2.5–3. For the same strength of irradiation, flows with temperature inversions reach a higher Mach number, due to the higher temperatures attained in the outermost layers.

Zoom In Zoom Out Reset image size

In order to make a crude estimate of the amount of energy available to be dissipated in shocks, we first compute the total amount of kinetic energy available in the flow up to pressure P, $E_{\rm kin}(P) =(1/2)\int _{r(P_0)}^{r(P)} dV \rho (v^2_\Theta + v^2_\Phi)$. Note that the dominant term comes from the zonal component ($v_\Theta$). Figure 5 shows E_kin(P) in the three representative Models C, W, and H. The cumulative amount from the whole flow varies by over three orders of magnitude between the C and the H model, with a weaker dependence on the shortwave opacity. Generally, we find that the stronger the irradiation, the deeper the penetration of the zonal winds.

Zoom In Zoom Out Reset image size

Together with the local kinetic energy density, at each point of the flow we track both the Mach number and the Richardson number. If either ${\cal M}>1$ or ${\cal R}i<1/4$, we assume that a fraction of 0.5 of the kinetic energy is dissipated/converted into heat. In strong shocks, this fraction is given by 4(γ − 1)/(γ + 1)², where γ is the adiabatic index (Draine & McKee 1993); for γ = 7/3, this amounts to 0.48. For shear instabilities, Li & Goodman (2010) estimate the dissipated fraction of kinetic energy to be ≈0.5–1. Hence, we take $E_{\rm dis}(P) = (1/2) f_{\rm dis}\int _{r(P_0)}^{r(P)} dV \rho (v^2_\Theta + v^2_\Phi)$, where f_dis = 0.5 if ${\cal M}>1$ or ${\cal R}i<1/4$ and f_dis = 0 otherwise. We remark once again that this is a rather crude estimate, motivated by the fact that current global three-dimensional circulation models are not yet able to follow the development, growth, and dissipation of instabilities and shocks.

The cumulative dissipated energy, E_dis(P), is shown in Figure 5 for the C, W, and H models. In the C model, E_dis(P) = 0 since, as discussed above, the criterion for the onset of the KH instability is never satisfied, and the velocities are always subsonic. In the W and H models, on the other hand, a fraction of about 10%–20% of the kinetic energy is available for dissipation, with the fraction being generally larger in models with temperature inversions, due to the stronger velocities (and their gradients) in the upper atmospheric layers, which are the ones more prone to becoming either KH unstable or shocked. For a given value of γ₀, dissipation occurs down to deeper levels as the strength of the irradiating flux increases, reflecting the penetration depth of the shocks (cf. Figure 4). However, even in the hottest models, shock dissipation never penetrates deeper than a few bars.

Previous work (Batygin & Stevenson 2010; Perna et al. 2010b) has shown that for Jupiter-like magnetic field strengths Ohmic dissipation can be a significant source of energy for the exoplanet. In particular, using a simulation specialized to the physical parameters and conditions of HD 209458b, Perna et al. (2010b) showed that Ohmic dissipation can penetrate down to levels deep enough to influence the radius evolution of the exoplanet. Since irradiation is the driver of the flow dynamics, currents (and hence Ohmic dissipation) are expected to depend sensitively on its strength. In the following, we investigate this dependence using our 16 baseline models. We assume the B-field to be an axisymmetric dipole where the rotational and dipolar axes are aligned. The dependence of Ohmic dissipation on the magnitude of the B-field has been explored (inclusive of an approximate feedback from magnetic drag), by Perna et al. (2010b).

To compute Ohmic dissipation within the context of our simulations, we follow the formalism of Liu et al. (2008) and Perna et al. (2010b). Under the assumption that zonal winds are stronger than the meridional motions, the only relevant component of the induction equation is the toroidal one. Under steady-state conditions, it was shown by Liu et al. (2008) that the resulting dominant component of the current induced by the zonal flow is the meridional one, which we write below in terms of the latitude Φ and the longitude Θ, and also incorporating the expression for an aligned dipole for the magnetic field:

$$\begin{eqnarray} {\cal J}_\Phi &=& \frac{B_0 r_{\rm max}^3 \sigma _e}{(r_{\rm max} - \tilde{z}) c} \int ^{\tilde{z}}_0\hspace{-3.61371pt}dz^\prime \left[\frac{2 \sin \Theta }{(r_{\rm max} - z^\prime)^2} \frac{\partial v_\Theta }{\partial z^\prime } + \frac{2 v_\Theta \sin \Phi }{(r_{\rm max} - z^\prime)^3}\right. \nonumber \\ &&+ \left.\frac{\cos \Phi }{(r_{\rm max} - z^\prime)^3} \frac{\partial v_\Theta }{\partial \Phi } + \frac{v_\Theta \sin \Phi }{(r_{\rm max} - z^\prime)^3}\right]. \end{eqnarray} \tag{ 17 }$$

Here, η = c²/(4πσ_e) is the electrical resistivity, and the electrical conductivity σ_e is computed using the same formalism as in Perna et al. (2010a, 2010b). The familiar radial coordinate, denoted by r, is equal to z + R_p where z is the vertical height measured from the bottom of the simulation domain. The top of the simulation domain is located at r_max = max(z + R_p). We choose to compute the current density as a function of the variable $\tilde{z}$, which is defined such that $\tilde{z}=0$ when r = r_max, i.e., it is the vertical height measured from the top of the model atmosphere. Again, $v_\Theta$ is the zonal velocity. The quantity B₀ is the normalization term in the dipole magnetic field, such that it is B₀ at the equator and 2B₀ at the poles. For our calculations, we choose B₀ = 1 G. Note that in the above equation we have neglected the contribution from an unknown boundary current.

The total Ohmic power dissipated between the bottom shell at pressure P₀ and a shell of material located at pressure P, is then readily computed as

$$\begin{eqnarray} &&P_J(P)=\int _{r(P)}^{r(P_0)}\hspace{-3.61371pt}dr^{\prime } {r^{\prime }}^2 \int _0^{2\pi }\hspace{-3.61371pt} d\Theta \int _{-\frac{\pi }{2}}^{\frac{\pi }{2}}\hspace{-3.61371pt} d\Phi \cos \Phi \frac{[{{\cal J}_\Phi (r^{\prime },\Theta,\Phi)}]^2}{\sigma _e(r^{\prime },\Theta,\Phi)}. \nonumber\\ \end{eqnarray} \tag{ 18 }$$

Figure 6 shows the cumulative Ohmic power as a function of the atmospheric depth, for the three representative Models C, W, and H, both with and without temperature inversions. The magnitude of P_J has a very strong dependence on the irradiating flux, varying at maximum by about 16–17 orders of magnitude between the C and the H models. While the increasing zonal wind speeds with T_irr do contribute to this effect, the dominant factor is the exponential dependence of the electrical conductivity on temperature.

Zoom In Zoom Out Reset image size

For the same T_irr, the maximum Ohmic power is larger for γ₀ = 2 (with a temperature inversion) than it is for γ₀ = 0.5 (no inversion), while the penetration depth on the other hand is larger for γ₀ = 0.5. This behavior can be understood as follows. For γ₀ = 2, the outer layers of the atmosphere are hotter than the inner layers, and hence the conductivity is much larger in the outer atmospheric regions. Since velocity gradients are also stronger in these regions, the net result is the generation of very large peripheral currents and hence large dissipation, which however saturates, even in the hottest models, close to the photosphere. For γ = 0.5 on the other hand, the outer layers are colder, and the low conductivity results in weak currents. Deeper in, the temperature rises, and so does the conductivity, resulting in an increase of the Ohmic power at deeper atmospheric levels. However, the maximum dissipated power saturates at smaller values than for the case of γ₀ = 2, since for γ₀ = 0.5 the region with higher conductivity corresponds to zonal flows with weaker gradients.

A comparison between the total power (i.e., integrated over the entire atmosphere) generated through shocks and shear instabilities, and through Ohmic dissipation, is shown in Figure 7. Since our simulations do not explicitly follow shocks and the development of shear instabilities, they cannot make predictions for the timescales over which dissipation associated with these phenomena occur. To compute the hydrodynamic power in Figure 7, we assumed that E_dis is dissipated over a rotational timescale (which is comparable in magnitude to the advection timescale). The two-dimensional simulations by Li & Goodman (2010) show that once shocks are formed they slow down over a very short time (shorter than the advection time). As long as the dissipation timescale is shorter than 2π/Ω_p, the power shown in Figure 7 can be considered as a good proxy for its corresponding time-averaged quantity over a rotational cycle.⁷

Zoom In Zoom Out Reset image size

Our simulations show that the Ohmic power is a much steeper function of irradiation than the hydrodynamic power, which we find to be dominated by shocks. For the latter, only velocities matter, while the dissipated Ohmic power depends on velocity gradients. But most importantly, as the temperature increases, so does the ionization fraction, and hence the conductivity, which results in stronger currents. However, we remark once again that the slope at high T_irr would be shallower if magnetic drag were included self-consistently in the simulations. Specifically, Menou (2012b) has shown via simple scaling laws that magnetic drag rapidly retards the atmospheric winds for strong stellar irradiation, and that, for a given magnetic field strength, there is a value of the irradiation/equilibrium temperature for which Ohmic dissipation peaks (and then declines as irradiation becomes more intense). Our results for the dissipated power should therefore be considered as upper limits.

Finally, we note that Heng (2012) provides a companion study to our paper by using semi-analytical models to elucidate the interplay between atmospheric scattering, absorption, and Ohmic dissipation with relevance to the possible presence of atmospheric clouds or hazes.

As discussed in Section 1, the anomalously large radii of a fraction of hot Jupiters remain a persistent problem for researchers. An "extra" source of internal energy is needed for the radii to be larger than predicted by standard evolutionary models. Guillot & Showman (2002) argue that if about 1% of the stellar insolation flux were deposited at pressures of tens of bars deep in the atmosphere, cooling would slow down sufficiently to explain the inflated radii of hot Jupiters. The evolutionary model presented by Guillot & Showman (2002) was tailored to the specific case of the exoplanet HD 209458b, with a radius of R = 1.35R_J (R_J is the Jupiter radius), and a convective–radiative boundary located at pressures of about 1 kbar (at its current age). Note that for exoplanets with comparable mass, the boundary layer between the convective/radiative zones is located at lower pressures as the planet is younger and larger (e.g., Guillot & Showman 2002; Burrows et al. 2003). The amount of internal heat needed to inflate the radius by a certain ΔR is also a strong function of the planet mass (Miller et al. 2009). The smaller the mass, the easier it is to inflate the planet. For example, Miller et al. (2009) show that a 20% increase in radius (with respect to R_J) can be achieved with a mere 10²⁴ erg s⁻¹ in a planet of mass 0.1 M_J, while the same radius increment would require about 10²⁸ erg s⁻¹ for a 10 M_J planet.

While several mechanisms could transport and dissipate energy into the interior, such as eccentricity damping (Bodenheimer et al. 2001), or forced turbulent mixing in the radiative layer (Youdin & Mitchell 2010), here we have investigated energy deposition resulting from the development of shear instabilities and shocks and from Ohmic dissipation. For all of these phenomena, our numerical simulations, while not modeling the dissipation itself and its feedback on the flow, do allow us to estimate their importance from an energetic point of view, as a function of the strength of the irradiating flux and of the atmospheric opacity.

We have found that shear instabilities can develop only very superficially, even in the hottest models, and hence they are unlikely to meaningfully influence the energy budget of the exoplanet. If the source of weak turbulence postulated by Youdin & Mitchell (2010) results from shear instabilities, then their proposed mechanism of forced turbulent mixing does not operate deeply enough. To further test the relevance of their mechanism in global circulation models requires the specification of an instability and its criterion for being invoked.

Shocks penetrate deeper, up to a few bars in the hottest models. This is a depth which could be interesting for influencing the thermal evolution of planets, and hence the degree by which their radii are able to shrink during the evolution. For the particular case of HD 209458b, the evolutionary models of Guillot & Showman (2002) showed that even at a pressure as shallow as 5 bars, dissipation of 10% of the absorbed stellar flux (corresponding to $\dot{E}=2.4\times 10^{28}$ erg s⁻¹ for this planet) is sufficient to retard cooling so that the planet contracts to its current size on a timescale ∼5 Gyr, comparable to its age. Inspection of the trend for the dissipated hydrodynamic power versus irradiating flux in Figure 7 shows that, for the strength of ${\cal F}_0$ corresponding to the case of HD 209458b, the hydrodynamic power is on the order of a few × 10²⁸ erg s⁻¹. Therefore, since the dissipated hydrodynamic power that our models predict falls in an interesting order-of-magnitude range (at least for the most irradiated objects), and its penetration depth also reaches levels close to the interesting ones identified by evolutionary models,⁸ we suggest that shock dissipation could indeed play a significant role in the thermal evolution of the planets, and hence influence their radii, although especially so in the most irradiated planets.

Our investigation of Ohmic dissipation has showed that for a B-field strength on the order of a Gauss, the total dissipated Ohmic power is generally lower than the hydrodynamic power except for the most irradiated hot Jupiters. However, it penetrates deeper, down to pressures of several tens of bars, if the atmosphere has no temperature inversion. In the evolutionary scenario of Guillot & Showman (2002), a deeper penetration depth can be "traded" for a smaller power. For example, in the case of HD 209458b, their evolutionary models could reproduce the exoplanetary radius at its current age both with 10% of the irradiating flux dissipated at 5 bars, or with 1% (i.e., 2.4 × 10²⁷ erg s⁻¹) dissipated at 21 bars. For a B-field strength of a few Gauss, this power is available in the most irradiated objects. Clearly, the importance of Ohmic dissipation for the thermal evolution of planets also depends strongly on the magnitude of the B-field, since the Ohmic power scales as⁹ B²₀, and, for fixed planet conditions, the penetration depth also increases with larger B₀ (Perna et al. 2010b). However, for a given B-field, since both the penetration depth and the power dissipated at a fixed pressure are rapidly increasing functions of the strength of the irradiating flux, the influence of Ohmic dissipation on radius inflation is expected to depend strongly on stellar irradiation.

Figure 8 summarizes the main observational constraints on the radius inflation and atmospheric circulation of hot gas giant exoplanets. The exoplanetary radius measured for transiting gas giants (M > 0.3 M_J) is shown as a function of the irradiation temperature. Different symbols indicate more/less massive exoplanets, exoplanets with or without a temperature inversion measured near the infrared photosphere, and exoplanets with or without efficient day–night temperature redistribution. The presence of a temperature inversion is inferred from the broadband infrared emission spectra measured during secondary eclipse with the Spitzer satellite. The efficiency of the day–night temperature redistribution is inferred from the comparison of the effective temperature in the Spitzer passbands with the equilibrium temperature. If the measured temperature is much larger than the equilibrium temperature, the redistribution of heat to the night side is likely to be low.¹⁰ In one case, HD 189733b, the heat redistribution is measured more directly, from the infrared phase curve of the exoplanet at 8 and 24 μm (Knutson et al. 2007, 2009).

Zoom In Zoom Out Reset image size

The exoplanet radii and irradiation data have been taken from the edited compilation at http://exoplanets.org (see references therein for primary sources). Information on the absence/presence of temperature inversions is adopted from Knutson et al. (2011 see references therein), while information on redistribution is taken from the compilation of Harrington (2011), updated with more recent data from Spitzer for HD 80606, XO-3, HAT-P-7, and WASP-12; we also include inferences on XO-4, HAT-P-6, and HAT-P-8 (Todorov et al. 2012). It is important to note that inferences on redistribution (presented in the lower panel of Figure 8) are based on data that stretches the capacity of the Spitzer instruments to their limits, and that the uncertainties associated with instrumental effects are still high. In at least two cases, Ups And b and HD 149026b, further measurements have shown the initial estimates for the surface temperature and redistribution efficiency to be incorrect. These two objects are not in our plot, Ups And b because it is not a transiting planet and therefore its radius is not measured, and HD 149026b because it is a core-dominated Saturn rather than a hot Jupiter, and hence its radius R = 0.65 ∼ R_J, is not determined by its thermal history, but by its heavy core.

Information on the angular offset of the hotspot is only sparsely available. Spitzer observations of HD 189733b show that the hottest region of the thermal phase curve is advected to the east of the substellar point (Knutson et al 2007) by 16° ± 6°; Knutson et al. (2009) inferred the one-dimensional brightness maps at 8 μm and 24 μm and found the offsets to be 30° ± 4° and 23° ± 7°, respectively. Measurements for a couple of other objects are less certain. Crossfield at al. (2010) report for Ups And b, which has a good redistribution, an offset of the one-dimensional brightness map to the east of either 57° ± 21° or 845 ± 23, using two different data sets for the analysis. Cowan et al. (2012), report, for WASP-12b, an offset to the east of 12° ± 6° from the Spitzer phase curve at 4.5 μm, and an offset of 53° ± 7° at 3.6 μm.

Figure 8 illustrates the strong correlation between irradiation and radius inflation. This indication that injection of a fraction of the irradiation into the exoplanet as internal entropy is responsible for the radius inflation (as first suggested by Showman & Guillot 2002) is reinforced by the observed correlation between radius inflation and the mass of the exoplanet (Fortney et al. 2011), and by the absence of radius inflation for cooler Jupiter-sized exoplanet candidates from the Kepler mission (Demory & Seager 2011). Hence, the observations strongly support explanations of radius inflation as a consequence of high stellar irradiation. The observed slope of the dependence, though, seems even steeper than that predicted by the injection of a constant fraction of the stellar flux as internal entropy as proposed by Showman & Guillot (2002), suggesting that the efficiency of the coupling mechanism increases with temperature. There is no significant correlation between radius inflation and orbital eccentricity (Pont 2009), nor with temperature inversion.

For the dissipation mechanisms that we have explored in this work, both of hydrodynamic and of magnetohydrodynamic origin, our simulations show a strong correlation with the strength of the irradiating flux. In particular, for a B-field on the order of a Gauss, we found that, for T_irr ≲ 2000 K, the dissipated Ohmic power is too small (and its penetration depth too shallow) to affect the radius evolution in any appreciable way. However, as T_irr approaches ∼2000 K, the Ohmic power grows noticeably, and its penetration depth, reaching down to several tens of bars, can then affect the evolution of the exoplanetary radius (Guillot & Showman 2002). Interestingly, as Figure 8 shows, radius inflation starts to become important as the irradiation temperature approaches the 2000 K range. We also note that variations of the magnetic field from exoplanet to exoplanet can result in a substantial scatter in the radius evolution, since the Ohmic power scales as B²₀, and this rapid growth with magnetic field strength dominates over the reduction to the power caused by the feedback effect of magnetic drag, at least for an interesting range of B-field strengths (Perna et al. 2010b). For a Jupiter-like magnetic field (∼14 G at the pole), the power produced by Ohmic dissipation can already affect planets with T_irr ∼ 1700–1800 K. The atmospheric (shortwave) opacity further contributes to the scatter, since the penetration depth of Ohmic dissipation is a non-negligible function of this variable.

In summary, and at a qualitative level, our simulations predict that both hydrodynamic dissipation (dominated by shocks) and Ohmic dissipation (for magnetic field strengths at least as large as several Gauss) can play an important role in altering the thermal evolution of hot Jupiters, slowing down cooling, and favoring the presence of planets of larger dimensions than predicted by evolutionary models without "extra" energy injection. A significant dissipated power (both hydrodynamic and Ohmic if B₀ ≳ 7–8 G) on the order of several percents of the stellar insolation, is predicted for irradiation temperatures T_irr ≳ 1700–1800 K. For both forms of dissipation, the stronger the irradiating flux, the larger the total dissipated power and the deeper the penetration depth, resulting in an increasingly larger inflation as the flux becomes more intense. At large insolation, even if the direct dissipation depth is rather shallow (as it is for shocks and for currents in highly inverted temperature profiles), vertical mixing can help to carry down heat to levels of tens of bars, hence maintaining the energy injection in the relevant regions for affecting the thermal evolution of the planet.

A quantitative prediction of the slope of the R_p(T_irr) curve would need a coupling between a specific dissipation mechanism and an evolutionary model for the exoplanet, which is beyond the scope of the present study. But even so, as discussed above, a large scatter would be expected due to exoplanets with different intrinsic properties (opacity, gravity, magnetic field strength, etc.).

Although the statistical information is still rather limited, inspection of Figure 8 further shows that the efficiency of day–night heat redistribution seems to diminish rapidly for T_irr above 2200 K, with the hottest gas giants showing dayside temperatures in the infrared much higher than the equilibrium temperature, an effect pointed out by Harrington (2011). Interestingly, our simulations (cf. Figure 2) also show that from irradiation temperatures around 2200–2500 K, the atmospheric flow pattern begins to evolve from a very good redistribution to a regime of no-redistribution. Exoplanets with temperature inversions evolve more quickly toward a no-redistribution regime by displaying, for the same irradiating flux, a larger day–night flux contrast. This is again supported by observations, as well as by other work (Fortney et al. 2008). Note that the core-dominated Saturn HD 149026b—not included in our observational plots since its radius is the result of a different evolutionary history than for the hot Jupiters—has a measured night–day flux contrast of 45% ± 19% (Knutson et al. 2009); this is consistent with our predictions (cfr. Figure 2) for redistribution at its temperature of T_irr ∼ 2500 K.

It may also be worthwhile to note that the drop in redistribution efficiency appears to occur around the temperature at which radius inflation saturates, although this is only a tentative indication at this point due to the low number of known cases, the large scatter at high temperatures, and the assumptions involved in inferring redistribution from infrared dayside fluxes in the Spitzer passbands. From a theoretical point of view, a saturation of radius growth at very high values of T_irr is expected due to feedback effects (e.g., shock dissipation and magnetic drag, both contributing to reduce the speeds of the fastest flows).

For the hotspot offset with respect to the substellar point, there is insufficient data statistics yet to compare theoretical and observational trends. However, we do note that as the observations show that the offset can be both rather small or very large, our theoretical modeling also shows a wide range of values, and further predicts a correlation between offset width and degree of redistribution. Although our modeling here has been for "typical" exoplanets and is not tailored to the properties of specific objects, it is, however, tempting to point out that with the irradiating flux level of HD 189733b, our predicted offset is in the 20°–40° range (cf. Figure 3), consistent with the range of values (16°–34°) obtained from the inferred one-dimensional brightness maps, at 8 and 24 μm, by Knutson et al. (2009). In the case of WASP-12b, the extremely high value of the irradiation temperature, T_irr ∼ 3500 K, is above the temperature range studied in our simulations. However, inspection of the offset trend at high T_irr in Figure 3 would suggest a value around a couple of tens of degrees at those high temperatures. The measured offset of 16° ± 4° at 4.5 μm appears compatible with this estimate. On the other hand, the much larger offset of 53° ± 7° at 3.6 μm does not. The difficulty in reconciling this large phase offset with the other properties of the object was already discussed by Cowan et al. (2012), though they also noted that the measurement could be affected by highly correlated residuals near the purported peak. If real, such a measurement might indicate a rather unusual opacity, and a detailed theoretical modeling would require a non-gray radiative transfer treatment (Fortney et al. 2006; Showman et al. 2009; Burrows et al. 2010). Similarly, for the case of Ups And b, the unusually large offset values reported by Crossfield et al. (2010), albeit with somewhat large error bars, cannot be simply accounted for with the opacity values adopted here. However, we note that even within the simple framework of the dual-band model, larger offsets can be produced for smaller values of both the optical and infrared opacities.