Structure and Evolution of Internally Heated Hot Jupiters

TL;DR

This study models the internal heating of hot Jupiters, showing that deeper and more intense heating can explain their observed large radii, with implications for understanding their thermal evolution.

Contribution

It systematically explores the effects of varying depth and intensity of internal heating on hot Jupiter radii without assuming specific heating mechanisms.

Findings

Deeper heating at 100 bars explains observed radii with only 1% of stellar irradiation.

Shallow heating at 1-10 bars is ineffective unless very intense.

Heating at 10^4 bars suppresses cooling as effectively as core heating.

Abstract

Hot Jupiters receive strong stellar irradiation, producing equilibrium temperatures of $1000 - 2500 \ \mathrm{Kelvin}$. Incoming irradiation directly heats just their thin outer layer, down to pressures of $\sim 0.1 \ \mathrm{bars}$. In standard irradiated evolution models of hot Jupiters, predicted transit radii are too small. Previous studies have shown that deeper heating -- at a small fraction of the heating rate from irradiation -- can explain observed radii. Here we present a suite of evolution models for HD 209458b where we systematically vary both the depth and intensity of internal heating, without specifying the uncertain heating mechanism(s). Our models start with a hot, high entropy planet whose radius decreases as the convective interior cools. The applied heating suppresses this cooling. We find that very shallow heating -- at pressures of $1 - 10 \ \mathrm{bars}$ -- does not significantly suppress cooling, unless the total heating rate is $\gtrsim 10\%$ of the incident stellar power. Deeper heating, at $100 \ \mathrm{bars}$, requires heating at only $1\%$ of the stellar irradiation to explain the observed transit radius of $1.4 R_{\rm Jup}$ after 5 Gyr of cooling. In general, more intense and deeper heating results in larger hot Jupiter radii. Surprisingly, we find that heat deposited at $10^4 \ \mathrm{bars}$ -- which is exterior to $\approx 99\%$ of the planet's mass -- suppresses planetary cooling as effectively as heating at the center. In summary, we find that relatively shallow heating is required to explain the radii of most hot Jupiters, provided that this heat is applied early and persists throughout their evolution.

Figures

Equations

$$\frac{dm}{dr}=4\pi r^2\rho ,$$

$$\frac{dP}{dm}=-\frac{Gm}{4\pi r^4},$$

$$\frac{dL}{dm}=\frac{d\left(L_{\mathrm{rad}}+L_{\mathrm{conv}}\right)}{dm}={ϵ}_{\mathrm{grav}}+{ϵ}_{\mathrm{extra}},$$

$$\frac{dT}{dm}=-\frac{GmT}{4\pi r^{4}P}\nabla .$$

$${\nabla }_{\mathrm{rad}}=\frac{3}{64\pi \sigma G}\frac{\kappa LP}{m{T}^4}.\text{vspace}0.01cm$$

$$\Gamma ={\int }_0^M{ϵ}_{\mathrm{extra}}dm,$$

$$\gamma \equiv \Gamma /L_{\mathrm{irr}}=10^{-5}-0.1$$

$$\frac{d{S}_c}{dt}\sim \frac{(\Gamma -L)}{T_{c}M}.$$

$$\Gamma {\tau }_{\mathrm{dep}}\sim \Gamma \frac{P_{\mathrm{dep}}\kappa }{g}\gtrsim L_{\mathrm{irr}}.$$

$$L_{\mathrm{cool}}={\int }_0^M-T\frac{dS}{dt}dm,$$

$$R=R_0+\Delta R_0\left[{\left(1+\eta \frac{\Gamma {\tau }_{\mathrm{dep}}}{L_{\mathrm{irr}}}\right)}^{\delta }-1\right],$$

Full text

††slugcomment: Accepted at ApJ

1. Introduction

Since the first transit detections of an extrasolar planet (Charbonneau et al., 2000; Henry et al., 2000), the close-in extrasolar giant planet, or “hot Jupiter,” population has proved to be enigmatic. It was recognized early that many of these hot Jupiters have radii larger than expected from standard models considering only cooling from a high-entropy initial state (for reviews see Fortney et al., 2010; Baraffe et al., 2010, 2014; Laughlin & Lissauer, 2015). A mass-radius diagram of the hot Jupiter sample is shown in Figure 1, with cooling models from Fortney et al. (2007) assuming varying amounts of incident stellar flux over-plotted. Approximately one-half of the observed hot Jupiters have radii above expectations from evolutionary models. Though including irradiation in these models causes a deep radiative zone which pushes the radiative-convective boundary (RCB) to higher pressures and thereby reduces the planetary cooling rate, this only affects radii by $\lesssim 20\%$ (Guillot et al., 1996; Arras & Bildsten, 2006; Fortney et al., 2007).

There is also a trend of increasing planet radius with increasing equilibrium temperature, found by Laughlin et al. (2011) and visible by eye in Figure 1. Further analysis shows that planets with equilibrium temperatures $T_{\mathrm{eq}}\lesssim 1000\ \mathrm{Kelvin}$ have radii which match the expectations of purely cooling models (Demory & Seager, 2011; Miller & Fortney, 2011). Re-inflated hot Jupiters have been proposed as potential additional evidence that stellar irradiation is what drives radius inflation (Lopez & Fortney, 2016), and there is a growing observed population of re-inflated hot Jupiters around post-main sequence stars (Grunblatt et al., 2016; Hartman et al., 2016). Hence, the mechanisms responsible for enlarged hot Jupiters are somehow correlated with incident stellar flux.

There are three classes of explanations for the radius anomaly of hot Jupiters: tidal mechanisms, modifications to our understanding of the microphysics of hot Jupiters, and incident stellar-flux driven mechanisms (Weiss et al., 2013; Baraffe et al., 2014). Tidal dissipation was the first proposed mechanism for explaining the bloated radii of many hot Jupiters (Bodenheimer et al., 2001), and as such has been followed up by a variety of studies (Jackson et al., 2008; Ibgui & Burrows, 2009; Miller et al., 2009; Ibgui et al., 2010; Leconte et al., 2010). However, tidal dissipation nominally requires the eccentricity of the planet to be pumped up by an external companion, as tidal dissipation damps eccentricity. Arras & Socrates (2010) proposed that thermal tides in the atmosphere of the planet itself can torque it away from synchronous rotation, enabling the interior of the planet to couple to the stellar gravitational tidal force and cause dissipation, increasing the viability of this mechanism. The second class of mechanisms, those which do not require internal heating, includes enhanced opacities (Burrows et al., 2007) and lowered internal heat transport due to double-diffusive convection (Chabrier & Baraffe, 2007; Leconte & Chabrier, 2012). Though certainly important for understanding the internal structure of gas giants, it is unclear if these microphysical mechanisms to increase the radius of hot Jupiters scale with the incoming stellar flux.

