OHMIC DISSIPATION IN THE INTERIORS OF HOT JUPITERS

First we review the general picture put forward by Batygin & Stevenson (2010) (see also Perna et al. 2010a). Due to strong irradiation from the host star, the hot Jupiter has a thermally ionized atmosphere, coupling the magnetic field and the atmospheric flow. The magnetic field could be either the intrinsic planetary magnetic field generated by a dynamo in the deep interior or the external field from the host star. In either case, the evolution of the magnetic field in the wind zone is governed by the induction equation

$$\begin{equation} \frac{\partial {\textbf {\em B} }}{\partial {t}}=-\nabla \times \eta (\nabla \times {\textbf {\em B} })+\nabla \times (\textbf {\em v} \times {\textbf {\em B}}), \end{equation} \tag{ 1 }$$

where ${\bm v}$ is the wind velocity and η is the magnetic diffusivity of the atmosphere. Assuming a steady state and that the planet magnetic field in the outer layers is well represented by a curl-free dipole field $\textbf {\em B}_{\rm dip}$, the induced magnetic field $\textbf {\em b}$ is given by

$$\begin{equation} \nabla \times \eta (\nabla \times {\textbf {\em b}})=\nabla \times (\textbf {\em v}\times {\textbf {\em B}_{\rm dip}}). \end{equation} \tag{ 2 }$$

If the wind is predominately a zonal flow $\textbf {\em v}=v_\phi \hat{\phi }$, the induced field is toroidal and will penetrate into the interior of the planet, with associated poloidal currents that close in the interior given by $\textbf {\em J}=(c/4\pi)\nabla \times {\textbf {\em b}}$. The internal toroidal field is given by

$$\begin{equation} \nabla \times (\eta \nabla \times {\textbf {\em b}})=0 \end{equation} \tag{ 3 }$$

below the wind zone where velocities are negligible. For a given poloidal current J, the local ohmic dissipation rate is P_ohm = J²/σ, where σ is the electrical conductivity, related to the magnetic diffusivity by η = c²/4πσ.

Some intuition about the solution can be obtained by considering a plane-parallel model, which we illustrate in Figure 1. There we divide the planet into three layers, representing the outermost isothermal layer, with pressures lower than 60 mbars; the wind zone, between 60 mbars and 10 bars; and the interior of the planet. The vertical field B_z is sheared by the wind with velocity v_y(z)exp (ikx). Focusing on the innermost layer representing the planet interior, and assuming constant conductivity, Equation (3) gives an induced field $\textbf {\em b}\propto \exp (-{kz}+ikx)\hat{y}$ and associated vertical current 4πj_z/c = ikb. Both the field and current decrease exponentially in the interior on a length scale 2π/k. When the variation of conductivity with depth is included, we must solve

$$\begin{equation} \frac{d^2b}{dz^2}-k^2b+\frac{d\eta }{dz}\frac{db}{dz}=0, \end{equation} \tag{ 4 }$$

in which case the thickness of the penetration depth depends on the length scale on which the conductivity varies.

Zoom In Zoom Out Reset image size

This simple model illustrates that the interior solution depends on three factors: the geometry of the shearing velocity in the wind zone (which is described by k in this simple model), the magnitude of the induced field at the base of the wind zone, and the profile of the electrical conductivity in the interior. The approach that we pursue in this paper is to consider the first two factors as boundary conditions on the interior. We choose to parameterize our models by specifying the toroidal magnetic field at a pressure of 10 bars, which we will refer to as B_ϕ0.

To calculate the distribution of currents inside the planet in detail, we consider a simple geometry with a dipole field and a zonal flow $\textbf {\em v}=v_0\sin {\theta }\,\hat{\phi }$. To solve Equation (3) in spherical coordinates, we write the induced toroidal field:

$$\begin{equation} {\bm B_{\phi }}=\frac{g(r)}{r}\sin {\theta }\cos {\theta }\hat{\phi }, \end{equation} \tag{ 5 }$$

in which case the radial dependence part of the field is given by

$$\begin{equation} g^{\prime \prime }(r)-\frac{d\ln \sigma }{dr}g^{\prime }(r)-l(l+1)\frac{g(r)}{r^2}=0, \end{equation} \tag{ 6 }$$

where l = 2 is the index of associated Legendre polynomial P¹_l. Having found g(r) and therefore B_ϕ(r), the currents are determined by Ampère's law $\textbf {\em J}=(c/4\pi)\nabla \times \textbf {\em B}$. The ohmic power in the interior is

$$\begin{equation} P=\int \frac{J^2}{\sigma }\,dV\sim 4\pi \int \frac{\langle {J}\rangle ^2}{\sigma }r^2\,dr, \end{equation} \tag{ 7 }$$

