ON THE ANOMALOUS RADII OF THE TRANSITING EXTRASOLAR PLANETS

The disparate radii of the transiting extrasolar planets have defied a straightforward explanation. With nearly 100 well-characterized examples now known, it is clear that where Jovian planets are concerned, variations in mass and stellar insolation are responsible only for a fraction of the observed range in planetary sizes.

Radius anomalies emerged with the discovery of the first transiting extrasolar planet, HD 209458b (Charbonneau et al. 2000). Structural evolutionary models computed by, e.g., Bodenheimer et al. (2001) and Guillot & Showman (2002) suggested that an H–He dominated planet with HD 209458b's mass, insolation, and age should have a radius of order R ∼ 1.1 R_Jup. This figure is startlingly at odds with the observed value, R = 1.38 ± 0.02 R_Jup (Southworth 2010).

A decade of additional discoveries has made it clear that HD 209458b's radius anomaly is by no means anomalous, and the status of discussion in the field is well covered in the recent reviews by Showman et al. (2010), Burrows & Orton (2009), and Baraffe et al. (2010). The extraordinary variation in observed radii in the planetary population is indicated by Figure 1. This diagram charts the measured radii and uncertainties for the 90 transiting planets with accurately measured masses against the planets' orbit-averaged effective temperatures given by

$$\begin{equation} T_{\rm eff}=\left({R_{\star }\over {2a}}\right)^{1/2} {T_{\star } \over {(1-e^{2})^{1/8}}} \,, \end{equation} \tag{ 1 }$$

where R_⋆ is the stellar radius, a is the semimajor axis, e is the orbital eccentricity, T_⋆ is the stellar effective temperature, and zero albedos have been assumed. Only planets with 0.1 M_Jup < M_pl < 10.0 M_Jup were used.³

Zoom In Zoom Out Reset image size

It is evident that the observationally well-characterized transiting planets exhibit a wide range of radii for a given mass, and a number of models have been advanced to explain this dispersion. Planets that are smaller than expected, with HD 149026b (Sato et al. 2005) providing the canonical example, have generally had their small radii attributed to a large fraction of heavy elements in their interiors. For planets such as HD 209458b, which are larger (and sometimes much larger) than expected, explanatory models must generally appeal to a cryptic source of heating at adiabatic depth. A number of heating mechanisms have been proposed, including dissipation stemming from tidal orbital circularization (Bodenheimer et al. 2001), "kinetic heating" in which wind energy is converted into heat (Guillot & Showman 2002), and Ohmic dissipation (Batygin & Stevenson 2010). This Letter shows that current observational evidence indicates that the radius anomaly, the difference between the observed radius and that predicted by theoretical models, displays a well-defined correlation with temperature (from Equation (1)). Further, this temperature dependence can be interpreted to favor a magnetohydrodynamic mechanism in which Ohmic heating and attendant magnetic braking of the velocity field at sufficient atmospheric ionization fraction play significant roles. Our approach relies on scaling arguments, as we believe that the extant observations are not yet sufficient to permit the discrimination between detailed models containing multiple free parameters.

Naively, one might expect that the radius of a mature gas-giant planet is primarily determined by its mass and by the amount of radiative energy that it receives from its star. This conjecture can be tested by evaluating model radii, $R_{m_{i}}$, of "baseline" evolutionary models for H–He composition planets spanning a range of masses and insolation, and comparing with those of known corresponding transiting planets, $R_{o_{i}}$. If the model has explanatory power for the aggregate of N known planets, then it should produce a statistically significant decrease in the quantity,

$$\begin{equation} \chi ^{2}_{m}={1\over {(N-N_{f})}}{\displaystyle \sum \limits _{i=1}^{N} {({R_{m_{i}}-R_{o_{i}})^{2}}\over {\sigma _{i}^{2}}} }, \end{equation} \tag{ 2 }$$

in comparison to χ²_null, obtained by replacing $R_{m_{i}}$ with R_av = 1.2 R_Jup, the average observed radius for transiting planets having 0.1 M_Jup < M_pl < 10.0 M_Jup. In the above equation, N_f = 2 is the number of free parameters (M_pl, T_eff) in the explanatory model.

