Structure and Evolution of Internally Heated Hot Jupiters
Abstract
Hot Jupiters receive strong stellar irradiation, producing equilibrium temperatures of . Incoming irradiation directly heats just their thin outer layer, down to pressures of 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 – does not significantly suppress cooling, unless the total heating rate is of the incident stellar power. Deeper heating, at bars, requires heating at only of the stellar irradiation to explain the observed transit radius of 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 – which is exterior to 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.
Subject headings:
methods: numerical  planets and satellites: gaseous planets  planets and satellites: atmospheres  planets and satellites: interiors  planets and satellites: individual (HD 209458b)1. Introduction
Since the first transit detections of an extrasolar planet (Charbonneau et al., 2000; Henry et al., 2000), the closein 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 highentropy initial state (for reviews see Fortney et al., 2010; Baraffe et al., 2010, 2014; Laughlin & Lissauer, 2015). A massradius 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 overplotted. Approximately onehalf 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 radiativeconvective boundary (RCB) to higher pressures and thereby reduces the planetary cooling rate, this only affects radii by (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 have radii which match the expectations of purely cooling models (Demory & Seager, 2011; Miller & Fortney, 2011). Reinflated 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 reinflated hot Jupiters around postmain 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 stellarflux 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 doublediffusive 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 highentropy 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 eastwest winds are driven by the large daytonight 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 followup 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 selfconsistent 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 vertical winds in hot Jupiter atmospheres transport energy downward. This energy is then dissipated in shear layers near the RCB through, for example, KelvinHelmholtz 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 twodimensional steadystate dynamical model that the largescale 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 timedependent simulations of hot Jupiter atmospheres have yet been performed to quantify the threedimensional 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 longtimescale 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 selfsimilar 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 radiativeconvective 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 radiativeconvective 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:
(1) 
(2) 
(3) 
(4) 
Equation (1) is the mass conservation equation giving the enclosed mass, , at radius, , with mass density . Equation (2) expresses hydrostatic equilibrium of the pressure, , with the gravitational constant, . Equation (3) states energy conservation, where the outgoing luminosity, , computed by MESA includes both radiative () and (where convectively unstable) convective () components. The relevant energy sources include the cooling term, , i.e. the loss of entropy, , that drives gravitational KelvinHelmholz contraction, and , the “extra” energy deposition described in detail below. Equation (4) is the energy transport equation, in terms of the logarithmic gradient, , where is temperature. MESA uses mixing length theory to calculate , the radiative gradient, (or , the convective gradient) in radiative (or convective) regions (respectively). Specifically, the Schwarzchild criterion sets the temperature gradient to the smaller of the adiabatic gradient or the radiative gradient
(5) 
In Equation (5), is the StefanBoltzmann constant and 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 , that assume a dustfree Solar composition (Freedman et al., 2008; Paxton et al., 2011, 2013).
2.2. Numerical Model
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 to compute an adiabatic starting planet model. The longterm evolution, at times , is independent of the initial radius and hence is also independent of the initial entropy (Arras & Bildsten, 2006). However, note that the KelvinHelmholtz 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 presentday transit radius of HD 209458b is . 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 to outgoing radiation is . 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 in Section 2.2.3, and our parameter choices in Section 2.2.4.
Irradiation
The irradiation of the exoplanet is included as an energy generation rate , which is applied where the outer mass column is less than the chosen (see Section 2.2.4 for specific parameter values). With as the incoming stellar flux, this heating drives a heat flux of , equivalent to the average of the incoming radiation over the planet’s spherical surface. This is the 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 after 100 Myr of evolution).
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 nearphotospheric levels to greater depths, we do not consider the impacts of specific heating mechanisms. Instead, we add an extra dissipation (computed each time step) that is taken to be a Gaussian with standard deviation , where 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 (distributed as described above), at the center of the planet, and at the
inner RCB. The integrated heating rates,
(6) 
are set to different fractions of the irradiation as
(7) 
for this grid of simulations. Heating rates are written as a fraction of the incident stellar power relevant for HD 209458b.
We choose an upper limit of because a conversion of starlight to kinetic energy, which is transported to depth via either hydrodynamic or magnetohydrodynamic mechanisms, is unrealistic. Thus we choose 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 can be motivated from the fraction of atmospheric kinetic energy dissipated to heat, which is unknown for hot Jupiters but is expected to be based on Earth (Peixoto & Oort, 1992; Guillot & Showman, 2002; Schubert & Mitchell, 2013). It is expected that for Ohmic dissipation (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 (Arras & Socrates, 2010) and may be as small as (Bodenheimer et al., 2001; Jackson et al., 2008). When we apply heating at the RCB, follows the location of the boundary. Heating at this boundary can be motivated by either downward mixing of heat by largescale 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).
Parameter Choices
We keep all external planetary parameters constant, running a suite of simulations varying only the integrated heating rates and deposition pressure . We use a mass (), composition (), 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 , corresponding to an opacity to incoming radiation, as in
Fortney et al. (2008); Guillot (2010); Owen & Wu (2016).
We use (where is Earth’s incident flux), which is equivalent to an equilibrium temperature 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 the irradiation power ). 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 declines, it eventually enters regime 1, where external irradiation exceeds the planet’s cooling luminosity (i.e. ). 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 ), is the most relevant for our study. The transition from regime 1 to 2 occurs naturally due to the gradual decline in 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 subregimes, 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 . Only hot Jupiters in regime 2(a) can reach an exact steady state with 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 which means the radiative luminosity along the adiabat at the location of applied heating.
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.
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 powerlaw in optical depth with a cutoff above a specified optical depth value. However, a comparison is still possible as the steepness of their powerlaw 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 nonexistent 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 selfsimilar 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 heating efficiency at bars case in Section 4) or never extends close to the center (e.g. our at 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 KelvinHelmholtz 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, , relative to the RCB depth, , at this moment. With “deep heating” and regime 2 begins as 2(b) [or 2(a) if heating is precisely at the center]. With “shallow heating” 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 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 temperaturepressure profile is given by , 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 , 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 timederivative of the characteristic internal entropy, assuming heating in the internal convective region:
(8) 
The internal entropy always decreases with time, as . However, if , the entropy decrease will be effectively zero. This leads to a steady state where 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 and corresponding decrease in . However, cannot be deeper than in regime 2(b), by definition, and cannot drop below the limit described above, i.e. . Furthermore, regime 2(b) should not simply evolve into regime 2(d) because with , 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 , 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 reinflation, 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 radius
First, there is a large increase in transit radius between the cases with heating at and . Secondly, the model transit radii are essentially the same for all of the cases with deeper heating, namely with heating either at , heating at the center, or heating at the RCB. Figure 4 shows that heating at still encloses of the mass of the planet. Even heating at , which is in the inner convective region, encloses 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
Structure: Varying Depth of Heating
The trends in transit radius with varying and 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 show almost the same radius inflation for a given heating rate . Secondly, there is a large () jump in radius between . 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 pressuretemperature profiles at the endstate of evolution for varying with a fixed normalized heating rate . These profiles show a clear bifurcation between shallow heating () and deep heating (). Simulations with shallow heating at have a deep outer radiative zone, extending to . By contrast, with deep heating at the outer radiative zone only extends to . 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 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 for simulations with the normalized heating rate . The two cases of and illustrate the transition from shallow to deep heating, respectively. The profiles in Figure 6 show that heating at triggers an outer convective region, with an inner radiative window below. Since , the same heating () at a deeper pressure forces a secondary convective region.
For heating at , the increase in 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):
(9) 
The quantity for and for , which is hence just shallow enough to not force a secondary convective zone. As a result, the simulations with match onto a lower internal adiabat, naturally explaining the smaller transit radius.
Evolution: Structural Quantities
Thus far we have only presented the endstate 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 and . Here we investigate the evolution of two slices of our grid, varying with fixed normalized heating rate . Figure 7 displays the timeevolution of photospheric radius, central entropy, the gravitational cooling luminosity
(10) 
and the pressure at the inner RCB for and .
The early evolution of these structural quantities is similar for both and . All models show an identical cooling phase over the first , which extends to with weaker heating of . 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 ) or deep ().
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, , the planet next reaches a quasiequilibrium 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 , and for and in this state, respectively. Thus the hot Jupiters are mostly in regime 2(b) during this quasiequilibrium 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 in which heating supplies the total interior luminosity. This behavior is evident in the sharp dropoff in , see the third row of Figure 7.
For the cases with deep atmospheric heating at , the quasiequilibrium 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 quasiequilibrium roughly corresponds to a sharp increase in , i.e. the opening of a radiative window which marks the transition from regime 2(b) to 2(c). When the radiative window opens, does not increase, but it stops declining rapidly with time. This more slowly evolving is sufficient to eventually cause a noticeable decrease in planetary entropy and radius.
The case with is notable because the evolution in radius and central entropy is nearly indistinguishable from the case with central heating. The behavior of is noticeably different for the two heating locations, with dropping to arbitrarily low levels at late times with applied central heating. However, with , stays just above , for both and . Crucially, this finite level of cooling corresponds to KelvinHelmholtz cooling times of , much longer than the system age. Thus, departures from a true equilibrium with no cooling are hard to discern.
A subtle point about the case is that the floor in at late times has a different origin in the and cases, despite their similar value of . For , the plateau in occurs when a radiative window opens. This window opens at , i.e. of the final age. However, the window is so deep, approaching a pressure of , 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 , a radiative window does not open. In this case the floor in corresponds to the minimum value of discussed in Section 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 luminositypressure profiles for a normalized heating rate and varying deposition pressure at fixed central entropy , which is approximately the equilibrium entropy for central heating at this value of . 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 , a secondary rise to the heating level , and a third to the irradiation level . The profiles in the outer two levels are almost the same for all , 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 (at fixed central entropy), due to the longer cooling timescales at greater depths.
To examine further the differences between the cases with and , we show the timeevolution of their luminosity and temperature profiles with 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 . In this regime, the cooling itself mostly occurs at the RCB, which corresponds to the timeevolution of shown in Figure 7. This cooling forces the RCB to move inward with time. In the case with deeper there is an outer convective zone and hence this cooling occurs at the inner RCB. Just above this inner RCB, there is a deep nearisothermal 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 for the cases with and are similar, the radius for 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 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 (equivalent to Equation (32) of Ginzburg & Sari, 2015) is
(11) 
where is the radius without extra heat deposition (), , the optical depth at which heat is deposited, and and 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 as motivated by their Equations (23) and (29).
The powerlaw exponent can be determined analytically from the use of a powerlaw 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 to calculate radius at a given age (see their Equation 23), while to calculate radius at a given internal luminosity.
Figure 10 compares Equation (11) for various and to our numerically determined radius relationship. First, note that as expected from Equation (11), the radius starts to deviate sharply from in our numerical solutions when . When , the heat source does not strongly perturb the radiativeconvective solution due to the relative dominance of the external irradiation. With , one can think of the heat source as an effective increase of the incident stellar power, which then perturbs the radiativeconvective solution.
In the regime with , the analytic prediction from Equation (11) is that the effect of heating increases with increasing and hence approximately increases as . In the case with , we find that as predicted by Ginzburg & Sari (2015) there is a nearly universal relation between radius and . We also find that our numerical radius relationship can be reproduced using a constant , which is expected given that is prescribed by the opacitypressuretemperature relationship (Ginzburg & Sari, 2015). Note that our effective is somewhat lower than the predicted by Ginzburg & Sari (2015) for comparison at equal ages. This is because of the use of full MESA opacities rather than prescribing a powerlaw opacity.
However, there is a deviation between our numerical models and the predictions of Ginzburg & Sari (2015) in the case of deep heating at . In this regime, we need significantly smaller values to match the radius for larger values of
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 is required to explain inflated hot Jupiter radii, assuming modest internal heating rates of the stellar irradiation. This pressure level is deeper than the photosphere yet is overall very shallow, enclosing of the mass of the planet. We show that modest heating at is deep enough to enhance the radius for two reasons. First, this heating lies within the convective interior for the initial 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.