The third class of mechanism involves transport of a fraction of the heat in the high-entropy irradiated region of the planet downward to the interior, where it then dissipates and modifies the entropy of the internal adiabat. These mechanisms are all necessarily linked to the vigorous atmospheric circulation in hot Jupiter atmospheres, where $\sim\mathrm{km/s}$ east-west winds are driven by the large day-to-night temperature contrasts in these atmospheres (Showman & Guillot, 2002; Cooper & Showman, 2005; Menou & Rauscher, 2009; Showman et al., 2009, 2010; Rauscher & Menou, 2010; Showman & Polvani, 2011; Perna et al., 2012; Mayne et al., 2014; Kataria et al., 2016; Komacek & Showman, 2016; Komacek et al., 2017). This class of mechanisms can additionally be split into two classes: the hydrodynamic and magnetohydrodynamic mechanisms. The latter, known as “Ohmic dissipation,” utilizes currents driven in the partially ionized atmosphere (threaded by a dipolar planetary magnetic field) which then resistively dissipate in the interior of the planet. This mechanism was proposed by Batygin & Stevenson (2010), with a variety of follow-up studies (Perna et al., 2010; Batygin et al., 2011; Heng, 2012; Huang & Cumming, 2012; Menou, 2012a, b; Rauscher & Menou, 2013; Wu & Lithwick, 2013; Rogers & Komacek, 2014; Rogers & Showman, 2014; Ginzburg & Sari, 2016) performed either supporting or refuting to varying degrees the Ohmic dissipation hypothesis. Models of purely kinematic magnetohydrodynamics normally attain dissipation rates large enough to explain radius anomalies, but self-consistent magnetohydrodynamic simulations do not.

Even for nominally simpler hydrodynamic dissipation mechanisms, several details remain poorly understood. The first of this class of mechanism proposed was downward kinetic energy transport (Guillot & Showman, 2002; Showman & Guillot, 2002). In this picture, the $\sim 10-100\ \mathrm{m/s}$ vertical winds in hot Jupiter atmospheres transport energy downward. This energy is then dissipated in shear layers near the RCB through, for example, Kelvin-Helmholtz instabilities. Another possible mechanism, the “Mechanical Greenhouse” (Youdin & Mitchell, 2010), involves downward transport of heat by forced turbulent mixing in the outer radiative zone of the planet. Building upon this, Tremblin et al. (2017) showed using a two-dimensional steady-state dynamical model that the large-scale circulation itself can produce enough downward entropy transport to explain the radius of HD 209458b. These are viable mechanisms to explain the hot Jupiter radius anomaly. However, no deep time-dependent simulations of hot Jupiter atmospheres have yet been performed to quantify the three-dimensional atmospheric circulation at levels where dissipation would strongly affect the interior entropy.

Though there have been a number of mechanisms proposed to explain the hot Jupiter radius anomaly, none of those discussed above have been shown to explain the entire sample of hot Jupiters. As a result, we focus in this work not on specific mechanisms but on how internal heating of any strength at any location affects the evolution and structure of hot Jupiters. Spiegel & Burrows (2013) showed that adding heat in the atmosphere of a hot Jupiter was more efficient at increasing pressure levels (see their Figure 4), with a markedly nonlinear relationship between transit radius and pressure of dissipation. They also showed that heating in the center of the planet enables a long-timescale radius equilibrium to be reached, which is not possible when heating is applied in the atmosphere alone. However, they did not examine the impact of heating at a given pressure level on the structural evolution of the hot Jupiter, as their numerical model decoupled the evolution of the radiative atmosphere and convective interior of the planet. As a result, the relationship between radius and deposited power from Spiegel & Burrows (2013) is determined only for jointly specified surface gravity and effective temperature, even though both depend on the deposited power itself.

It is possible, as has been shown by the self-similar analytic model of Ginzburg & Sari (2015), that outer convection zones can be forced due to strong heating. Their analytic model shows that the presence or absence of an outer convective zone is determined by a comparison between the combination of the strength and depth of deposited heating and the incoming stellar flux. The analytic model of Ginzburg & Sari (2015) also predicts the transition, in terms of heating strength and depth, where heating in the outer radiative zone begins to have a substantial effect on the resulting transit radius. In this work we test these predictions using a detailed numerical evolutionary model, which also enables us to consider the regime where heating occurs deeper than the inner radiative-convective boundary. This regime is not directly considered in Ginzburg & Sari (2015). Another possibility which we take into account is that dissipation occurs at a given structural location within the planet, for example the radiative-convective boundary, which is not fixed in pressure with time.

Here we use the stellar & planetary structure code MESA (Paxton et al., 2011, 2013, 2015) to compute the evolution of a hot Jupiter with different amounts of heating deposited at different depths or structural locations. These models allow a detailed understanding of the structural and evolutionary pathways of planets subject to different types of heating. This paper is organized as follows. Section 2 describes the numerical setup used, including choices for heating profiles. Section 3 describes the general theoretical framework that we use to interpret our results. Section 4 displays our results for how heating of varying strength and deposition depth affects the structure and evolution of hot Jupiters. We compare these results to the analytic theory of Ginzburg & Sari (2015) in Section 5, and delineate conclusions in Section 6.

2. Methods

2.1. Planetary Structure Equations

We use MESA to solve the following equations of “stellar” structure (Chandrasekhar, 1939; Kippenhahn et al., 2012), applied to giant planets:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Equation (1) is the mass conservation equation giving the enclosed mass, $m$, at radius, $r$, with mass density $\rho$. Equation (2) expresses hydrostatic equilibrium of the pressure, $P$, with the gravitational constant, $G$. Equation (3) states energy conservation, where the outgoing luminosity, $L$, computed by MESA includes both radiative ($L_{\mathrm{rad}}$) and (where convectively unstable) convective ($L_{\mathrm{conv}}$) components. The relevant energy sources include the cooling term, $\epsilon_{\mathrm{grav}}=-TdS/dt$, i.e. the loss of entropy, $S$, that drives gravitational Kelvin-Helmholz contraction, and $\epsilon_{\mathrm{extra}}$, the “extra” energy deposition described in detail below. Equation (4) is the energy transport equation, in terms of the logarithmic gradient, $\nabla\equiv d\mathrm{ln}T/d\mathrm{ln}P$, where $T$ is temperature. MESA uses mixing length theory to calculate $\nabla=\nabla_{\mathrm{rad}}$, the radiative gradient, (or $\nabla\simeq\nabla_{\mathrm{ad}}$, the convective gradient) in radiative (or convective) regions (respectively). Specifically, the Schwarzchild criterion sets the temperature gradient $\nabla=d\mathrm{ln}T/d\mathrm{ln}P$ to the smaller of the adiabatic gradient $\nabla_{\mathrm{ad}}$ or the radiative gradient

[TABLE]