Our baseline models were published by Bodenheimer et al. (2003), hereafter BLL. As described in BLL, radius estimates, $R_{m_{i}}$, were computed with a Henyey-type planetary structure calculation, and the reader is referred to that paper for details regarding the input physics and assumptions. The model radii were tabulated at 4.5 Gyr for a grid of assumed M_pl and T_eff. We ignore the small dependence of radius on age for mature planets and use bilinear interpolation to obtain an estimate for $R_{m_{i}}$ at given M_pl and T_eff. Our estimates are drawn from BLL's model sequence of core-free solar-composition planets with no anomalous energy sources.

For the 90 transiting planets, we find χ²_m = 23.5 and χ²_null = 32.6. Not surprisingly, this result indicates that the baseline structural models can explain some, but by no means all, of the observed variation in planetary radii. As an alternative to bilinear interpolation between table values, it can be useful to have a simple fitting relation, R_pl(M_pl, T_eff). Defining m = log₁₀(M_pl/M_Jup) and t = T_eff/1000, we find that the two-dimensional polynomial fitting function,

$$\begin{eqnarray} R_{\rm pl}/R_{\rm Jup} &=& 1.08417+0.0940857\,m-0.242831\,m^{2} \nonumber \\ && +0.0947349\,m^{3} +0.0387851\,t+0.00243981\,mt \nonumber \\ && -0.0244656\,m^{2}t +0.0130659\,m^{3}t+0.0240409\,t^{2} \nonumber \\ && -0.0419296\,mt^{2}+0.00693348\,m^{2}t^{2}\nonumber \\ && +0.00302157\,m^{3}t^{2}, \end{eqnarray} \tag{ 3 }$$

provides an acceptable approximation to the BLL baseline structural models throughout the region where 0.1 M_Jup < M_pl < 10.0 M_Jup and 100 < T_eff < 2500.⁴

For each planet, we define a radius anomaly, ${\cal R}_{i}=R_{m_{i}}-R_{o_{i}}$, and look for correlations between $\cal R$ and other measured quantities (such as T_eff, T_⋆, M_⋆, etc.). Many authors (e.g., Enoch et al. 2011) have noticed that planetary radii tend to swell dramatically with increasing insolation. Figure 2 illustrates the significant correlation between $\cal R$ and planetary T_eff. Using a bootstrap replacement method (Press et al. 1992), we find a best-fit power-law dependence,

$$\begin{equation} {\cal R}\propto T_{\rm eff}^{1.4\,\pm \,0.6}. \end{equation} \tag{ 4 }$$

Zoom In Zoom Out Reset image size

Among the various mechanism that have been invoked to explain the radius anomalies, we expect that both Ohmic heating (Batygin & Stevenson 2010) and kinetic heating (Guillot & Showman 2002) should show a positive correlation between $\cal R$ and T_eff. We can ask, furthermore, whether the measured exponent, α = 1.4 ± 0.6, is consistent with either of these proposed mechanisms.

In the treatment of Batygin & Stevenson (2010), energy deposition in the planet is approximated by integrating over the resistivity in each mass element:

$$\begin{equation} {\dot{E}}= \int \!\! \int \!\! \int {{\bf J}^{2}\over {\sigma _r(r)}} dV, \end{equation} \tag{ 5 }$$

where J is the current density and σ_r(r) is the conductivity in the near-surface (but still adiabatic) layers where the resistive heating occurs. For order-of-magnitude purposes, they ignore the radial dependence, σ_r(r) ∼ σ_r. Further, it is assumed that the relevant currents are induced by the planet's intrinsic dipole field, B, so that

$$\begin{equation} {\bf J}=\sigma _r({\bf v}\times {\bf B}) \, \quad {\rm and} \quad \dot{E}\propto \sigma _r |{\bf v}|^{2}|{\bf B}|^{2}. \end{equation} \tag{ 6 }$$