where 〈J〉 = 〈J²_r + J_θ²〉^1/2 is the effective angle-averaged current at radius r. To get some feeling for the dependence of the field and current on depth, we can consider a power-law dependence σ∝r^α. In that case, the solution is B_ϕ∝r^β with

$$\begin{equation} \beta =\frac{(\alpha -1)+\sqrt{(1+\alpha)^2+24}}{2}, \end{equation} \tag{ 8 }$$

and 〈J〉 ∝ r^{β − 1}. For more complex wind geometries, which involve l > 2, the solution is ${\bm B_{\phi} }\,{\propto}\, r^lP_l^1({\rm cos}\theta)$ for constant conductivity. For example, this indicates that more zonal jets in the wind zone imply a shallower penetration depth for the induced field.

In fact, as we argue in the next section (Section 2.2), the internal heating is dominated by the lowest densities, since the local heating rate is inversely proportional to the electrical conductivity, which increases rapidly with increasing pressure. This means that the current J can be taken as a constant without making a significant error in the heating profile. We have confirmed this by comparing constant current solutions with detailed solutions of Equation (6).

For the constant current case, we compute the current as

$$\begin{equation} J={c\over 4\pi }{B_{\phi 0}\over R_J}, \end{equation} \tag{ 9 }$$

where R_J = 7 × 10⁹ cm is the radius of Jupiter and B_ϕ0 = B_ϕ(r, p = 10 bars). This means that there is a direct mapping between the chosen value of B_ϕ0, the radial current inside the planet J_r, and the local heating rate, taken to be J²_r/σ(r) per unit volume. Note that we do not take into account the averages over angle in Equation (7) or the true radius of the planet, and so our value of B_ϕ0 for a given amount of internal heating could be a factor of a few of the toroidal field in a model that self-consistently includes both the wind zone and the interior. Instead, our parameter B_ϕ0 should be interpreted as a measure of the internal heating (given by Equation (9) and then J²_r/σ locally).

It is useful to estimate the expected magnitude of ohmic heating and how it scales with parameters such as planet mass M. The first step is to use Equation (2) to estimate the expected strength of the induced field by dimensional analysis,

$$\begin{equation} \frac{b}{B_{\rm dip}}=R_M=\frac{4\pi \sigma {H}v}{c^2}, \end{equation} \tag{ 10 }$$

or

$$\begin{equation} b=B_{\rm dip}\left(\frac{\sigma }{10^7\,{\rm s}^{-1}}\right)\left(\frac{H}{0.01\,R_J}\right)\left(\frac{v}{1\,{\rm km}\,{\rm s}^{-1}}\right), \end{equation} \tag{ 11 }$$

where v is an average wind speed and σ is a typical value of electrical conductivity in the layer,³ and we take the vertical length scale to be the pressure scale height H.

In this paper, we take B_dip = 10 G as a standard value. Typical dipole field strengths for hot Jupiters have been estimated from scalings with planet parameters (see Trammell et al. 2011 for a detailed summary and discussion). Sánchez-Lavega (2004) argued that the field is generated by the dynamo action in the metallic region as in Jupiter (Stevenson 1983), with the field strength closely related to the rotation of the planet, B ∼ (ρΩη)^1/2. This predicts that the field on typical hot Jupiters should be a factor of a few smaller than that on the Jupiter, with a typical value of equatorial magnetic field B_eq  ∼  5 G. However, Christensen et al. (2009) argue that the field instead scales with the heat flux escaping from the conductive core at large enough rotation rate, giving B ∼ (ρF²_core)^1/3. This gives a field strength an order of magnitude larger than estimated with the previous method.

Given an induced field b, we can then estimate the expected ohmic power per unit mass, P_m = 〈J〉²/(ρσ). Approximating the angle-averaged current at the base of the wind zone using the constant current case as Equation (9) gives

$$\begin{eqnarray} P_m &=& 10^{-1}\ \mathrm{erg\ g^{-1}\ s^{-1}} \left(\frac{B_{\phi 0}}{10\, \mathrm{G}}\right)^2 \left(\frac{\sigma }{10^{6}\ \mathrm{s}^{-1}}\right)^{-1} \nonumber\\ &&\times\left(\frac{\rho }{10^{-4}\ \mathrm{g\ cm^{-3}}}\right)^{-1}. \end{eqnarray} \tag{ 12 }$$

Since σ increases dramatically in the core of the planet, heating at low density dominates. If the conductivity profile scales as exp (− r/H), for example, where we take the length scale as the pressure scale height H, the total power in the interior is