In Equation (5), $\sigma$ is the Stefan-Boltzmann constant and $\kappa$ is the local opacity. The set of Equations (1) - (4) is closed with a thermodynamic equation of state (Saumon et al., 1995; Paxton et al., 2013), and with tabulated opacities, needed for $\nabla_{\rm rad}$, that assume a dust-free Solar composition (Freedman et al., 2008; Paxton et al., 2011, 2013).

2.2. Numerical Model

2.2.1 Setup

The initial models for our MESA evolutionary calculations are computed as described in Paxton et al. (2013). We choose a mass relevant for HD 209458b and assume an initial radius of $2.3\ R_{\mathrm{Jup}}$ to compute an adiabatic starting planet model. The long-term evolution, at times $\gtrsim 10\ \mathrm{Myr}$, is independent of the initial radius and hence is also independent of the initial entropy (Arras & Bildsten, 2006). However, note that the Kelvin-Helmholtz contraction timescale of a planet is inversely related to both its radius and luminosity. This means that if we instead chose a larger starting radius (and hence entropy) the planet would have a more rapid early evolution but after this initial cooling phase would end up on the same evolutionary track as a planet with a smaller initial radius. For comparison, the present-day transit radius of HD 209458b is $\approx 1.35\ R_{\mathrm{Jup}}$. We evolve the initial model with external irradiation and internal heating, as described below.

Equations (1)-(4) are solved using the Henyey method (Henyey et al., 1959; Bodenheimer et al., 2007; Kippenhahn et al., 2012), including automatic mesh refinement (Paxton et al., 2011). The outer boundary pressure and temperature are fixed at the location where the optical depth $\tau$ to outgoing radiation is $2/3$. We ensure that the Henyey residuals are small enough to not affect structure for all of our evolution calculations. We describe the specified external irradiation in Section 2.2.2, our choices for $\epsilon_{\mathrm{extra}}$ in Section 2.2.3, and our parameter choices in Section 2.2.4.

2.2.2 Irradiation

The irradiation of the exoplanet is included as an energy generation rate $\epsilon_{\mathrm{extra}}=F_{\star}/(4\Sigma_{p})$, which is applied where the outer mass column $\Sigma$ is less than the chosen $\Sigma_{p}$ (see Section 2.2.4 for specific parameter values). With $F_{\star}$ as the incoming stellar flux, this heating drives a heat flux of $F_{\star}/4$, equivalent to the average of the incoming radiation over the planet’s spherical surface. This is the $F_{\star}-\Sigma_{p}$ irradiation routine (Paxton et al., 2013), which has also been used by Owen & Wu (2016); Valsecchi et al. (2015). This irradiation method agrees with the detailed “grey irradiated” solutions of Guillot (2010), as shown in Paxton et al. (2013) (see their Figure 3, which shows agreement within $\sim 1-2\%$ after 100 Myr of evolution).

2.2.3 Internal Heating

In this work, we consider the impacts of heating at different depths and structural locations within a given hot Jupiter on the resulting evolution of the planet. Though this heating can be interpreted as due to dynamical processes which deposit heat from near-photospheric levels to greater depths, we do not consider the impacts of specific heating mechanisms. Instead, we add an extra dissipation $\epsilon_{\mathrm{extra}}$ (computed each time step) that is taken to be a Gaussian with standard deviation $0.5H$, where $H=p/(\rho g)$ is a pressure scale height. This heating profile is similar to that applied by Spiegel & Burrows (2013), but as we are using a global planetary structure code (instead of a detailed atmosphere model matched onto an adiabat) we can calculate the impact of heating near or below the RCB.

We consider heating deposited at pressures of $P_{\mathrm{dep}}=1-10^{4}\ \mathrm{bars}$ (distributed as described above), at the center of the planet, and at the inner RCB. The integrated heating rates,

[TABLE]

are set to different fractions of the irradiation as

[TABLE]

for this grid of simulations. Heating rates are written as a fraction of the incident stellar power $L_{\mathrm{irr}}=2.4\times 10^{29}\ \mathrm{erg}\ \mathrm{sec}^{-1}$ relevant for HD 209458b.

We choose an upper limit of $\gamma=0.1$ because a $100\%$ conversion of starlight to kinetic energy, which is transported to depth via either hydrodynamic or magnetohydrodynamic mechanisms, is unrealistic. Thus we choose $10\%$ conversion as an upper limit, consistent with previous studies. In the case of hydrodynamic dissipation, the upper end of these values for the normalized heating rate $\gamma$ can be motivated from the fraction of atmospheric kinetic energy dissipated to heat, which is unknown for hot Jupiters but is expected to be $\sim 1-10\%$ based on Earth (Peixoto & Oort, 1992; Guillot & Showman, 2002; Schubert & Mitchell, 2013). It is expected that for Ohmic dissipation $\gamma\lesssim 1-5\%$ (Perna et al., 2010; Batygin et al., 2011; Rauscher & Menou, 2013; Ginzburg & Sari, 2016) and may be up to two orders of magnitude smaller (Rogers & Showman, 2014; Rogers & Komacek, 2014), and for tidal dissipation $\gamma\lesssim 10\%$ (Arras & Socrates, 2010) and may be as small as $\sim 0.1\%$ (Bodenheimer et al., 2001; Jackson et al., 2008). When we apply heating at the RCB, $P_{\mathrm{dep}}=P_{\mathrm{RCB}}$ follows the location of the boundary. Heating at this boundary can be motivated by either downward mixing of heat by large-scale eddies or communication between atmospheric motions and the deep interior resulting in shear instabilities at this interface (Guillot & Showman, 2002; Showman & Guillot, 2002; Youdin & Mitchell, 2010).

2.2.4 Parameter Choices

We keep all external planetary parameters constant, running a suite of simulations varying only the integrated heating rates $\Gamma$ and deposition pressure $P_{\mathrm{dep}}$. We use a mass ($0.69\ M_{\mathrm{Jup}}$), composition ($Y=0.24,Z=0.02$), and irradiation flux relevant for HD 209458b, taken from Guillot & Showman (2002). As in Guillot & Showman (2002), our model does not include a heavy element core. As a result, our resulting radii are upper limits for the assumed heavy element composition.

For our irradiation routine, we choose $\Sigma_{p}=250\ \mathrm{g}\ \mathrm{cm}^{-2}$, corresponding to an opacity $\kappa_{\mathrm{vis}}=4\times 10^{-3}\ \mathrm{cm}^{2}\ \mathrm{g}^{-1}$ to incoming radiation, as in Fortney et al. (2008); Guillot (2010); Owen & Wu (2016). We use $F_{\star}=1.0012\times 10^{9}\ \mathrm{erg}\ \mathrm{cm}^{-2}\ \mathrm{s}^{-1}=1200F_{\oplus}$ (where $F_{\oplus}$ is Earth’s incident flux), which is equivalent to an equilibrium temperature $T_{\mathrm{eq}}=1450\ \mathrm{Kelvin}$ assuming full longitudinal redistribution of heat. We evolve our modeled planets to 5 Gyr, a typical age for main sequence systems and a nominal stopping point to compare with observed transit radii (e.g. Huang & Cumming, 2012; Wu & Lithwick, 2013).