We want to estimate the temperature dependence, $\dot{E}(T_{\rm eff})$, implied by the Batygin–Stevenson model. First, we need an estimate of how wind kinetic energy at adiabatic depth scales with temperature. One possibility is that there is a constant partitioning between thermal energy and kinetic energy, an assumption that conforms with the kinetic heating hypothesis of Guillot & Showman (2002), which posits that a fraction, η, of the total flux received by a planet is converted into kinetic energy by atmospheric pressure gradients, and that this energy is dissipated at depth by Kelvin–Helmholtz instabilities. In this case, one has v ∝ T²_eff. Alternately, one can examine three-dimensional hydrodynamical models for irradiated giant planets that experience a range of insolations. These models suggest a much weaker dependence of wind speed on T_eff at large optical depths where energy dissipation can effectively couple to the planetary structure. For example, the Lewis et al. (2010) model for Gliese 436b (M_pl = 24 M_⊕, T_eff = 707 K), which adopts a solar metallicity atmosphere, finds v ∼ 100 ms⁻¹ at P = 10 bar depth. This wind velocity is nearly identical to the average wind velocity at 10 bar depth found by Showman et al. (2009) for a solar metallicity atmospheric model of HD 189733 (M = 1.2 M_Jup, T_eff = 1202 K). For our rough analysis, we choose v ∝ T^1/2_eff, bearing in mind that the temperature's power-law index might reasonably fall anywhere in the range from 0 to 2.

Second, we follow Batygin & Stevenson (2010) and assume that the Elsasser Number, Λ = σ_iB²/2ρΩ, is roughly constant across the aggregate of transiting planets. (The quantity σ_i is the conductivity appropriate to the interior regions where the B field is assumed to be generated by dynamo action, ρ is the density of the planet, and Ω is its angular spin frequency.) To the precision of our analysis, this relation holds for the bounding cases of Jupiter and the Sun; Jupiter has a magnetic field roughly 10 times the solar value and has a spin period that is roughly 100 times shorter. Both bodies have similar densities and are highly conductive in their interiors. With Λ and σ_i assumed constant, we can remove the magnetic field dependence, yielding $\dot{E}\propto \sigma _r \rho \Omega T_{\rm eff}$. If we further assume synchronous rotation (which should hold for hot Jupiters on circular orbits), then 1/Ω ∝ a^3/2 and a ∝ T⁻²_eff (via Equation (1)), giving $\dot{E}\propto \sigma _r T_{\rm eff}^{4}$.

We need to identify how the dissipation-layer conductivity, σ_r, depends on temperature. The conductivity is directly tied to the ionization fraction, which in turn is determined by the Saha equation

$$\begin{equation} {n^{+}_j n_e \over n_j - n^{+}_j} = \left({m_e kT \over 2 \pi \hbar ^2} \right)^{3/2} {2 g^{+}_{j}\over {g_{j}}} \exp[ - I_j /kT], \end{equation} \tag{ 7 }$$

where the j chemical species typically include Na, K, Li, Rb, Fe, Cs, and Ca, and T is the local temperature. The quantity I_j is the ionization potential (for the first electron to be ionized) and ranges from 4.3 eV (for K) to 7.9 eV (for Fe), whereas g⁺_j and g_j are the degeneracies of states for the ions and neutrals, respectively, of species j. Note that n_e is the total electron number density (since all electrons can interact with each species of ion). In the relevant temperature regime, ionization fractions are low, so that n⁺_j ≪ n_j. Summing over all species and simplifying yields

$$\begin{equation} n_e = \left({m_e kT \over 2 \pi \hbar ^2} \right)^{3/4} \bigg[ \sum _j f_j n_j \exp[ - I_j / kT] \bigg]^{1/2}, \end{equation} \tag{ 8 }$$

where f_j = 2g⁺_j/g_j.

Next, we note that the electrical conductivity, σ_r, appropriate to the dissipative region has the form

$$\begin{equation} \sigma _r = {n_e \over n} {e^2 \over m_e A} \left({\pi m_e \over 8 k T} \right)^{1/2} \,, \end{equation} \tag{ 9 }$$

where A is a weighted cross section, see, e.g., Tipler & Llewellyn (2002). The quantity of interest is the index p defined by

$$\begin{equation} { p }= {T \over \sigma _r} {\partial \sigma _r \over \partial T} = - {1 \over 2} + {T \over n_e} {\partial n_e \over \partial T} \,. \end{equation} \tag{ 10 }$$

The index p allows us to estimate the temperature dependence of the conductivity, σ_r. Using our expression for n_e, we obtain