$$\begin{eqnarray} P_{\rm ohm, total}&\approx &4{\pi }P_{\rm top}{R}^2\rho {H} \end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray} &=&10^{23}\,\mathrm{erg\ s^{-1}}\,\left(\frac{B_{\phi 0}}{10\ \mathrm{G}}\right)^2\,\left(\frac{\sigma _t}{10^{6}\ \mathrm{s}^{-1}}\right)^{-1}\left(\frac{H}{0.01\ R_J}\right)\left(\frac{R}{R_J}\right)^2\nonumber \\ &=& 3\times 10^{22}\ \mathrm{erg\ s^{-1}}\,\left(\frac{B_{\phi 0}}{10\ \mathrm{G}}\right)^2\,\left(\frac{\sigma _t}{10^{6}\ \mathrm{s}^{-1}}\right)^{-1}\left(\frac{T}{1500\ {\rm K}}\right)\nonumber\\ &&\times \left(\frac{R}{R_J}\right)^4 \left(\frac{M}{M_J}\right)^{-1}, \end{eqnarray} \tag{ 14 }$$

where the subscript t indicates a quantity evaluated in the outermost regions of the interior just below the wind zone. Note that $H=\mathcal {R}T/g$, where we adopt $\mathcal {R}=3.64\times 10^{7}\ \mathrm{erg\ g^{-1}\ K^{-1}}$ for a hydrogen-molecule-dominated composition with helium fraction Y = 0.25. As we mentioned previously, because the total ohmic power is dominated by the heating at low pressure, it is not very sensitive to the radial profile of the current. The scaling in Equation (13) implies that P_ohm∝1/M when the radius and conductivity of the planet do not strongly depend on mass, a scaling that we find in our numerical solutions. More massive planets have less ohmic power for a given B_ϕ0 and temperature T_iso at the base of the wind zone.

We now make models of gas giants with given mass M and central entropy S and use them to calculate the ohmic dissipation in the planet, but without including the effect of ohmic heating on the planet structure. This allows us to calculate the ohmic heating within the convection zone and, by including this ohmic power in Equation (16), the corresponding slowing of the cooling rate.

We present the detailed microphysics in our planet model in the Appendix. Here, we only note the differences between our opacity and conductivity profiles and those used in other works. Our opacity profile uses a different extrapolation method in the intermediate pressure range between 10³ and 10⁵ bars from Paxton et al. (2011), who take the opacity for log R > 8 (where R = ρ/T³₆ is used in the opacity tables) to be a constant set by the value at log R = 8. As far as we are aware, opacity calculations for this intermediate pressure range have not been carried out. For planets undergoing a large amount of internal heating, the convection zone boundary moves into this region, and so knowing the opacity there is important for understanding the location of the convection zone boundary. As we describe later, this controls the contraction rate of the cooling planet. Our potassium conductivity profile (Equation (A2)) is the same as used by Perna et al. (2010a) but is different from Batygin & Stevenson (2010), who use an electron–neutron cross section of π(7.2 × 10⁻⁹ cm)² = 1.6 × 10⁻¹⁶ cm² rather than 10⁻¹⁵ cm² (Draine et al. 1983) and a slightly different thermal averaging factor for the velocity. The resulting difference is that the conductivity of Batygin & Stevenson (2010) is a factor of nine times larger than our conductivity.

The planet model is calculated by integrating outward from the center, following the convective adiabat until it intersects a radiative zone extending inward from the surface. When calculating the cooling curve of an irradiated gas giant, a reasonable approximation is to take the outer radiative zone of the planet to be isothermal. However, because ohmic heating is very sensitive to pressure, a small error in the determination of the convective boundary results in a much larger error in the ohmic power. To illustrate this, we have calculated models with an isothermal radiative layer and with a radiative layer that is in thermal equilibrium and carries a constant luminosity equal to the luminosity from the convection zone.

For the isothermal case, we integrate

$$\begin{eqnarray} \frac{dm}{dr}&=&4\pi \rho {r^2} \end{eqnarray} \tag{ 18 }$$

$$\begin{eqnarray} \frac{dp}{dr}&=&-\rho \frac{Gm}{r^2} \end{eqnarray} \tag{ 19 }$$

$$\begin{eqnarray} \frac{dT}{dr} &=&\frac{T}{P}\nabla \frac{dp}{dr} \end{eqnarray} \tag{ 20 }$$

outward from the center, taking ∇ = ∇_ad for T > T_iso (adiabatic interior) and ∇ = 0 for T < T_iso (isothermal layer). In the non-isothermal case, we take

$$\begin{eqnarray} \nabla &=&\min {(\nabla _{\rm rad},\nabla _{\rm ad})} \end{eqnarray} \tag{ 21 }$$

$$\begin{eqnarray} \nabla _{\rm rad}&= &\frac{3{\kappa }L}{16{\pi }c{\it GM}}\frac{p}{aT^4}. \end{eqnarray} \tag{ 22 }$$