3. Cooling Regimes of Hot Jupiters

3.1. Definition of Regimes

We here describe the main cooling regimes that occur in the evolution of giant gaseous planets, subject to different levels of irradiation and different amounts and depths of internal heating. Understanding these regimes and the evolution between them facilitates interpretation of our results in Section 4. Figure 2 shows a schematic (which will be referred to throughout this section) of the various possible structural regimes a gas giant may lie in (top panels), along with example evolutionary paths between these regimes (bottom panels).

First, consider regime 0, which applies for planets which have a cooling luminosity that exceeds any external irradiation (i.e. the cooling luminosity $L_{\mathrm{cool}}\gg$ the irradiation power $L_{\mathrm{irr}}$). Such a planet has a fully convective interior and radiates from a RCB at the planet’s photosphere. Regime 0 applies to weakly irradiated planets or to irradiated planets with high initial entropies. Our strongly irradiated models are never in regime 0, and if we were to start with higher entropies this phase would be very brief.

As an irradiated regime 0 planet cools and $L_{\mathrm{cool}}$ declines, it eventually enters regime 1, where external irradiation exceeds the planet’s cooling luminosity (i.e. $L_{\mathrm{irr}}\gg L_{\mathrm{cool}}$). For regime 1, we also require that the cooling exceeds any extra heat that is input from e.g. turbulent, Ohmic or tidal dissipation. Thus in regime 1, this extra heat does not yet significantly affect evolution. As the planet cools in regime 1, the outer radiative zone recedes below the photosphere, i.e. the RCB increases in pressure with time as the entropy of the convective interior decreases, because the majority of planetary cooling is from the convective interior (Arras & Bildsten, 2006). Planets without any extra heat will remain in regime 1. Between the convective interior and radiative exterior, intermediate convective zones and radiative windows may develop due to abrupt changes in the EOS and/or opacity. We ignore this detail for now, by noting that the key issue is the cooling rate from the innermost RCB, which connects to the central entropy. We will shortly address intermediate convective zones triggered by extra heating, which is not relevant in regime 1.

Regime 2, where the planet’s extra heating exceeds the cooling rate (i.e. the deposited power $\Gamma\gg L_{\mathrm{cool}}$), is the most relevant for our study. The transition from regime 1 to 2 occurs naturally due to the gradual decline in $L_{\rm cool}$ with time as the internal entropy decreases. During regime 2, the extra heating reduces the cooling luminosity as either luminosity “replacement” or luminosity “suppression.” To try to be exhaustive, we describe four sub-regimes, which we label:

• •

2(a): The limiting case of heating at the very center of the planet, or at the boundary between a solid core and the gaseous envelope.

• •

2(b): Heating in the convective interior at an intermediate radius.

• •

2(c): Heating outside the convective interior which triggers an outer convective zone above – and a corresponding radiative window below – the heating level.

• •

2(d): Heating outside the convective interior in a fully radiative zone.

Regime 2(a) represents the simplest case of heating at the deepest possible level of the interior. Here the added heat simply supplies, or “replaces,” some of the planet’s luminosity as $L_{\rm cool}=L_{\rm RCB}-\Gamma$. Only hot Jupiters in regime 2(a) can reach an exact steady state with $L_{\rm cool}=0$ within stellar main sequence lifetimes. However, planets in other regimes can (and do) have a cooling time that is longer than the age of the system, with negligible entropy loss.

Regime 2(b) is similar to 2(a) except heating is no longer exactly at the center. The heating is still within the convective interior and can still be thought of as replacing some of the RCB luminosity to decrease cooling. In this regime, the convective region below the heating layer emits cooling luminosity at a finite (but perhaps negligible) rate. A lower bound on the initial cooling luminosity is set by $L_{\mathrm{rad,ad}}(P_{\mathrm{dep}})$ which means the radiative luminosity along the adiabat at the location of applied heating. 111This estimate assumes a single depth of heat deposition, and could be refined for more broadly distributed heating. However, this refinement is not needed for a basic understanding. This minimum initial amount of cooling can be quite low for large $P_{\mathrm{dep}}$.

In regimes 2(c) and 2(d), the heating is outside the convective interior. In both regimes, heating pushes the RCB of the convective interior deeper in the planet, which suppresses cooling. The difference between the regimes is that in 2(c) the heating is deep and intense enough to trigger convection above – and a radiative window below – where heat is deposited. In regime 2(c), the outer boundary of the convective interior (below the radiative window) is pushed to greater depths, which significantly reduces cooling compared to regime 2(d). Ginzburg & Sari (2015) analytically derived the conditions for triggering an outer convective zone, explaining the resulting strong reduction in planetary cooling.

3.1.1 Comparison with the Evolutionary Stages of Ginzburg & Sari (2016)

Before discussing the evolution between regimes, it is useful to compare our framework to the analysis of Ginzburg & Sari (2016), who considered the stages that a planet with applied internal heating evolves through (see their Appendix A and Figure 7). Note that Ginzburg & Sari (2016) consider a different heating profile: a power-law in optical depth with a cutoff above a specified optical depth value. However, a comparison is still possible as the steepness of their power-law nicely corresponds to different cases of deep vs. shallow heating that we study here.

Stage 1 of Ginzburg & Sari (2016) corresponds to our regime 0 of a fully convective planet that cools independent of heating or irradiation. For hot Jupiters, we again emphasize that this phase is either non-existent or very brief, only occurring in the initial evolution. Stage 2 of Ginzburg & Sari (2016) is the same as our regime 1. We define this regime in terms of heating and cooling rates, and Ginzburg & Sari (2016) express the limits equivalently in terms of optical depths in their self-similar model. Stage 3 of Ginzburg & Sari (2016) corresponds to our regimes 2(c) and 2(d), which in turn corresponds to the case I and II (respectively) that Ginzburg & Sari (2016) describe for the depth of heating. Our regime 2(a) corresponds to the special case of central heating (in case I of Ginzburg & Sari, 2016) in either stage 3 (while cooling is proceeding) or stage 4 (once steady state is reached) of Ginzburg & Sari (2016).

Along with all these similarities, there are two subtle but notable differences. We do not describe the general case of their stage 4, a steady state in which the deep interior is fully radiative and isothermal, all the way to the center. The overlap with the limiting case of stage 4 mentioned above is for heating that occurs precisely at the center. We do not describe the general case of stage 4 in Ginzburg & Sari (2016) simply because main sequence stellar lifetimes are not long enough for hot Jupiters to reach this stage. Specifically we find that even with deep heating, a radiative window may never open (see e.g. our $\gamma=10^{-2}$ heating efficiency at $P_{\rm dep}=10^{4}$ bars case in Section 4) or never extends close to the center (e.g. our $\gamma=10^{-3}$ at $P_{\rm dep}=10^{4}$ bars case). This is because by the time the radiative window opens (if it does), the deep heating has reduced the cooling luminosity to very low values. As a result, the Kelvin-Helmholtz timescale is much longer than main sequence lifetimes, and this timescale steadily increases as the radiative window deepens.