$$\begin{eqnarray} { p } &=& {1 \over 4} + {T \over 2} \bigg[ \sum _j f_j n_j \exp [ - I_j / kT] \bigg]^{-1} \nonumber \\ && \times {\partial \over \partial T} \bigg[ \sum _j f_j n_j \exp[ - I_j / kT] \bigg]. \end{eqnarray} \tag{ 11 }$$

After differentiation, this equation can be written in the (apparently) simpler form

$$\begin{equation} { p }= {1 \over 4} + { { \langle I \rangle }\over 2 kT}, \end{equation} \tag{ 12 }$$

where we have defined

$$\begin{eqnarray} { \langle I \rangle } & \equiv & \bigg[ \sum _j f_j n_j \exp[ - I_j / kT] \bigg]^{-1} \nonumber \\ && \times \bigg[ \sum _j I_j f_j n_j \exp[ - I_j / kT] \bigg]. \end{eqnarray} \tag{ 13 }$$

Since kT ≪ I_j, the exponential suppression factors dominate the behavior of this function. As a result, we can consider the largest term to dominate (which will be the term with the smallest ionization potential I_*):

$$\begin{equation} { p }\approx {1 \over 4} + {I_\ast \over 2 kT}. \end{equation} \tag{ 14 }$$

For T = T_eff = 2000 K, the quantity 2kT = 0.345 eV. If I_* = 4 eV, then p ≈ 12. It is clear that $\dot{E}(T_{\rm eff})$ depends sensitively on T_eff, and that for a given choice of I_*, the index p increases sharply with T_eff (for I_* = 4 eV, p = 15.7 at 1500 K, and p = 9.5 at 2500 K).

A similar result is obtained by Batygin & Stevenson (2010), who present a simple but accurate model for the conductivity as a function of radius, where σ ∝ exp [ − (r − R)/H], where H is the conductivity scale height (which is about twice the thermal scale height). Since H ∼ T, the index p ≈ (r − R)/H, which is of order 10.

Next, we need to establish the dependence of the radius anomaly, $\cal R$, on the internally generated power, $\dot{E}$. This depends on the planetary mass and the size of the assumed core—it is easier to inflate a low-mass planet without a core. For a typical hot Jupiter, such as HAT-P-13, models based on the BLL prescription indicate that by increasing $\dot{E}$ more than 10-fold, from 3.5 × 10²⁵ erg s^-1 to 5.6 × 10²⁶ erg s^-1, one inflates the radius from R = 1.2 R_Jup to R = 1.36 R_Jup (Batygin et al. 2009). This single example suggests that $\cal R$ depends on $\dot{E}$ raised to a small fractional power, x. Gu et al. (2004) have computed sequences of evolutionary models (again using a method similar to that described in BLL) which systematically explore how planetary radii depend on internally generated power. We obtained a fit to their sequences and find that an index x ∼ 1/6 applies over a wide range of masses for planets having log(R_pl/R_Jup) < 0.3. This regime covers nearly all of the known transiting cases. Using p = 12, we arrive at a predicted dependence

$$\begin{equation} {\cal R}\propto T^{\alpha } \, {\approx }\, T_{\rm eff}^{{\,(12+4)/{6}}} \approx T_{\rm eff}^{\,2.666}\,\,\,\,, \end{equation} \tag{ 15 }$$

for which the index α is more than one power of temperature steeper than suggested by the currently observed aggregate of transiting planets (${\cal R} \propto T^{1.4}$).

Caution regarding the quantitative aspect of the dependence given by Equation (15) is in order. Order-of-magnitude scaling arguments can at best give only rough insight into the physical balance that determines the radii of the hot Jupiters. As mentioned above, there are considerable uncertainties associated with the wind speeds at adiabatic depth, and the treatment leading to Equation (14) is highly idealized. We have assumed, furthermore, that the Ohmic heating rate is directly proportional to the conductivity, σ_r, in the weather layer. If the ionization fraction is large enough, a more fully self-consistent treatment must account for the back-reaction onto the wind velocity, v, generated by the Lorentz force