We calculate a model for a given M and S by integrating outward from the center and inward from the surface to a matching pressure, p = 30 kbars. For the outward integration, we choose the central pressure p_c and cooling rate dS/dt (or equivalently cooling time t_S = S/|dS/dt|). For the inward integration, we start at a pressure of 10 bars and set the temperature there to be T_iso. We then integrate inward, choosing the luminosity L and radius R. A multi-dimensional Newton–Raphson method is used to find the correct choices of (p_c, t_S, R, L) that result in m, r, T, and L agreeing to within 1% at the matching pressure.

As an example, Figure 2 compares the entropy and temperature profiles for models with an isothermal and a non-isothermal radiative zone. In the isothermal case, we choose T_c = 3 × 10⁴ K, p_c = 2 × 10⁷ bars, and T_iso = 1500 K, which gives an M = 0.96 M_J, R = 1.25 R_J planet with core entropy S = 7.98 and convective zone boundary p_conv = 62.76 bars. The luminosity from the interior is L = 1.28 × 10²⁶ erg s⁻¹, giving t_S = 5.78 Gyr. With a radiative zone, we obtain the same mass, radius, and entropy with a convective zone boundary p_conv = 131.7 bars. The luminosity from the interior is L = 7.7 × 10²⁵ erg s⁻¹, and t_S = 10.5 Gyr to cool.

Zoom In Zoom Out Reset image size

In each case, the magnetic field structure in the planet interior is obtained by solving Equation (6) using the conductivity profile in the planet (shown in Figure 17), and then the ohmic heating profile is determined. For an induced magnetic field B_ϕ0 = 10 G at the bottom of the wind zone (p = 10 bars), we find P_ohm(p ⩾ p_conv) = 8 × 10²³ erg s⁻¹ for the isothermal model and P_ohm(p ⩾ p_conv) = 6.8 × 10²² erg s⁻¹ in the non-isothermal case. While the cooling time for the planet changes by a factor of two between the two models, the ohmic power changes by more than an order of magnitude. Therefore, it is crucial to locate the convective boundary accurately when calculating the ohmic power inside the convection zone.

In Figure 3, we compare the ohmic power calculated in this way, which includes the correct radial distribution of current, with the ohmic power calculated by assuming a constant radial current, independent of depth. The agreement is excellent (within a factor of two) except at the highest pressures in the central regions of the planet.

Zoom In Zoom Out Reset image size

The fact that the ohmic heating per unit mass rises rapidly to lower densities (Figure 3) suggests that the heating in the regions lying between the wind zone and the convection zone boundary will be larger than the heating in the convective interior. We include the ohmic heating in the radiative layer by allowing L to vary throughout the radiative zone with

$$\begin{equation} {dL\over dr}={J^2\over \sigma }. \end{equation} \tag{ 23 }$$

We do not include ohmic heating at pressures less than 10 bars. Instead, we specify the temperature T_iso at p = 10 bars and integrate inward. Of course, there could be significant ohmic heating within the wind zone at p < 10 bars, but we absorb this into the boundary condition. Note that this means that early in the lifetime of the planet, when the entropy is large enough that p_conv < 10 bars, the models here revert back to our previous models with no feedback. However, at those early times, ohmic heating is generally not yet important. Also note that in the models without feedback, the strength of the induced field B_ϕ does not influence the internal structure of the planet, whereas here a larger B_ϕ results in more heating in the radiative layer, which can push the convective boundary deeper.

In Table 1, we compare the models with feedback to our earlier models without feedback. The internal structures are shown in Figure 2. The planet radius does not vary much between different models. The biggest difference is in the position of convective zone boundaries (marked by black vertical bars in Figure 2), which results in a difference in the cooling luminosity and the ohmic heating in both the convective zone and atmosphere. Note that this means that the cooling history of a planet using these three approaches would be different, especially the time at which ohmic heating begins to become important for evolution. We use this feedback model for all the calculations carried on below.

The luminosity and ohmic power are shown as a function of entropy in Figure 4. As noted in particular by Arras & Bildsten (2006), Equation (15) shows that L∝M at fixed entropy, and we see that scaling in Figure 4. On the other hand, the ohmic power decreases with increasing M, as discussed in Section 2.2. We find that the decrease in P_ohm is well described by P_ohm∝R^2.4/M, which has a shallower dependence on R than in Equation (13) because of the dependence of the conductivity term on mass, which compensates the R⁴ term. In our feedback model, for a lower mass planet the higher atmospheric ohmic heating pushes the convective zone boundary slightly deeper, resulting in a higher conductivity at the top of the convection zone. Overall, lower mass planets generally have higher ohmic power deposited in the convection zone.

Zoom In Zoom Out Reset image size