The second difference is that Ginzburg & Sari (2016) do not explicitly describe our regime 2(b), during which the internal heating is both significant and within the convective interior. Our regime 2(b) is distinct from the closest related stages of Ginzburg & Sari (2016); internal heating is not relevant in their stage 2 and heating is disconnected from the convective interior in their stage 3. We find that regime 2(b) is significant for understanding hot Jupiter evolution. For one thing, some planets with deep heating will spend their entire late stage cooling (i.e. after regime 1, which is independent of internal heating) in regime 2(b), never opening a radiative window. Second, for planets that do evolve from regime 2(b) to 2(c) by opening a radiative window, the suppression of cooling is much more complete during regime 2(b). In these cases, the final radius correlates positively with the (logarithmic) fraction of time spent in regime 2(b). Our numerical results will more clearly show the importance of regime 2(b) in interesting regions of parameter space. To our knowledge, a simple analytic theory that includes the equivalent of our regime 2(b) has not yet been published.

3.2. Evolution

To better understand these regimes, it helps to consider how a planet enters regime 2, which is when the heating first becomes significant. A crucial issue is the depth that characterizes heat deposition, $P_{\mathrm{dep}}$, relative to the RCB depth, $P_{\mathrm{RCB}}$, at this moment. With “deep heating” $P_{\mathrm{dep}}>P_{\mathrm{RCB}}$ and regime 2 begins as 2(b) [or 2(a) if heating is precisely at the center]. With “shallow heating” $P_{\mathrm{dep}}<P_{\mathrm{RCB}}$ and regime 2 begins as 2(d).

The next issue is how regimes evolve as the planet cools. In considering this issue, we will explain how regime 2(c) arises. The bottom panels of Figure 2 illustrate possible evolutionary pathways in regime 2, assuming that the heat input remains steady at a fixed depth. In regime 2(d), planetary cooling pushes $P_{\mathrm{RCB}}$ to greater depths. The planet remains in regime 2(d), as depicted in the bottom left panel. Note that continued cooling, and the associated global radius decrease, will not cause a planet to leave regime 2(d). To see this, note that the temperature-pressure profile is given by $d\mathrm{ln}(T)/d\mathrm{ln}(P)=\mathrm{min}(\nabla_{\mathrm{ad}},\nabla_{\mathrm{rad}})$, and that both the adiabatic and radiative gradients are explicitly independent of radius (see Equation 5). In principle a fully radiative end state to regime 2(d) is possible, but in practice the timescale to reach this state is extremely long.

Evolution in regime 2(a) is depicted in the bottom central panel of Figure 2. The planet never leaves regime 2(a) if the central heating remains steady. Cooling and entropy loss can proceed until the planet reaches steady state with $L_{\mathrm{cool}}=0$, as described above. In practice, however, the cooling time simply becomes greater than the age of the system. Following Arras & Bildsten (2006), we can use the entropy equation to estimate the time-derivative of the characteristic internal entropy, assuming heating in the internal convective region:

[TABLE]

The internal entropy $S_{c}$ always decreases with time, as $L=\Gamma+L_{\mathrm{cool}}>\Gamma$. However, if $\Gamma\gg L_{\mathrm{cool}}$, the entropy decrease will be effectively zero. This leads to a steady state where $L_{\mathrm{cool}}\rightarrow 0$ and all of the outgoing luminosity from the RCB is supplied by the deposited heat.

Cooling evolution in regime 2(b) is the most complex, as illustrated in the bottom right panel of Figure 2. Cooling in regime 2(b) leads to an increase in $P_{\mathrm{RCB}}$ and corresponding decrease in $L_{\mathrm{cool}}$. However, $P_{\mathrm{RCB}}$ cannot be deeper than $P_{\mathrm{dep}}$ in regime 2(b), by definition, and $L_{\mathrm{cool}}$ cannot drop below the limit described above, i.e. $L_{\mathrm{rad,ad}}(P_{\mathrm{dep}})$. Furthermore, regime 2(b) should not simply evolve into regime 2(d) because with $\Gamma\gg L_{\mathrm{cool}}$, the heating is sufficient to trigger convection as the cooling drops. Instead, continued cooling in regime 2(b) leads to the opening of a radiative window below $P_{\mathrm{dep}}$, i.e. a transition to regime 2(c) as depicted in Figure 2. Once in regime 2(c), cooling proceeds qualitatively as in 2(d), except with an outer convective layer. Hence, the inner RCB retreats with time and the planet remains in regime 2(c).

In summary, in regime 2 with steady heating a planet that enters regime 2(a), (c), or (d) should remain in that specific regime, while regime 2(b) can transition to regime 2(c). This basic understanding is reflected in our numerical evolutionary models. Clearly more complicated behavior is possible if the heating is very broadly distributed or evolving in time. These more complicated cases are not our focus, but they can probably be understood with a combination of the regimes described here. The notable exception is re-inflation, which we do not attempt to describe here because the processes by which a planetary interior gain entropy remain uncertain. Next, in Section 4 we describe our numerical results, using the discussion here as a backbone for understanding the various structural regimes present and the evolution between them.

4. Results

4.1. Radii of Internally Heated Hot Jupiters

Figure 3 shows the transit radius222To calculate the transit radius from the photospheric radius, we use the isothermal limit of Guillot (2010) (their Equation 60), setting the ratio of visible to infrared opacities $\kappa_{v}/\kappa_{th}=0.4$. after 5 Gyr of evolution for our grid of simulations varying $P_{\mathrm{dep}}$ and $\Gamma$. Generally, we find that the final radius increases monotonically with increasing $P_{\mathrm{dep}}$ and $\Gamma$, as expected from previous numerical (Spiegel & Burrows, 2013) and analytic (Ginzburg & Sari, 2015) work. However, there two key interesting features imprinted on top of this trend.

First, there is a large increase in transit radius between the cases with heating at $P_{\mathrm{dep}}=10$ and $100\ \mathrm{bars}$. Secondly, the model transit radii are essentially the same for all of the cases with deeper heating, namely with heating either at $P_{\mathrm{dep}}=10^{4}\ \mathrm{bars}$, heating at the center, or heating at the RCB. Figure 4 shows that heating at $100\ \mathrm{bars}$ still encloses $>99.99\%$ of the mass of the planet. Even heating at $10^{4}\ \mathrm{bars}$, which is in the inner convective region, encloses $\approx 99\%$ of the mass. Hence, heating at very shallow regions can greatly affect evolution, in some cases having a similar effect to heating at the center of the planet. We will focus on these two key features in most of the discussion on internal structure and evolution that follows.