$$\begin{equation} \rho {d{\bf v}\over {dt}} = {1\over {c}} \sigma _r({\bf v}\times {\bf B})\times {\bf B}\,. \end{equation} \tag{ 16 }$$

Indeed, Perna et al. (2010a, 2010b) have suggested that the bulk Lorentz force acts as a drag term on the velocity field. In this event, higher temperatures increase ionization, generating Lorentz forces that brake the velocity field. Reduced wind speeds lower the rate of Ohmic dissipation, which in turn would decrease the temperature power-law index α. A complete understanding of this process requires a full three-dimensional MHD approach. While such simulations are not beyond current computational capabilities, the lack of detailed knowledge of the atmospheric conditions likely renders them somewhat premature at this time.

Transiting extrasolar planets tend to have radii that are larger than those predicted by structural models parameterized by M_pl and T_eff. The resulting radius anomalies, ${\cal R}$, are strongly correlated with effective temperature, with ${\cal R}\propto T_{\rm eff}^{1.4}$ providing a maximum reduction in variance. The simple scaling arguments of the previous section suggest that Ohmic heating may play a significant role and that further progress may be made as self-consistent MHD-mediated weather models are constructed.

A very rapid increase in $\dot{E}$ with temperature poses difficulties for some posited explanations for the radius anomalies. For example, the Guillot & Showman (2002) kinetic heating mechanism proposes that a constant fraction, η ∼ 0.01, of the total flux received by a planet is converted into kinetic energy by atmospheric gradients, and that this energy is dissipated at adiabatic depth. While the flow patterns on the surfaces of strongly irradiated extrasolar planets are still a matter of debate, an $\dot{E}\propto T_{\rm eff}^{4}$ scaling is implied by the kinetic heating hypothesis, leading to ${\cal R}\propto T^{2/3}$, which (to >1σ confidence) is weaker than the observed dependence. This argument does not imply that kinetic heating is absent, rather, the implication is that it is unlikely to be the primary mechanism that determines the radius anomalies.

Finally, we note that substantial variance in the observed radii remains after the ${\cal R}\propto T_{\rm eff}^{1.4}$ scaling has been removed. Inclusion of the best-fit βT^α_eff dependence into the model radii increases N_f from 2 to 4, while causing χ²_m as measured with Equation (2) to decrease from χ²_m = 23.5 to $\chi ^{2}_{m^{\prime }}=14.8$. Clearly, processes other than Ohmic heating are also contributing. For example, our base-line models unrealistically assume solar-composition planets. Both Jupiter and Saturn are substantially super-solar, and the small radii of planets such as HD 149026b and CoRoT 8b can only be understood if elements heavier than H/He contribute more than 50% of the planetary mass. Figure 3 indicates that there is evidence for a significant correlation between the residual radius anomalies and host star metallicities. This correlation is generally interpreted as evidence that metal-rich protoplanetary disks lead to planets with larger core masses. Furthermore, a fraction of the known transiting planets have eccentric orbits, meaning that interior heating from tidal orbital circularization must be ongoing in some of the observed planets. For modest eccentricities, and at a given time-average tidal quality factor, Q, tidal heating rates scale with R⁵_p (Murray & Dermott 1999). This strong dependence could be a determining factor in producing the residual radius anomalies for many of the observed planets. Furthermore, even if a planet's orbital eccentricity is currently zero, it may still be inflated as a consequence of thermal inertia from a tidal heating episode that occurred in the past (Ibgui & Burrows 2009).

Zoom In Zoom Out Reset image size

The number of transiting planets with well-determined properties will increase rapidly in the coming years. An expanded catalog of planets, in conjunction with improvements to the structural models, will allow empirical quantities, such as the power-law index "α" that we have used here, to be determined with increased confidence. In addition, follow-up campaigns from platforms such as JWST will improve our understanding of the atmospheric conditions on short-period planets, and we can look forward to a time when the construction of detailed, fully self-consistent three-dimensional MHD climate models are a sensible response to the quality and depth of the observational data sets.

The authors acknowledge useful discussions with Konstantin Batygin, and we are grateful to an anonymous referee for a prompt and insightful report. This work was funded by NASA grant NNX08AY38A and NASA/Spitzer/JPL grant 1368434 to G.L., and by NASA grant NNX07AP17G to F.C.A.