Combined with the mass–luminosity dependence, the decrease of ohmic power with mass means that ohmic heating becomes important for lower mass planets at a much higher entropy than for more massive planets. The value of entropy at which ohmic heating becomes important depends on the boundary-induced field B_ϕ0. Turning this around, for an observed planet with measured radius and mass, we can infer the entropy and therefore derive a limit on the wind zone B_ϕ0 required for ohmic heating to be providing a significant part of the luminosity in that object. We carry out this procedure in Section 5, but we first describe our calculations of the time evolution of planets with ohmic heating.

Both Batygin et al. (2011) and Menou (2012) write down simplified models for the wind zone dynamics including the effects of magnetic drag. In both cases, following Perna et al. (2010a), the magnetic drag force is assumed to be $\textbf {\em J}\times \textbf {\em B}/c$ per unit volume, with the current $\textbf {\em J}$ set by a balance between the shearing of the magnetic field by the fluid and ohmic diffusion of magnetic field lines against the fluid motion, $\textbf {\em J}=\sigma \textbf {\em v}\times \textbf {\em B}/c$ (Section 2). However, the dynamical balance in the two models is quite different. Menou (2012) writes the force balance for the equatorial flow as (see also Showman et al. 2010)

$$\begin{equation} 0=-\frac{v_{\phi }^2}{R_P}+\frac{\mathcal {R}{\Delta }T_{\rm horiz}\Delta \ln {p}}{R_P}-\frac{v_{\phi }B_r^2}{4\pi \rho \eta }. \end{equation} \tag{ 24 }$$

The first two terms represent a balance between the advective term and the horizontal driving from the day–night temperature difference ΔT_horiz. This balance is thermal-wind-like in that the horizontal pressure gradients require a vertical gradient in the fluid velocity v_ϕ over a vertical pressure scale Δln p. The final term represents the magnetic drag force, again integrated over a vertical scale Δln p. Batygin et al. (2011), on the other hand, consider the meridional circulation induced by magnetic drag on the azimuthal flow, so that, for example, the latitudinal force balance is fv_y = v_ϕ/τ_L, where f = 2Ωsin θ is the Coriolis parameter and τ_L is the magnetic drag timescale. Their solution represents a thermal wind balance involving the equator-pole temperature gradient, modified by magnetic drag.

In both cases, magnetic drag limits the fluid velocity at high temperatures, where the large degree of ionization and therefore large electrical conductivity result in strong coupling of the fluid and magnetic field. Balancing the second and third terms in Equation (24) gives v_ϕ∝η when magnetic drag dominates, and therefore the magnetic Reynolds number R_M = v_ϕH/η becomes almost constant, varying only slowly with temperature. Similarly, Equation (16) of Batygin et al. (2011) has two possible limits, either v_ϕ∝η when the lateral temperature gradient is large, in which case R_M becomes almost constant at large T_iso, or v_ϕ∝η² when the drag timescale is comparable to the rotation period, while the lateral gradient of temperature is still small, in which case R_M∝η declines rapidly at large T_iso.

A dynamical model including both day–night driving and meridional circulation with magnetic drag is not yet available. For our purposes, we have implemented the model of Menou (2012) as described by Equation (24), with the day–night temperature difference given by

$$\begin{equation} {\Delta }T_{\rm horiz}=\left\lbrace \begin{array}{@{} ll}\frac{T_{\rm day}}{2} \left(\frac{\tau _{\rm adv}}{\tau _{\rm rad}}\right) & \tau _{\rm adv}<\tau _{\rm rad}\\ \ms \frac{T_{\rm day}}{2} & \tau _{\rm adv}>\tau _{\rm rad} \end{array}.\right. \end{equation} \tag{ 25 }$$

In Equation (25), T_day is the dayside-averaged temperature considering a dilution factor of 0.5, T⁴_irr = 2T_day⁴ = 4T⁴_eq, and the advective and radiative timescales are

$$\begin{equation} \tau _{\rm adv}=\frac{R_P}{v_{\phi }} \end{equation} \tag{ 26 }$$

$$\begin{equation} \tau _{\rm rad}=\frac{C_pp}{g{\sigma _{SB}}T_{\rm day}^3}. \end{equation} \tag{ 27 }$$

Note that these timescales are evaluated at the outermost pressure, which, following Menou (2012), is taken to be 60 mbars, the estimated location of the thermal photosphere. This is the reason for adopting the thermal timescale appropriate for an optically thin region, so that τ_rad∝p in Equation (27); in deeper, optically thick layers, τ_rad has an extra factor of the optical depth τ, leading to τ_rad∝p² (e.g., Figure 3 of Showman et al. 2008). The magnetic drag term is also evaluated at p = 60 mbars; this term is integrated over height, but since σ decreases with increasing pressure (for an isothermal layer), the dominant contribution to the integral is from the lower limit on pressure, and so η and ρ are evaluated there. This is an important difference from Batygin et al. (2011), who evaluated their magnetic drag timescale at p = 10 bars, which gives a drag time an order of magnitude longer than we find here. Based on that estimate, Batygin et al. (2011) concluded that the drag timescale was always much longer than a rotation period.