4.2. Structure & Evolution

4.2.1 Structure: Varying Depth of Heating

The trends in transit radius with varying $P_{\mathrm{dep}}$ and $\gamma$ shown in Figure 3 can be understood by examining the internal structure. Our goal is to understand the two specific features introduced above. First, all cases with $P_{\mathrm{dep}}\gtrsim 10^{4}\ \mathrm{bars}$ show almost the same radius inflation for a given heating rate $\gamma$. Secondly, there is a large ($\sim 25\%$) jump in radius between $P_{\mathrm{dep}}=10-100\ \mathrm{bars}$. As we show below, the larger radius with deeper heating corresponds to the opening of an internal radiative window, i.e. cooling in regime 2(c) instead of regime 2(d), as introduced in Section 3.

Figure 5 shows pressure-temperature profiles at the end-state of evolution for varying $P_{\mathrm{dep}}$ with a fixed normalized heating rate $\gamma=10^{-2}$. These profiles show a clear bifurcation between shallow heating ($P\leq 10\ \mathrm{bars}$) and deep heating ($P\geq 100\ \mathrm{bars}$). Simulations with shallow heating at $P_{\mathrm{dep}}=1-10\ \mathrm{bars}$ have a deep outer radiative zone, extending to $P_{\mathrm{RCB}}>10^{3}\ \mathrm{bars}$. By contrast, with deep heating at $P_{\mathrm{dep}}\geq 100\ \mathrm{bars}$ the outer radiative zone only extends to $P\sim 10\ \mathrm{bars}$. The shallow cases have cooled in regime 2(d) while the deep cases end in regime 2(c), or 2(a) with heating directly at the center. With deeper heating, the transition to a steeper adiabatic profile at lower pressures clearly allows the planet to reach a higher central temperature and thus central entropy and transit radius. Ultimately though, it is the suppression of cooling by the deeper heating that allows this higher internal entropy at late times.

Note that the outer convective zones with $100\ \mathrm{bars}\leq P_{\mathrm{dep}}\leq 10^{4}\ \mathrm{bars}$ do not match onto the same internal adiabat. Deep heating forces outer convective zones with corresponding internal radiative windows, as found in both numerical (Guillot & Showman, 2002; Batygin et al., 2011; Wu & Lithwick, 2013) and analytic (Ginzburg & Sari, 2015) studies. Radiative windows require the internal adiabat to have lower entropy than the outer convective zone (Guillot et al., 1994). Nevertheless, the central entropy still remains much larger than for the shallow heating cases which fail to trigger outer convective zones.

Figure 6 shows profiles of $\nabla,\nabla_{\mathrm{ad}},\nabla_{\mathrm{rad}}$ for simulations with the normalized heating rate $\gamma=10^{-2}$. The two cases of $P_{\mathrm{dep}}=10$ and $100\ \mathrm{bars}$ illustrate the transition from shallow to deep heating, respectively. The profiles in Figure 6 show that heating at $P_{\mathrm{dep}}=100\ \mathrm{bars}$ triggers an outer convective region, with an inner radiative window below. Since $\nabla_{\mathrm{rad}}\propto LP$, the same heating ($L$) at a deeper pressure forces a secondary convective region. For heating at $P_{\mathrm{dep}}=10\ \mathrm{bars}$, the increase in $\nabla_{\mathrm{rad}}$ is too small to trigger convective instability.

The amount of heat required to force a detached outer convective zone can be quantified as (Ginzburg & Sari, 2015):

[TABLE]

The quantity $\Gamma\tau_{\mathrm{dep}}\sim 10^{2}L_{\mathrm{irr}}$ for $P_{\mathrm{dep}}=100\ \mathrm{bars}$ and $\Gamma\tau_{\mathrm{dep}}\sim L_{\mathrm{irr}}$ for $P_{\mathrm{dep}}=10\ \mathrm{bars}$, which is hence just shallow enough to not force a secondary convective zone. As a result, the simulations with $P_{\mathrm{dep}}\leq 10\ \mathrm{bars}$ match onto a lower internal adiabat, naturally explaining the $\sim 25\%$ smaller transit radius.

4.2.2 Evolution: Structural Quantities

Thus far we have only presented the end-state of evolution, without examining any temporal changes in structural quantities. However, examining the cooling history of our models is crucial for understanding trends with varying $P_{\mathrm{dep}}$ and $\Gamma$. Here we investigate the evolution of two slices of our grid, varying $P_{\mathrm{dep}}$ with fixed normalized heating rate $\gamma=10^{-2},10^{-3}$. Figure 7 displays the time-evolution of photospheric radius, central entropy, the gravitational cooling luminosity

[TABLE]

and the pressure at the inner RCB for $\gamma=10^{-2}$ and $10^{-3}$.

The early evolution of these structural quantities is similar for both $\gamma=10^{-2}$ and $10^{-3}$. All models show an identical cooling phase over the first $\sim 10\ \mathrm{Myr}$, which extends to $\sim 50\ \mathrm{Myr}$ with weaker heating of $\gamma=10^{-3}$. This uniform cooling phase is called regime 1 in Section 3. This phase ends when the gravitational cooling luminosity falls below the applied heating rate, explaining why this cooling phase lasts longer with weaker heating. After this point, the evolution diverges based on whether the heating is shallow (i.e. at $1-10\ \mathrm{bars}$) or deep ($\geq 100\ \mathrm{bars}$).

When heating is shallow, the planet continuously cools over time in regime 2(d) (i.e. with heating in the radiative exterior). The central entropy and thus planetary radius decrease smoothly with time. The RCB gradually retreats deeper into the planet, which causes the smooth decline in cooling rate.

For cases with deeper heating, $P_{\mathrm{dep}}\geq 100\ \mathrm{bars}$, the planet next reaches a quasi-equilibrium state, during which the planetary radius and entropy remain (nearly) constant over some period of time. The bottom panel of Figure 7 shows that during this equilibrium the interior RCB lies outside (at lower pressure than) the heating zone. For the case of $P_{\mathrm{dep}}=100\ \mathrm{bars}$, $P_{\mathrm{RCB}}\sim 15$ and $70\ \mathrm{bars}$ for $\gamma=10^{-2}$ and $10^{-3}$ in this state, respectively. Thus the hot Jupiters are mostly in regime 2(b) during this quasi-equilibrium phase.

The central heating case, i.e. regime 2(a), stays in the near equilibrium state. This central heating case rapidly approaches a true equilibrium with $L_{\mathrm{cool}}\rightarrow 0$ in which heating supplies the total interior luminosity. This behavior is evident in the sharp drop-off in $L_{\mathrm{cool}}$, see the third row of Figure 7.