Heating applied at any depth yields nearly identical levels of radius inflation. Since of the mass of a hot Jupiter lies below , 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 to explain radius inflation, even relatively shallow hydrodynamic inflation mechanisms can explain the transit radii of many hot Jupiters.
Footnotes
 slugcomment: Accepted at ApJ
 This 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.
 To 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 .
 Note that the needed becomes tiny () 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 . The case with heating at the center of the planet does, however, have the same , showing that the opacitypressuretemperature relationship still determines the power of the relationship between and radius even with heating at the very center of the planet.
References
 Arras, P. & Bildsten, L. 2006, The Astrophysical Journal, 650, 394
 Arras, P. & Socrates, A. 2010, The Astrophysical Journal, 714, 1
 Baraffe, I., Chabrier, G., & Barman, T. 2010, Reports on Progress in Physics, 76, 30
 Baraffe, I., Chabrier, G., Fortney, J., & Sotin, C. Protostars and Planets VI, ed. H. Beuther, R. Klessen, C. Dullemond, & T. Henning (Tucson, AZ: University of Arizona Press)
 Batygin, K. & Stevenson, D. 2010, The Astrophysical Journal Letters, 714, L238
 Batygin, K., Stevenson, D., & Bodenheimer, P. 2011, The Astrophysical Journal, 738, 1
 Bodenheimer, P., Laughlin, G., Rozyczka, M., & Yorke, H. 2007, Numerical Methods in Astrophysics (Boca Raton, FL: Taylor & Francis)
 Bodenheimer, P., Lin, D., & Mardling, R. 2001, The Astrophysical Journal, 548, 466
 Burrows, A., Hubeny, I., Budaj, J., & Hubbard, W. 2007, The Astrophysical Journal, 661, 502
 Chabrier, G. & Baraffe, I. 2007, The Astrophysical Journal Letters, 661, 81
 Chandrasekhar, S. 1939, An Introduction to Stellar Structure (Chicago, IL: University of Chicago Press)
 Charbonneau, D., Brown, T., Latham, D., & Mayor, M. 2000, The Astrophysical Journal, 529, L45
 Cooper, C. & Showman, A. 2005, The Astrophysical Journal Letters, 629, L45
 Demory, B. & Seager, S. 2011, The Astrophysical Journal Supplement Series, 197, 12
 Fortney, J., Baraffe, I., & Militzer, B. Exoplanets, ed. S. Seager (Tucson, AZ: University of Arizona Press)
 Fortney, J., Lodders, K., Marley, M., & Freedman, R. 2008, The Astrophysical Journal, 678, 1419
 Fortney, J., Marley, M., & Barnes, J. 2007, The Astrophysical Journal, 659, 1661
 Freedman, R., Marley, M., & Lodders, K. 2008, The Astrophysical Journal Supplement Series, 174, 504
 Ginzburg, S. & Sari, R. 2015, The Astrophysical Journal, 803, 111
 —. 2016, The Astrophysical Journal, 819, 116
 Grunblatt, S. K., Huber, D., Gaidos, E. J., Lopez, E. D., Fulton, B. J., Vanderburg, A., Barclay, T., Fortney, J. J., Howard, A. W., & Isaacson, H. T. 2016, The Astronomical Journal, 152, 1
 Guillot, T. 2010, Astronomy and Astrophysics, 520, A27
 Guillot, T., Burrows, A., Hubbard, W., Lunine, J., & Saumon, D. 1996, The Astrophysical Journal Letters, 459, L35
 Guillot, T., Gautier, D., Chabrier, G., & Mosser, B. 1994, Icarus, 112, 337
 Guillot, T. & Showman, A. 2002, Astronomy and Astrophysics, 385, 156
 Han, E., Wang, S., Wright, J., Feng, Y., Zhao, M., Fakhouri, O., Brown, J., & Hancock, C. 2014, Publications of the Astronomical Society of the Pacific, 126, 827
 Hartman, J. D., Bakos, G. Á., Bhatti, W., Penev, K., Bieryla, A., Latham, D. W., Kovács, G., Torres, G., Csubry, Z., de ValBorro, M., Buchhave, L., Kovács, T., Quinn, S., Howard, A. W., Isaacson, H., Fulton, B. J., Everett, M. E., Esquerdo, G., Béky, B., Szklenar, T., Falco, E., Santerne, A., Boisse, I., Hébrard, G., Burrows, A., Lázár, J., Papp, I., & Sári, P. 2016, The Astronomical Journal, 152, 182
 Heng, K. 2012, The Astrophysical Journal Letters, 748, L17
 Henry, G., Marcy, G., Butler, R., & Vogt, S. 2000, The Astrophysical Journal, 529, L41
 Henyey, L., Wilets, L., Bohm, K., LeLevier, R., & Levee, R. 1959, The Astrophysical Journal, 129, 628
 Huang, X. & Cumming, A. 2012, The Astrophysical Journal, 757, 47
 Ibgui, L. & Burrows, A. 2009, The Astrophysical Journal, 700, 1921
 Ibgui, L., Burrows, A., & Spiegel, D. 2010, The Astrophysical Journal, 713, 751
 Jackson, B., Greenberg, R., & Barnes, R. 2008, The Astrophysical Journal, 681, 1631
 Kataria, T., Sing, D., Lewis, N., Visscher, C., Showman, A., Fortney, J., & Marley, M. 2016, The Astrophysical Journal, 821, 9
 Kippenhahn, R., Weigert, A., & Weiss, A. 2012, Stellar Structure and Evolution, 2nd edn. (New York: Springer)
 Komacek, T. & Showman, A. 2016, The Astrophysical Journal, 821, 16
 Komacek, T., Showman, A., & Tan, X. 2017, The Astrophysical Journal, 835, 198
 Laughlin, G., Crismani, M., & Adams, F. 2011, The Astrophysical Journal Letters, 729, L7
 Laughlin, G. & Lissauer, J. Treatise on Geopysics, 2nd edn., ed. G. Schubert (Elsevier)
 Leconte, J. & Chabrier, G. 2012, Astronomy and Astrophysics, 540, A20
 Leconte, J., Chabrier, G., Baraffe, I., & Levrard, B. 2010, Astronomy and Astrophysics, 516, A64
 Lopez, E. & Fortney, J. 2016, The Astrophysical Journal, 818, 4
 Mayne, N., Baraffe, I., Acreman, D., Smith, C., Browning, M., Amundsen, D., Wood, N., Thuburn, J., & Jackson, D. 2014, Astronomy and Astrophysics, 561, A1
 Menou, K. 2012a, The Astrophysical Journal, 745, 138
 —. 2012b, The Astrophysical Journal Letters, 754, L9
 Menou, K. & Rauscher, E. 2009, The Astrophysical Journal, 700, 887
 Miller, N. & Fortney, J. 2011, The Astrophysical Journal Letters, 736, L29
 Miller, N., Fortney, J., & Jackson, B. 2009, The Astrophysical Journal, 702, 1413
 Owen, J. & Wu, Y. 2016, The Astrophysical Journal, 817, 107
 Paxton, B., Bildsten, L., Dotter, A., Herwig, F., Lesaffre, P., & Timmes, F. 2011, The Astrophysical Journal Supplement Series, 192, 3
 Paxton, B., Cantiello, M., Arras, P., Bildsten, L., Brown, E., Dotter, A., Mankovich, C., Montgomery, M., Stello, D., Timmes, F., & Townsend, R. 2013, The Astrophysical Journal Supplement Series, 208, 4
 Paxton, B., Marchant, P., Schwab, J., Bauer, E., Bildsten, L., Cantiello, M., Dessart, L., Farmer, R., Hu, H., Langer, N., Townsend, R., Townsley, D., & Timmes, F. 2015, The Astrophysical Journal Supplement Series, 220, 15
 Peixoto, J. & Oort, A. 1992, Physics of Climate (New York: American Institute of Physics)
 Perna, R., Heng, K., & Pont, F. 2012, The Astrophysical Journal, 751, 59
 Perna, R., Menou, K., & Rauscher, E. 2010, The Astrophysical Journal, 724, 313
 Rauscher, E. & Menou, K. 2010, The Astrophysical Journal, 714, 1334
 —. 2013, The Astrophysical Journal, 764, 103
 Rogers, T. & Komacek, T. 2014, The Astrophysical Journal, 794, 132
 Rogers, T. & Showman, A. 2014, The Astrophysical Journal Letters, 782, L4
 Saumon, D., Chabrier, G., & Horn, H. V. 1995, The Astrophysical Journal Supplement Series, 99, 713
 Schubert, G. & Mitchell, J. Comparative Climatology of Terrestrial Planets, ed. S. Mackwell, A. SimonMiller, J. Harder, & M. Bullock (Tucson, AZ: University of Arizona Press)
 Showman, A., Cho, J., & Menou, K. Exoplanets, ed. S. Seager (Tucson, AZ: University of Arizona Press)
 Showman, A., Fortney, J., Lian, Y., Marley, M., Freedman, R., Knutson, H., & Charbonneau, D. 2009, The Astrophysical Journal, 699, 564
 Showman, A. & Guillot, T. 2002, Astronomy and Astrophysics, 385, 166
 Showman, A. & Polvani, L. 2011, The Astrophysical Journal, 738, 71
 Spiegel, D. & Burrows, A. 2013, The Astrophysical Journal, 772, 76
 Tremblin, P., Chabrier, G., Mayne, N., Amundsen, D., Baraffe, I., Debras, F., Drummond, B., Manners, J., & Fromang, S. 2017, The Astrophysical Journal, 843, 30
 Valsecchi, F., Rappaport, S., Rasio, F., Marchant, P., & Rogers, L. 2015, The Astrophysical Journal, 813, 101
 Weiss, L., Marcy, G., Rowe, J., Howard, A., Isaacson, H., Fortney, J., Miller, N., Demory, B., Fischer, D., Adams, E., Dupree, A., Howell, S., Kolbl, R., Johnson, J., Horch, E., Everett, M., Fabrycky, D., & Seager, S. 2013, The Astrophysical Journal, 768, 14
 Wu, Y. & Lithwick, Y. 2013, The Astrophysical Journal, 763, 13
 Youdin, A. & Mitchell, J. 2010, The Astrophysical Journal, 721, 1113