We solve Equation (24) for v_ϕ as a function of T_eq and find the corresponding value of the induced field B_ϕ from Equation (10) for different values of the dipole field B_dip. For Δln p = 0.9, we reproduce the results of Menou (2012) (see his Figure 1), but we also consider larger values of Δln p. Menou (2012) models the weather layer with a modest vertical extension around 1 pressure scale height. We also solve the equation with Δln p = 3 for typical values in hot Jupiter atmosphere as reported by Showman et al. (2010), and Δln p = 5 for a wind zone extending to p ∼ 10 bars, for comparison with Batygin & Stevenson (2010) and Batygin et al. (2011). The effect of varying Δln p on B_ϕ is shown in the left panel of Figure 10. For numerical convenience, we fit the B_ϕ–T_eq relation with the following:

$$\begin{equation} \frac{1}{B_{\phi }(T_{\rm eq})}=\frac{1}{B_{\rm adv}}+\frac{1}{B_{\rm drag}}, \end{equation} \tag{ 28 }$$

where

$$\begin{eqnarray} B_{\rm adv} &=& 2.8\times 10^6\ {\rm G}\ T_{\rm eq} \exp \left(-{2.53\times 10^4\over T_{\rm eq}}\right) \nonumber\\ &&\times \left({\Delta \ln {p}\over 3}\right)^{1/2} \left({B_r\over 10\ {\rm G}}\right) \end{eqnarray} \tag{ 29 }$$

and

$$\begin{equation} B_{\rm drag}=1125\ {\rm G}\ \left({T_{\rm eq}\over 1000\ {\rm K}}\right)\left({\Delta \ln p\over 3}\right)\left({B_r\over 10\ {\rm G}}\right)^{-1}. \end{equation} \tag{ 30 }$$

This reproduces B_ϕ to within ≈10% for T_eq in the range 1100–2200 K.

Zoom In Zoom Out Reset image size

The transition to the regime where R_M is approximately constant occurs when τ_drag exceeds τ_adv, where the magnetic drag timescale is τ_drag = 4πρη/B²_r = η/v²_A, where v_A is the Alfvén speed, and again the timescale is evaluated at the top of the wind zone (p = 60 mbars here). These timescales are plotted as a function of T_eq in Figure 10. Changing the wind zone thickness from Δln p = 0.9 to Δln p = 5 moves the transition temperature from T_eq ≈ 1400 K to 1700 K.

To use this value of B_ϕ as a boundary condition for our evolutionary models, we must relate the temperature T_eq at low pressure to the temperature T_iso at p = 10 bars. This relation depends on the details of energy transport in the wind zone, including the effects of ohmic heating, and needs to be studied further. Here, we adopt the atmospheric temperature profile from Guillot (2010) and keep in mind the uncertainty in the relation between T_eq and T_iso when interpreting our results below. The relation from Guillot (2010) is (see his Equation (29))

$$\begin{equation} T^4=\frac{3T_{\rm irr}^4}{4}f\left[\frac{2}{3}+\frac{1}{\gamma \sqrt{3}}\right], \end{equation} \tag{ 31 }$$

where γ is the ratio between visible and infrared opacities and f = 1/2 for a dayside average or f = 1/4 for an average over the whole surface. Choosing γ = 0.4, as appropriate for a planet like HD 409658b (see, e.g., Figure 1 of Hubeny et al. 2003) and a dayside average f = 0.5, we obtain T_iso = 0.94T_irr = 1.33T_eq. In the following section, we will use this relation to infer the appropriate value of T_iso from the T_eq of observed planets. We note here that we do not have a good knowledge of what the γ parameter would be for most of the observed planets. While γ parameter could vary in a very large numerical range, the ratio between T_iso and T_eq only changes within a factor of a few [$0.99(\gamma \rightarrow \inf)<(T_{\rm {iso}}/T_{\rm {eq}})<3.05(\gamma =0.01)$]. Since the observed properties of planets give T_eq and thus the boundary-induced field, changing the value T_iso/T_eq is equivalent to shift inside the plane T_iso − B_ϕ given by Figure 8. Generally, a smaller T_iso/T_eq is favored to inflate the planet with the same B_ϕ.

In Figures 11–13, we compare our results with the observed properties of transiting planets taken from the TEPcat transiting planet catalog,⁴ which gives the planet mass, radius, and equilibrium temperature T_eq = T_{⋆, eff} (R_⋆/2a)^1/2, where T_{⋆, eff} is the stellar effective temperature. As the ages of most stars are unknown or highly uncertain, we assume an age of 3 Gyr when comparing with the observed planets.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