For the cases with deep atmospheric heating at $10^{2}\ \mathrm{bars}\leq P_{\mathrm{dep}}\leq 10^{4}\ \mathrm{bars}$, the quasi-equilibrium state eventually comes to an end and the planet begins cooling and contracting again. The bottom row of Figure 7 shows that the departure from quasi-equilibrium roughly corresponds to a sharp increase in $P_{\mathrm{RCB}}$, i.e. the opening of a radiative window which marks the transition from regime 2(b) to 2(c). When the radiative window opens, $L_{\mathrm{cool}}$ does not increase, but it stops declining rapidly with time. This more slowly evolving $L_{\mathrm{cool}}$ is sufficient to eventually cause a noticeable decrease in planetary entropy and radius.

The case with $P_{\mathrm{dep}}=10^{4}\ \mathrm{bars}$ is notable because the evolution in radius and central entropy is nearly indistinguishable from the case with central heating. The behavior of $L_{\mathrm{cool}}$ is noticeably different for the two heating locations, with $L_{\mathrm{cool}}$ dropping to arbitrarily low levels at late times with applied central heating. However, with $P_{\mathrm{dep}}=10^{4}\ \mathrm{bars}$, $L_{\mathrm{cool}}$ stays just above $10^{23}\mathrm{erg}\ \mathrm{s}^{-1}$, for both $\gamma=10^{-2}$ and $10^{-3}$. Crucially, this finite level of cooling corresponds to Kelvin-Helmholtz cooling times of $\sim 10^{3}\ \mathrm{Gyr}$, much longer than the system age. Thus, departures from a true equilibrium with no cooling are hard to discern.

A subtle point about the $P_{\mathrm{dep}}=10^{4}\ \mathrm{bars}$ case is that the floor in $L_{\mathrm{cool}}$ at late times has a different origin in the $\gamma=10^{-3}$ and $10^{-2}$ cases, despite their similar value of $L_{\mathrm{cool}}$. For $\gamma=10^{-3}$, the plateau in $L_{\mathrm{cool}}$ occurs when a radiative window opens. This window opens at $\sim 500\ \mathrm{Myr}$, i.e. $\sim 10\%$ of the final age. However, the window is so deep, approaching a pressure of $10^{5}\ \mathrm{bars}$, that cooling is very slow. This example emphasizes that the opening of the radiative window [i.e. entering regime 2(c)] does not necessarily correspond to measurable radius contraction over stellar lifetimes. For $\gamma=10^{-2}$, a radiative window does not open. In this case the floor in $L_{\mathrm{cool}}$ corresponds to the minimum value of $L_{\mathrm{rad,ad}}(P_{\mathrm{dep}})$ discussed in Section 3.

4.2.3 Evolution: Luminosity and Temperature Profiles

To understand in detail how the pressure at which heat is deposited changes the internal cooling rate, we examine runs at fixed central entropy. Though these model planets have the same central entropy and radius, we are sampling different times in their evolution and hence different internal cooling rates. Figure 8 shows luminosity-pressure profiles for a normalized heating rate $\gamma=10^{-2}$ and varying deposition pressure at fixed central entropy $S_{c}/(k_{b}N_{A})=9$, which is approximately the equilibrium entropy for central heating at this value of $\gamma$. Luminosity is an integrated quantity from the center, and as a result it monotonically increases outward. There are three main features in this luminosity profile, from the center outwards: an initial rise to the level of $L_{\mathrm{cool}}$, a secondary rise to the heating level $\Gamma$, and a third to the irradiation level $L_{\mathrm{irr}}$. The profiles in the outer two levels are almost the same for all $P_{\mathrm{dep}}$, but the internal cooling rates (shown by the dashed lines) vary by four orders of magnitude across the sample. This is because if heating is deep enough to affect the internal structure, radiative losses effectively occur from just below the heating layer. The cooling rate then decreases with increasing $P_{\mathrm{dep}}$ (at fixed central entropy), due to the longer cooling timescales at greater depths.

To examine further the differences between the cases with $P_{\mathrm{dep}}=10$ and $100\ \mathrm{bars}$, we show the time-evolution of their luminosity and temperature profiles with $\gamma=10^{-2}$ in Figure 9. During the free cooling phase, which corresponds to regime 1, the cooling luminosity is much larger than the integrated heating rate and the heating does not affect the temperature structure. In regime 2, the planet has cooled such that the integrated heating rate is larger than the internal cooling rate. As a result, the luminosity and temperature profiles are essentially fixed in time at levels shallower than the deposition pressure $P_{\mathrm{dep}}$. In this regime, the cooling itself mostly occurs at the RCB, which corresponds to the time-evolution of $L_{\mathrm{cool}}$ shown in Figure 7. This cooling forces the RCB to move inward with time. In the case with deeper $P_{\mathrm{dep}}=100\ \mathrm{bars}$ there is an outer convective zone and hence this cooling occurs at the inner RCB. Just above this inner RCB, there is a deep near-isothermal region, which corresponds to the location of a radiative window. This radiative window occurs just below the heating level, and the temperature at the top of the radiative window is determined by the integrated heating rate. Hence, even though the integrated cooling rates after $5\ \mathrm{Gyr}$ for the cases with $P_{\mathrm{dep}}=10$ and $100\ \mathrm{bars}$ are similar, the radius for $P_{\mathrm{dep}}=100\ \mathrm{bars}$ is larger. The larger radius is due to the much higher temperature at the bottom of the innermost radiative zone, which leads to a hotter internal adiabat.

5. Comparison with Analytic Theory

In this section, we compare our numerical models to the analytic theory of Ginzburg & Sari (2015). They related the quantity $\Gamma\tau_{\mathrm{dep}}$ to the radius of a hot Jupiter given specified external parameters (equilibrium temperature, planet mass, composition), as it is linked to the central temperature. Specifically, their prediction for radius given $\Gamma\tau_{\mathrm{dep}}$ (equivalent to Equation (32) of Ginzburg & Sari, 2015) is

[TABLE]

where $R_{0}$ is the radius without extra heat deposition ($\approx 1.1\ R_{\mathrm{Jup}}$), $\Delta R_{0}=0.3R_{\mathrm{Jup}}$, $\tau_{\mathrm{dep}}$ the optical depth at which heat is deposited, and $\eta$ and $\delta$ fitting parameters. Note that this prediction is for the radius of a planet with a given internal luminosity, while our numerical simulations calculate the radius at a given age. Ginzburg & Sari (2015) make no explicit prediction for radius at a given age, instead predicting radius for a given cooling luminosity to compare with the numerical calculations of Spiegel & Burrows (2013). Despite the fact that Ginzburg & Sari (2015) make no explicit prediction for the radius as a function of time, here we test the utility of their scaling $\Delta R\propto T_{\mathrm{c}}\propto\left(1+\Gamma\tau_{\mathrm{dep}}/L_{\mathrm{irr}}\right)^{\delta}$ as motivated by their Equations (23) and (29).