First, Figure 11 shows the effect of increasing B_ϕ0 at fixed T_iso on the planet radius. To help compare with the data, we divide the observed sample into two groups with either T_eq > 1500 K (green points) or <1500 K (red points) and show theoretical curves for either T_eq = 1316 K or 1692 K (these two temperatures correspond to T_iso = 1750 and 2250 K, respectively). We see that for the low-T_eq group, an induced field of 10–100 G can explain most of the observed radii, while the high-T_eq planets need at least B_ϕ0 = 1000 G to match the observed radii. It is clear that a higher induced magnetic field is needed to explain a given radius at higher equilibrium temperature.

Next, we use the wind zone model described in Section 5.1 to calculate B_ϕ0 as a function of T_eq, assuming canonical values B_r = 10 G and Δln p = 3. In the top panel of Figure 12, we show the radius as a function of T_eq, with the colors representing three different bins in planet mass. There exists a clear correlation between the radius and T_eq, in both the observations and the models. In addition, we see that the amount of inflation is also strongly dependent on the planet mass. Planets within the mass bin 0.3–0.6 M_J agree quite well with our ohmic heating model. However, ohmic heating clearly cannot explain planets with mass ∼1 M_J and large inflated radii ≳ 1.4 R_J. Ohmic heating can help to increase the radius (for comparison the dashed line shows models with no ohmic heating), but not enough to match the observed value. This is a consequence of the increased power needed to maintain the radius of a massive planet at a particular value, as well as the reduced ohmic heating power at larger masses.

In the lower panel of Figure 12, we show the radius as a function of mass. As in Figure 11, we divide the data into two temperature ranges and show the model results for two representative temperatures, now using the wind zone model to specify the value of B_ϕ0 for each temperature. The parameter region where ohmic heating has the largest effect is high-temperature, low-mass planets. The radii of the low-temperature group (T_eq < 1500 K) can almost all be explained without ohmic heating. For the high-temperature group, ohmic heating can explain the observed radii of low-mass planets, but most of the radii of the high-temperature group lie well above the models, especially at large planet masses ≳ 1 M_J.

In Figure 13, we show the ratio between observed and predicted radii R_obs/R_pred against T_eq (upper panel) and against M (lower panel) for each observed planet. In this case, we use the observed values of M and T_eq to calculate the evolution of the planet, and, in the absence of stellar ages, we take R_pred to be the radius at 3 Gyr. In the upper panel, we see that for T_eq ≳ 1600 K, there are many planets whose radii lie above the predicted values. The slow increase of B_ϕ0 with T_iso at large T_iso due to the magnetic drag term results in a much weaker dependence of R_pred on T_eq than observed, and most outliers lie at the highest temperatures. In the lower panel, we see that the majority of the unexplained objects (R_obs > R_pred) are at larger masses M ≳ 0.7 M_J. We note that the choice of estimating the planet radius at 3 Gyr is not critical for the above picture. Constraining ourself within the time range of 1–5 Gyr, the predicted radius only varies within several percent.

To look in more detail at the effect of our assumed wind zone model on how successfully we are able to reproduce the observed radii, in Figure 14 we show the results in the B_ϕ0–T_eq parameter space. For each observed planet, we first make a model with no ohmic heating, varying the internal entropy S at the measured M and T_eq until we match the measured radius R_p. Then we calculate the value of B_ϕ0 required in that model for the ohmic power in the convection zone P_ohm to be 30% of the planet's luminosity. We color code the data points according to whether they are successfully explained by our time evolutions, i.e., whether they have R_pred larger or smaller than R_obs in Figure 13. These two groups of data points lie on either side of the B_ϕ0–T_eq relation from the wind zone model (solid curve). This shows that the approach of using a structural model with no ohmic heating (we refer to this as a "no feedback" model in Section 3) to estimate the critical magnetic field is a good approximation of our detailed time-evolution models including feedback.

Zoom In Zoom Out Reset image size

Comparing the red points in Figure 14 with the solid curve gives a sense of how far short the ohmic heating model falls in explaining the most inflated planets. For example, HAT-P-32 is about a factor of 3–4 above the curve, so that the heating rate (∝B²) needs to be increased by about an order of magnitude to explain the observed radius. It is interesting that most of the unexplained objects lie within a factor of three in terms of B_ϕ0 of the wind zone model. Figure 14 helps to show what changes to the wind zone model would explain more of the observed objects. We have assumed the relation T_iso = 1.33T_eq (from Equation (31)); a larger factor between T_iso and T_eq would move the solid curve to the left, allowing ohmic heating to explain the radii of low-T_eq planets such as WASP-06. The dashed and dotted curves show the effect of changing B_r. Increasing B_r from 10 to 100 G does increase B_ϕ0 at low temperatures, but it reduces B_ϕ0 at high temperatures where magnetic drag is enhanced. A larger depth Δln P would help to reduce the number of discrepant objects since both B_adv and B_drag increase with Δln P (Equations (29) and (30)).

We discuss the microphysics input in our gas giant models here. We adopt the equation of state from Saumon et al. (1995) with helium fraction Y = 0.25. In order to maximize the planet radius, we do not include a solid core or elements heavier than helium.

The radiative opacity is taken from Freedman et al. (2008), and in the core we include thermal conduction by electrons from Potekhin & Chabrier (2010). The transition from radiative to conducive opacity occurs at a pressure that is greater than the maximum pressure of 300 bars covered by the Freedman et al. (2008) tables. In the intermediate regime, we assume the scaling κ∝p^0.5. The opacity profile for our standard model is shown in Figure 15, overplotted with opacity taken from MESA using the same planet structure (Paxton et al. 2011).

Zoom In Zoom Out Reset image size

The electrical conductivity has contributions from alkali metal ionization in the outer layers and hydrogen in the interior. In the upper atmosphere of hot Jupiters, the conductivity is set by the ionization of alkali metals. For potassium, which has the lowest ionization potential,⁵ the Saha equation (Balbus & Hawley 2000; Perna et al. 2010a) gives the ionization fraction x_k = n_e/n as

$$\begin{eqnarray} x_k&=&\left[\frac{f_k}{n}\left(\frac{m_ek_BT}{2\pi \hbar ^2}\right)^{3/2}e^{-4.35 {\rm eV}/k_BT}\right]^{\frac{1}{2}} \\ \nonumber &=&1.03\times 10^{-3}\ T_3^{5/4}e^{-25.19/T_3} \left(\frac{f_K}{10^{-7}}\right)^{1/2} \left(\frac{p}{1\, {\rm bar}}\right)^{-1/2}, \end{eqnarray} \tag{ A1 }$$

where f_K is the number fraction of potassium. Although potassium dominates, we also include the contribution of Na, Mg, and Fe in the ionization balance to sum up the total ionization fraction. The ionization fraction of each alkali metal is computed separately by assuming a balance independent of the presence of other elements. We do not include the contribution of Al in the calculation, which is likely condensed out (Lodders 1999). But our results are not very sensitive to elements beyond potassium. This is illustrated in Figure 16, which shows the contributions to the ionization level from K, Na, and Al at different pressures. Once the ionization fraction is determined, the conductivity is σ = n_ee²/m_eν, where the collision frequency of electron–neutral collisions is ν_en = n_n〈σv〉_e given by Draine et al. (1983) as

$$\begin{equation} \langle \sigma {v}\rangle _e=10^{-15}\left(\frac{128k_BT}{9\pi {m_e}}\right)^{1/2}\ {\rm cm}^3\ {\rm s}^{-1}. \end{equation} \tag{ A2 }$$

The conductivity is then

$$\begin{equation} \sigma =8.8\times 10^{-2}\ {\rm S\ m^{-1}}\ \left({x\over 10^{-7}}\right)\left({T\over 1500\ {\rm K}}\right)^{-1/2}. \end{equation} \tag{ A3 }$$

Zoom In Zoom Out Reset image size

In the deeper part of the planet, the hydrogen is ionized by high pressure and the conductivity is dominated by electron–proton collisions. In the fully degenerate limit, ν_epd = 4e⁴m_eΛ/3πℏ² = 1.8 × 10¹⁶ s⁻¹. In the non-degenerate limit, ν_epnd = 6.4 × 10²³ s⁻¹ρY_eT^−3/2, in which Y_e is the electron fraction. We interpolate between the two limits to give an estimation of the total contribution: ν⁻²_ep = ν_epd⁻² + ν⁻²_epd. We also include the conductivity at intermediate densities as given by Liu et al. (2006). Before the hydrogen molecule is fully ionized, the band-gap of hydrogen will diminish with increasing pressure to a level where there is a significant contribution to the conductivity. Liu et al. (2006) give this as

$$\begin{equation} \sigma _s= \sigma _0 \exp \left(\frac{-E_g(\rho)}{k_BT}\right), \end{equation} \tag{ A4 }$$

where between 0.2 Mbars and 1.8 Mbars, E_g = 20.3–64.7ρ, where E_g is in eV and ρ is in mol cm⁻³, and σ₀ = 3.4 × 10²⁰ s⁻¹ exp (− 44 ρ). The overall conductivity is constructed as σ = σ_s + n_e e²/m_eν = σ_s + 1.52 × 10³²ρY_e/ν. The collisional frequency ν is the sum of electron–neutral and electron–proton collisions. A typical conductivity profile and the contribution of different components are shown in Figure 17.

Zoom In Zoom Out Reset image size