The power-law exponent $\delta$ can be determined analytically from the use of a power-law opacity and relationship between radiation energy density and optical depth (see Equation 31 of Ginzburg & Sari, 2015), but here we treat it as a free parameter due to the use of full opacity tables. Ginzburg & Sari (2015) predicted that $\delta\approx 0.19$ to calculate radius at a given age (see their Equation 23), while $\delta\approx 0.15$ to calculate radius at a given internal luminosity.

Figure 10 compares Equation (11) for various $\eta$ and $\delta$ to our numerically determined radius-$\Gamma\tau_{\mathrm{dep}}$ relationship. First, note that as expected from Equation (11), the radius starts to deviate sharply from $R_{0}$ in our numerical solutions when $\Gamma\tau_{\mathrm{dep}}/L_{\mathrm{irr}}\geq 1$. When $\Gamma\tau_{\mathrm{dep}}<L_{\mathrm{irr}}$, the heat source does not strongly perturb the radiative-convective solution due to the relative dominance of the external irradiation. With $\Gamma\tau_{\mathrm{dep}}\geq L_{\mathrm{irr}}$, one can think of the heat source as an effective increase of the incident stellar power, which then perturbs the radiative-convective solution.

In the regime with $\Gamma\tau_{\mathrm{dep}}\geq L_{\mathrm{irr}}$, the analytic prediction from Equation (11) is that the effect of heating increases with increasing $\tau^{\delta}_{\mathrm{dep}}$ and hence approximately increases as $P^{\delta}_{\mathrm{dep}}$. In the case with $P_{\mathrm{dep}}\leq 100\ \mathrm{bars}$, we find that as predicted by Ginzburg & Sari (2015) there is a nearly universal relation between radius and $\Gamma\tau_{\mathrm{dep}}$. We also find that our numerical radius-$\Gamma\tau_{\mathrm{dep}}$ relationship can be reproduced using a constant $\delta\approx 0.1$, which is expected given that $\delta$ is prescribed by the opacity-pressure-temperature relationship (Ginzburg & Sari, 2015). Note that our effective $\delta$ is somewhat lower than the $\delta\approx 0.19$ predicted by Ginzburg & Sari (2015) for comparison at equal ages. This is because of the use of full MESA opacities rather than prescribing a power-law opacity.

However, there is a deviation between our numerical models and the predictions of Ginzburg & Sari (2015) in the case of deep heating at $P_{\mathrm{dep}}>100\ \mathrm{bars}$. In this regime, we need significantly smaller $\eta$ values to match the radius for larger values of $P_{\mathrm{dep}}$333Note that the $\eta$ needed becomes tiny ($\sim 10^{-6}$) in the case with heating at the very center of the planet (not shown). This case was not examined in detail by Ginzburg & Sari (2015), who instead focused on $\tau_{\mathrm{dep}}\lesssim\tau_{\mathrm{RCB}}$. The case with heating at the center of the planet does, however, have the same $\delta$, showing that the opacity-pressure-temperature relationship still determines the power of the relationship between $\Gamma\tau_{\mathrm{dep}}$ and radius even with heating at the very center of the planet.. Thus, the effects of increased heating depths on radius inflation are modest for $P_{\mathrm{dep}}>100\ \mathrm{bars}$. This finding is consistent with our main conclusions, which further show a complete independence of radius inflation on $P_{\mathrm{dep}}\geq 10^{4}\ \mathrm{bars}$. Note that this deviation between simulations and theory is not necessarily a disagreement, as Equation (11) is not directly applicable deeper than the radiative-convective boundary. In this case where $\tau_{\mathrm{dep}}>\tau_{\mathrm{RCB}}$, it is required that $\Gamma>L_{\mathrm{cool}}$ in order for heating to significantly affect the radius (Guillot & Showman, 2002; Wu & Lithwick, 2013; Ginzburg & Sari, 2015). However, Ginzburg & Sari (2015) did not develop an analytic theory to compare with numerical results with $P_{\mathrm{dep}}>100\ \mathrm{bars}$, because these results did not exist at the time. We thus show that their model is very useful for shallow heating, but should not be used for a quantitative prediction of radius at a given age for heat deposited deeper than $100\ \mathrm{bars}$.

6. Conclusions

To better understand the transit radii of hot Jupiters, we studied the evolution of a giant planet subject to intense stellar irradiation at the surface and internal heating. To do so, we used MESA to compute a grid of evolutionary models in which we varied the amount and depth of internal heating using the parameters of a typical hot Jupiter, HD 209458b. To interpret our numerical results, we developed a framework to understand the different cooling regimes of internally heated hot Jupiters. Based on previous work, it is known that the radii of hot Jupiters is larger for more intense or deeper internal heating, provided that this heating is applied throughout the planets’ evolution. Our results are broadly consistent with this expectation, but show that this trend with heating depth, while monotonic, is very uneven. Specifically, we find that:

A minimum heating depth of $100\ \mathrm{bars}$ is required to explain inflated hot Jupiter radii, assuming modest internal heating rates $\lesssim 1\%$ of the stellar irradiation. This pressure level is deeper than the photosphere yet is overall very shallow, enclosing $99.99\%$ of the mass of the planet. We show that modest heating at $\gtrsim 100\ \mathrm{bars}$ is deep enough to enhance the radius for two reasons. First, this heating lies within the convective interior for the initial $10-100\ \mathrm{Myr}$ of evolution. Cooling is thereby significantly suppressed during this crucial early stage. Second, at later times a radiative window opens below the heating layer. The cooling rate of the planet, set at the base of this radiative window, is very low. This enables the planet to retain a much larger radius than for shallower heating. 2. 2.

Heating applied at any depth $\gtrsim 10^{4}\ \mathrm{bars}$ yields nearly identical levels of radius inflation. Since $\approx 99\%$ of the mass of a hot Jupiter lies below $10^{4}\ \mathrm{bars}$, it is remarkable that deeper heating – even heating at the center – produces a similar cooling history. The reason for this depth independence is that cooling timescales from these deep pressures exceed the several Gyrs lifetime of the planet.

Since our models are agnostic as to the mechanism of internal heating, they can be used to constrain a range of hot Jupiter heating mechanisms. Most notably, because we find that at minimum modest heating must occur at $\sim 100\ \mathrm{bars}$ to explain radius inflation, even relatively shallow hydrodynamic inflation mechanisms can explain the transit radii of many hot Jupiters.

We thank Sivan Ginzburg, Bill Hubbard, Adam Showman, and the Steward Observatory Planet Theory group for insightful discussions. We thank the anonymous referee for a thorough review, which improved the manuscript. We thank the MESA team for making this valuable tool publicly available. This research has made use of the Exoplanet Orbit Database and the Exoplanet Data Explorer at exoplanets.org. T.D.K. acknowledges support from NASA headquarters under the NASA Earth and Space Science Fellowship Program Grant PLANET14F-0038. A.N.Y. acknowledges support from the NASA ATP program through grant NNX16AB26G.