MNRAS 469, 1768–1782 (2017) doi:10.1093/mnras/stx831
Advance Access publication 2017 April 5
Tidal heating and stellar irradiation of hot Jupiters
Adam S. Jermyn,1‹ Christopher A. Tout1 and Gordon I. Ogilvie2
1Institute of Astronomy, University of Cambridge, Madingley Rd, Cambridge CB3 0HA, UK
2Department of Applied Mathematics and Theoretical Physics, University of Cambridge, Madingley Rd, Cambridge CB3 0WA, UK
Accepted 2017 March 31. Received 2017 March 30; in original form 2016 November 9
ABSTRACT
We study the interaction between stellar irradiation and tidal heating in gaseous planets
with short orbital periods. The intentionally simplified atmospheric model we employ makes
the problem analytically tractable and permits the derivation of useful scaling relations. We
show that many tidal models provide thermal feedback, producing interior radiative zones
and leading to enhanced g-mode dissipation with a wide spectrum of resonances. These
resonances are dynamically tuned by the thermal feedback, and so represent a novel form of
thermomechanical feedback, coupling vibrational modes to the very slow thermal evolution
of the planet. We then show that stellar irradiation allows the heat produced by these modes to
be trapped at depth with high efficiency, leading to entropy increase in the central convective
region, as well as expansion of the planet’s radius sufficient to match observed swelling.
We find that thermally driven winds play an essential role in this process by making the
thermal structure of the atmosphere spherically symmetric within a few scale heights of the
photosphere. We characterize the relationship between the swelling factor, the orbital period
and the host star and determine the time-scale for swelling. We show that these g modes
suffice to produce bloating on the order of the radius of the planet over Gyr time-scales
when combined with significant insolation and we provide analytic relations for the relative
magnitudes of tidal heating and insolation.
Key words: planets and satellites: interiors – planet–star interactions.
1 IN T RO D U C T I O N
In recent data from Kepler and ground-based follow-ups, there is
evidence for a large population of hot Jupiters that are substan-
tially inflated relative to their degenerate radii (Hartman et al. 2012;
Hellier et al. 2012; Weiss et al. 2013). The radii and periods of
known members of this population and of the broader Jupiter-sized
population are shown in Fig. 1 (Rein 2012). There is an apparent
split in the observed population around periods of 10 d, such that
planets with longer periods are generally not inflated while those
with shorter periods are often substantially inflated. Importantly,
planets at or above 2RJ must be inflated relative to their degener-
ate radii, otherwise their implied masses would make them stars
(Stevenson 1991). Importantly, Jupiter has approximately the max-
imum radius for an unheated gas giant, so planets at 2RJ must be
bloated regardless of their mass. In order to achieve this level of
expansion, the central convection zone must be heated considerably
relative to what would be expected as a result of the residual heat of
formation (Lopez & Fortney 2016), and there is evidence of planets
re-inflating after cooling down (Hartman et al. 2016). Complicating
 E-mail: adamjermyn@gmail.com
this is the thermodynamic requirement that heat flows only from hot
to cold, not in the reverse fashion, and in the absence of an internal
heat source there is no means for blocked heat transport to heat
the interior. This, combined with the expectation that temperature
increases towards the core of the planet, means that any change in
the temperature at depth must be due to heat generated at or deeper
than the point of interest.
A variety of mechanisms have been suggested to gener-
ate deep heating, with ohmic processes (Batygin, Stevenson &
Bodenheimer 2011; Spiegel & Burrows 2013) and tidal dissipa-
tion (Miller, Fortney & Jackson 2009; Youdin & Mitchell 2010;
Socrates 2013) among the more popular models. Given that heat
cannot flow from the surface of the planet into its core, the stellar
flux is often neglected. However, somewhat surprisingly, the ob-
served radii correlate strongly with the incident stellar flux, so that
this flux may play a role in the inflation process (Lopez & Fort-
ney 2016). Confounding this analysis is the fact that stellar flux is
not independent of the orbital period. So a theory of hot Jupiter
inflation must separately handle the effects of orbital period and
incident flux, particularly when dealing with tidal heating.
We investigate the effects of stellar flux on the structure of an
internally heated hot Jupiter, making few assumptions about the
nature or profile of the heating and considering the effects of wind
C© 2017 The Authors
Published by Oxford University Press on behalf of the Royal Astronomical Society
Tidal heating of hot Jupiters 1769
Figure 1. Known population of short-period planets with radii near that of
Jupiter.
redistribution. We show that the stellar flux acts to modulate the
rate at which heat escapes from the planet. We then investigate the
feedback that this heating produces on the thermal structure of the
planet and show that a wide variety of realistic heating profiles
gives rise to interior radiative zones. These zones migrate within
the planet on thermal time-scales, giving a broad and dynamically
tuned spectrum of g-mode resonances that dissipate heat tidally in
the planet. We then show that these modes suffice to produce the
observed bloating. Finally, we predict the relation between stellar
flux, orbital period and planetary radius.
The new thermomechanical feedback mechanism we propose,
shown schematically in Fig. 2, underscores the importance of con-
sidering planets as dynamical objects with complex behaviours cou-
pling wildly different time-scales. Vibrational effects with periods
ranging from seconds to days can have a tremendous impact on
thermal evolution over millions of years, and that thermal evolution
in turn feeds back into the vibrational modes, creating a dynami-
cally tuned spectrum that can ultimately determine the large-scale
structure of the planet.
2 ISOTROPIC P L ANETA RY ST RUCTURE
We discuss the structure of a planet isotropically illuminated by
flux Fe from its host star. Fig. 3 shows the orbital configura-
tion of the planet–star system, and Fig. 4 shows the thermody-
namic structure of the planet with the relevant variables defined
schematically.
For the deep interior of the planet we adopt the analytic brown
dwarf structure of Stevenson (1991):
ψ ≡ kBT
EF
= 8 × 10−6μ2/3e
(
ρ
g cm−3
)−2/3 (
T
K
)
, (1)
R0 = 2.8 × 109 cm
(
M
M
)−1/3
μ−5/3e , (2)
r = R0
(
1 + ψ + ψ
2
1 + ψ
)
, (3)
Figure 2. Schematic of the proposed thermomechanical feedback mecha-
nism. The upper convective layer (blue), radiative layer (beige) and inner
convective layer (yellow) are shown as concentric shells. The boundaries of
the radiative layer are moving inwards at different rates, allowing the zone
to resize. Profiles of g modes (dark blue) are shown along the equator and
schematically depicted at other latitudes.
Figure 3. Orbital configuration of the planet and its host star.
Figure 4. Thermodynamic structure of the planet. The unperturbed struc-
ture is shown on the left, while the heated (perturbed) structure is shown on
the right.
MNRAS 469, 1768–1782 (2017)
1770 A. S. Jermyn, C. A. Tout and G. I. Ogilvie
P = 1013 erg cm−3 μ−5/3e
(
ρ
g cm−3
)5/3 (
r
R0
)
(4)
and
∇a ≡ ∂ ln T
∂ lnP
∣∣∣∣
s
= 2
5
, (5)
where ψ is the electron degeneracy parameter, ρ is the density, T is
the temperature, R0 is the degenerate radius, r is the radial coordi-
nate, P is the pressure, kB is the Boltzmann constant, EF is the Fermi
energy, M is the mass of the planet and ∇a is the adiabatic temper-
ature gradient. These relations effectively parametrize a γ = 5/3
adiabatic atmosphere, accounting for electron degeneracy at high
pressures. We assume solar composition in this paper, so that the
mean molecular weight of electrons μe ≈ 1.15. In addition, we take
R to be the radius of the planet. For convenience, we define the
parameters
R ≡ R/R0 (6)
and
M ≡ M/MJ, (7)
where MJ = 1.838 × 1030 g is the mass of Jupiter. We expect the
gas line opacity to dominate in hot atmospheres, so we use this as
fiducial and define (Stevenson 1991)
κ0 ≡ 10−2 cm2 g−1. (8)
We connect the top of the convection zone to the photosphere
with a radiative zone at transition pressure Pt. The photospheric
temperature is given by
Tph =
(
F
σ
)1/4
, (9)
where F is the total flux leaving the planet’s atmosphere and σ is
the Stefan–Boltzmann constant. This flux may be divided into two
components, as
F = Fi + Fe, (10)
where Fi is the heat arriving from the planet’s interior and Fe is
the heat arriving from the host star. When Fe is large relative to
the flux that would escape through the planet’s natural cooling, the
photospheric temperature is determined entirely by Fe, such that
Tph =
(
Fe
σ
)1/4
. (11)
The escaping flux Fi from the planet’s core is just the flux that
escapes from the convection zone. This may be generated by grav-
itational contraction, radioactive decay or by a decrease in the in-
terior entropy; here we generally assume the last of these to be the
dominant source of interior flux.
The radiative gradient is given by
∇r = 3κPL64πGMσT 4 , (12)
where L is the luminosity driving the gradient. At the convective–
radiative boundary we have
∇a = ∇r = 3κPtLi64πGMσT 4t
, (13)
where
Li = 4πR2Fi. (14)
If we take κ to be a power law in both P and T of the form
κ = κ0T aP b, (15)
then
∇a = 3κ0P
1+b
t Li
64πGMσT 4−at
. (16)
Above the transition we have
∇r = ∇a
(
P
Pt
)1+b (
T
Tt
)a−4
= d ln T
d lnP
. (17)
Integrating from the photosphere to the transition yields
1 −
(
Tph
Tt
)4−a
= 4 − a
1 + b∇a
(
1 −
(
Pph
Pt
)1+b)
. (18)
The photosphere pressure is generally much lower than the transi-
tion pressure so
Tph
Tt
=
(
1 − 4 − a
1 + b∇a
) 1
4−a
, (19)
which has a solution if and only if
1 − 4 − a
1 + b∇a ≥ 0. (20)
If b >−1 and a > 4 or b <−1 and a < 4 or b > 3−2a5 and a < 4, this
is satisfied. There are other conditions under which it is satisfied,
but these are the most relevant common cases. The exponents are
generally of order unity and the result is raised to a small power so
Tph ≈ Tt. This agrees with other analyses, which have found that
in Jupiter-like planets, Tph < Tt < 1.5Tph (Ginzburg & Sari 2015).
The precise temperature ratio depends on the nature of the opacity
function, so we simply take Tt = 21/4Tph (the Eddington closed grey
body) as representative. If this holds and the flux from the star is
dominant, then
Fi = 16σg∇aT
4
t
3κPt
=
(
32 g∇a
3κPt
)
Fe. (21)
Eliminating μe between equations (1) and (4), we may write Pt in
terms of Tt at r ≈ R and R0, such that
Pt = 1013 erg cm−3
(
R
R0
)(
Tt
K
)5/2 (
ψ
8 × 106
)−5/2
. (22)
From equation (3) and the definition of R0 (equation 6) we derive
ψ = 1
4
⎛
⎝ R
R0
− 2 +
√(
R
R0
)2
+ 4
(
R
R0
)
− 4
⎞
⎠ (23)
= 1
4
(
R− 2 +
√
R2 + 4R− 4
)
, (24)
where we have taken the positive root because ψ > 0 and R > R0.
Note that this is always of order unity, so to good approximation the
majority of the variation in Pt comes from the R and Tt dependence
in equation (22).
If the convective–radiative transition occurs at a shallow point in
the atmosphere the corresponding column density is just
t ≈ Pt
g
= 4 × 109 g cm−2R3M−5/3
(
Tt
K
)5/2 (
ψ
8 × 10−4
)−5/2
,
(25)
MNRAS 469, 1768–1782 (2017)
Tidal heating of hot Jupiters 1771
again with M = M/MJ. Eliminating Tt in favour of Fe and using
T = 5777 K we find
t ≈ 2 × 106 g cm−2R3M−5/3
(
Fe
F
)5/8
ψ−5/2. (26)
Inserting this result into equation (21) we obtain
Fi
Fe
= 32g∇a
3κPt
= 2 × 10−4M5/3
(
Fe
F
)−5/8 (
κ
κ0
)−1
ψ5/2
R3 . (27)
As one final manipulation, we wish to put our equations in terms
of the stellar luminosity and orbital radius. The stellar luminosity is
related to the external flux Fe by
4πR2Fe = πR
2L
4πa2orbit
, (28)
where L is the stellar luminosity, aorbit is the orbital radius of the
planet. The factor of πR2 on the right-hand side is just the cross-
section of the planet as seen from the star, while the factor of 4πR2
on the left-hand side reflects the definition of Fe as an average over
the surface of the planet. So
Fe = L16πa2orbit
. (29)
Comparing with the Sun we find
Fe
F
= 1
4
(
L
L
)(
aorbit
R
)−2
. (30)
Thus
Fi
Fe
= 5 × 10−4M5/3
(
L
L
)−5/8 (
aorbit
R
)5/4 (
κ
κ0
)−1
ψ5/2
R3 .
(31)
Importantly, the exponent on the luminosity is greater than −1. This
means that while the ratio of escaping to incident flux decreases
with increasing stellar flux, the total escaping flux increases. This
conclusion is dependent primarily on how strongly the ratio Tt/Tph
varies with Fe, which in turn depends on the form of the opacity
function. In particular, it does not generally hold at extremely high
temperatures where the gas line opacity ceases to dominate and
Kramers-like rules take over. For brown dwarfs and hot Jupiters,
however, this variation is small and should not pose a problem. It is
also useful to compute the transition column density:
t ≈ Pt
g
= 8 × 105 g cm−2R3M−5/3
×
(
L
L
)5/8 (
aorbit
R
)−5/4
ψ−5/2. (32)
This is small enough that the shallow approximation is not bad.
3 A N G U L A R T E M P E R ATU R E D I S T R I BU T I O N
The planets under consideration are generally highly insolated. This
can lead to significant temperature differences between the day and
night sides, particularly if the planet is tidally locked. In this section
we show that winds suffice to make the thermal structure of the
atmosphere spherically symmetric at depth even when there is a
large temperature difference at the photosphere. This allows us to
treat the structure of the planet as spherically symmetric where tidal
effects are most prominent.
Consider a wind driven from one side of the planet to the other
along isobars with characteristic velocity v. Suppose further that
the character of this wind changes in the vertical direction over
distances of order the pressure scale height h and that it changes in
the horizontal direction over distances of order the planet’s radius.
The specific force due to shear in the vertical direction is
Fv = νv ∂v
∂r
, (33)
where νv is the viscosity for a circumferential flow shearing in the
vertical direction. The corresponding power dissipated is
Pv = ∂v
∂r
· Fv = νv
(
∂v
∂r
)2
. (34)
Likewise, the force due to shear in the horizontal direction is
Fh = νh ∂v
∂ξ
, (35)
where νh is the viscosity for a circumferential flow shearing in the
other circumferential direction and ξ is a coordinate along the flow.
The corresponding power dissipated is
Ph = ∂v
∂ξ
· Fh = νv
(
∂v
∂ξ
)2
. (36)
The total power dissipated is then
P = Pv + Ph ≈ v2
[ νv
h2
+ νh
r2
]
, (37)
where we have approximated the velocity derivatives with the ve-
locity magnitude and the relevant scale heights, the pressure scale
height h in the vertical direction and the radius r in the horizontal
direction. We have also simplified the viscosity from a rank-4 tensor
to two scalars, so this relation ought only to be interpreted as an
order of magnitude of the power.
To determine v, we now match this power to the work that the
wind may extract as a heat engine. We are interested in cases where
the temperature difference between the two sides is large so
the efficiency of the heat engine is of order unity even if diffu-
sive losses make it irreversible. We may neglect diffusive losses
because we have taken the microscopic thermal diffusivity to be
small on the relevant scales. So we may write the specific rate of
work as
W = cpv · ∇T ≈ cpvT
πr
, (38)
where
cp = 5Rgas/2, (39)
for a monatomic ideal gas, is the specific heat at constant pressure.
This is simply the specific heat that is transported from one side of
the planet to the other. In our case
T ≡ Tday − Tnight (40)
and
T ≡ 1
2
(
Tday + Tnight
)
, (41)
so T refers to the average temperature while T refers to the tem-
perature difference. By definition, T/T ≤ 2. In the most extreme
case this gives
W ≈ cpvT
πr
≈ cpT v
πr
(
T
T
)
≈ 5c
2
s v
2γπr
, (42)
where
cs =
√
γRgasT
μ
(43)
MNRAS 469, 1768–1782 (2017)
1772 A. S. Jermyn, C. A. Tout and G. I. Ogilvie
is the adiabatic sound speed and μ is the mean molecular weight.
Equating the rate of work and power gives
c2s =
2γ
5
πrv
[ νv
h2
+ νh
r2
]
. (44)
To proceed further we must examine the forms of νv and νh. The na-
ture of the viscosity differs between stably stratified and buoyantly
unstable zones, so we must determine which of these are relevant
and treat them separately.
We begin with radiative zones. In a stably stratified region the two
viscosities differ because of Richardson stabilization, an effect that
limits the scale of turbulence in the vertical direction by means of a
buoyant restoring force (Galperin, Sukoriansky & Anderson 2007).
A straightforward prescription for the viscosities in this context is
νh ≈ vr, (45)
νv ≈ v2
(
α + νh
gh(∇a − ∇)
)
, (46)
where α is the microscopic thermal diffusivity (Mathis, Palacios
& Zahn 2004). Generally we expect α to be small compared to
νh because horizontal radiative transfer is inefficient, so we may
neglect α and write
νv = v2
(
νh
gh(∇a − ∇)
)
= v
3r
gh(∇a − ∇) . (47)
By the Schwarzschild criterion ∇ < ∇a in a stably stratified zone.
In general we expect radiative transport to be efficient far from the
zone boundaries, so we take ∇ 	 ∇a in most of such a zone. Using
this we write
νv = v
3r
gh∇a . (48)
Now making use of
gh = g
∣∣∣∣ drd lnP
∣∣∣∣ = gP
∣∣∣∣ drdP
∣∣∣∣ = gPgρ = Pρ = γ−1c2s , (49)
we find
νv = v
3rγ
c2s ∇a
. (50)
Inserting equations (45) and (50) into equation (44) gives
c2s =
2γ
5
πrv
[
v3rγ
c2s h
2∇a +
v
r
]
. (51)
This may be rearranged to
5
2πγ
= γ∇a
(
v
cs
)4 ( r
h
)2
+
(
v
cs
)2
. (52)
Solving equation (52) gives
(
v
cs
)2
= ∇ah
2
2γ r2
⎡
⎣−1 ±
√
1 + 10r
2
π∇ah2
⎤
⎦. (53)
The positive branch is the one of interest, because we have implicitly
taken v > 0 in writing it as a magnitude. In the upper regions of the
planet’s atmosphere r 
 h so
v
cs
≈ h
r
√
5
2πγ
. (54)
Using equation (49) the rate at which heat is transported may be
written as
ε =W ≈ 5c
2
s v
2πγ r
≈ c
3
s h
r2
(
5
2πγ
)3/2
≈ γ c
5
s
gr2
(
5
2πγ
)3/2
. (55)
The region of interest is shallow so gr2 ≈ GM and
ε ≈ γ c
5
s
GM
(
5
2πγ
)3/2
. (56)
The depth, as measured by column density i, over which the
winds make the flux distribution spherically symmetric is
i = Fe
ε
≈ GMσT
4
ph
γ c5s
(
2πγ
5
)3/2
. (57)
Evaluating the sound speed at the photosphere gives
i ≈ 3 × 103 g cm−2M
(
T
×103 K
)3/2
, (58)
where mp is the proton mass. For comparison, the photosphere is at
a depth of
ph ≈ κ−1 = 102 g cm−2
(
κ
κ0
)−1
. (59)
Thus the temperature distribution becomes spherically symmetric
deeper than the photosphere but shallower than the convective tran-
sition. So we need not worry about the viscosity in convection
zones.
This remains valid as long as the planet rotates slowly relative
to v/R, such that the characteristic scale of circumferential motion
remains R and is not reduced by Coriolis effects. At short periods,
where this condition is most in danger, the anisotropy is very large,
such that v ≈ hcs/r, and the surface temperature should be quite
high because of insolation, such that cs ≈ 3 × 105 cm s−1. In this
regime, the period of a Jupiter-radius planet must be at least 30 d
with h/r ≈ 10−2 or 3 d with l/r ≈ 10−1 for the Coriolis effect to be
negligible. Even at the shortest known periods of just under a day, the
correction term is not too great and does not alter the conclusion that
the temperature distribution becomes spherical above the convection
zone, so we continue to use this approximation with the knowledge
that it becomes worse as the period diminishes.
4 H E AT E D T H E R M A L ST RU C T U R E
In this section we work on time-scales long compared to the ad-
justment of radiative or convective zones to thermal perturbations
but short compared to the characteristic thermal time-scale of the
planet. This is the instantaneous equilibrium approximation. This
separation of scales exists because the thermal time-scale of the
planet is set by the thermal content of the core, whereas the radia-
tive and convection regions of interest are shallow zones with much
less mass and at much lower temperatures.
The equations governing the luminosity of the planet as a function
of mass coordinate are
∂L
∂m
= ε(m) − cp ∂S
∂t
, (60)
L(0) = 0 (61)
and
L(M) = Li, (62)
MNRAS 469, 1768–1782 (2017)
Tidal heating of hot Jupiters 1773
where ε(m) is the specific energy generation by tides, radioactive
decay and ohmic processes and the mass coordinate m corresponds
to the spherical shell containing mass m. Note that mechanical ex-
pansion and contraction can generate energy, but in this coordinate
system that generation provides no net contribution because it does
not alter the specific entropy.
Thermodynamic consistency imposes the condition that heat trav-
els from hot regions to cool ones. Assuming that T increases towards
the core of the planet, this means that L(m) ≥ 0 everywhere. In a
convective atmosphere, the thermal gradient is almost independent
of the luminosity. This follows because the luminosity is determined
by the superadiabaticity of the thermal gradient, rather than by the
gradient. When convection is efficient, the convective zone is nearly
isentropic so L ∝ (∇ − ∇a)3/2 (Kippenhahn & Weigert 1990) and
the atmosphere achieves significant scaling of luminosity with only
small changes to ∇. As a result, the conditions on L cannot generally
be satisfied. This means that radiative zones are generically needed
as interfaces between convective regions. More formally, we work
in the limit of perfectly efficient convection, such that
T (P )P−∇a = constant (63)
We also make the assumption that the convective turnover time for
any region of interest is much shorter than the time-scale over which
thermal quantities change, such that convection may be assumed to
enforce an instantaneous adiabatic law.
Now suppose that we perturb a planet by injecting luminos-
ity L somewhere below the radiative–convective boundary. For
L 	 Li, we may solve equation (60) by simply reducing the lumi-
nosity escaping from deeper regions of the planet. That is, Li goes
unchanged but the luminosity in regions deeper than the injection
depth is reduced by L. In the limit of very efficient convection (or
large opacity), this adjustment holds until L ≈ Li. For L > Li
the adjustment still occurs, with the deep luminosity falling to the
radiative luminosity at the adiabatic gradient, the minimum needed
to maintain convection. The difference is that in this case there is
an excess of luminosity reaching the convective–radiative transition
and this must be accounted for. At the boundary we must have
∇a = ∇r = 3κPtLi64πGMσT 4t
, (64)
which must remain satisfied when we perturb Li so
 lnPt − 4 ln Tt +  ln κ +  lnLi = 0. (65)
If the transition temperature is similar to the photospheric temper-
ature and if the radius does not change substantially owing to the
perturbation ln L = ln (Li + Le) = 4ln Tt. Because Fe 
 Fi
and Le is fixed,
4 ln Tt ≈ Li
Le
	 1. (66)
So we may neglect the change in Tt and find
 lnPt +  ln κ +  lnLi = 0. (67)
We generally expect that, at a fixed temperature, κ increases as P
increases. As a result, Pt must fall to satisfy this relation, so either
the entropy of the central adiabat must rise or the adiabatic law must
be broken somewhere in the planet. The central entropy cannot rise
unless either heat is being added at the core or the photosphere is
hotter than the core, because heat cannot be forced to move up the
temperature gradient. Neither of these are generally the case so the
adiabatic law must be broken. As a result the planet must form an
interior radiative zone.
Figure 5. Perturbed (red) and unperturbed (black) pressure–temperature
profiles.
To characterize these radiative zones, let P1 be the transition
pressure between the surface radiative zone and the new convection
zone, P2 the transition pressure between this zone and the interior
radiative zone and P3 the transition pressure between this zone
and the central adiabat. The perturbed and unperturbed pressure–
temperature structures are shown in Fig. 5. Let Tj, mj, κ j and Lj
be the corresponding temperature, mass coordinate, opacity and
luminosity at each transition. The new convection zone is adiabatic,
so
T
1/∇a
1
P1
= T
1/∇a
2
P2
. (68)
Assuming that m ≈ M, the condition (64) for transition between
radiative and convective zones gives
P1κ1L1T
−4
1 = P2κ2L2T −42 = P3κ3L3T −43 . (69)
Finally, recalling that the central adiabat is fixed and that the new
adiabat contains a point at the fixed temperature Tt, we find
T
1/∇a
2
P2
= T
1/∇a
3
P3
(
Pt,i
Pt,f
)
, (70)
where the subscripts i and f refer to the initial unperturbed and final
perturbed system, respectively. Equation (70) thus expresses the
entropy difference between the two adiabats. Note that with these
subscript definitions,
Li,f ≡ L1. (71)
Again with κ ∝ TaPb we may write equations (65)–(71) as a system
of linear equations in the logarithms of temperature and pressure.
Solving this system yields
ln
T1
T2
= ∇a
w
ln
L1
L2
, (72)
ln
T2
T3
= ∇a
w
ln
Li,fL2
Li,iL3
, (73)
ln
P1
P2
= 1
w
ln
L1
L2
(74)
MNRAS 469, 1768–1782 (2017)
1774 A. S. Jermyn, C. A. Tout and G. I. Ogilvie
and
ln
P2
P3
= 1
w
ln
Li,fL2
Li,iL3
+ 1
1 + b ln
Li,f
Li,i
, (75)
where
w ≡ (4 − a)∇a − (1 + b). (76)
The transition temperature T1 ≡ Tt is known from the unperturbed
state so with equations (72) and (73) we may determine the remain-
ing temperatures. Likewise the unperturbed transition pressure Pt,i is
known from the unperturbed state. The perturbed transition pressure
P1 is related to the unperturbed by equation (67) so, with equations
(74) and (75) we may determine the remaining pressures.
In equilibrium, the luminosities are related by
L1 = L2 +
∫ m1
m2
(m) dm (77)
and
L2 = L3 +
∫ m2
m3
(m) dm. (78)
With these we can compute the luminosity ratios. A consequence
of equation (74) is that the new convective zone is maintained by
heat generation in between P1 and P2, or equivalently between m1
and m2, because this is what allows for L1 = L2.
The minimum luminosity required for convection may be calcu-
lated from equation (64) as
Lmin = Li ∇a∇r . (79)
Both sides of this equation are functions of pressure. We generally
expect that ∇ r rises quickly towards the interior of the planet as
convection becomes more efficient so Lmin is a small fraction of Li.
This is actually guaranteed by equation (20) so we expect that Lmin
is suppressed relative to Li by a power law in P and may calculate
L3 =
∫ m3
0
(m) dm + Li ∇a∇r(Pinject) , (80)
where Pinject is the pressure inside which minimal luminosity is
injected.
5 EX PA N S I O N
The expansion associated with changing the temperature profile of
the planet is given by
V =
∫ M
0

(
ρ−1
)
dm. (81)
In the limit where R/R is small, P/P is small at fixed m, so

(
ρ−1
) ≈ ∂T
∂ρ
∣∣∣∣
P

(
T −1
) ≈ ρ−1 T
T
≈ ρ−1 ln T . (82)
Substituting this into equation (81) we find
V ≈
∫ M
0
ρ−1 ln T dm ≈
∫ R
0
4πr2 ln T dr, (83)
where the coordinate r refers to the unheated system. The integration
proceeds up to R as an approximation, once more in the limit where
R/R is small. When this is the case and when the majority of the
heating occurs near the surface at r ≈ R this may be approximated
by
R ≈
∫ R
R−δR
 ln T dr. (84)
Now we may approximate ln T as ln (T2/T1). The pressure depth
over which this approximation (rather than ln T ≈ 0) is valid
is ln P ≈ ln (P1/Pt,i). This corresponds to a physical depth of
h ln (P1/Pt,i), because h is the characteristic scale of the thermal
properties of the planet and hence sets the scale of the radiative
zone that forms. So we may approximate equation (84) in terms of
the heating parameters by
R ≈ h ln P1
Pt,i
ln
T2
T1
. (85)
With equation (67) we find
R = − h∇a(1 + b) (w) ln
Li,f
Li,i
ln
L1
L2
. (86)
For very deep zones, the relevant scale height is that near the base of
the zone rather than the top, because the majority of the contribution
to the integral comes from this region. This may be taken into
account by noting that the scale height at the base of the radiative
zone is given by
h = − dr
d lnP
= kBT3
μmpg
. (87)
Inserting equations (72) and (73) we have
h = kBT1
μmpg
(
L2i,f
Li,iL3
)− ∇aw
. (88)
Making use of T1 ≈ Tph, equations (11) and (28) give us
h = kB
μmpg
(
L
4πσa2orbit
)1/4 ( L2i,f
Li,iL3
)− ∇aw
(89)
≈ 0.2RJ
(
L
L
)1/4 (
M
M
)−1/6 ( τorbit
10 d
)−4/3
×
(
M
MJ
)−1 (
R
RJ
)2 ( L2i,f
Li,iL3
)−∇a/w
, (90)
where τ is the orbital period. The expansion is therefore
R ≈ 0.2RJ∇a(1 + b)w
(
L
L
)1/4 (
M
M
)−1/6 ( τorbit
10 d
)−4/3 ( M
MJ
)−1
×
(
R
RJ
)2( L2i,f
Li,iL3
)−∇a/w
ln
Li,f
Li,i
ln
L2
L1
. (91)
A factor of a few from the luminosity term is therefore sufficient to
substantially inflate the planet at short orbital periods.
6 G MO D ES
The existence of internal radiative zones raises the possibility that
g modes may contribute to tidal heating. This is particularly inter-
esting because, if g-mode dissipation is the dominant form of tidal
heating,  is actually a function of the thermal structure of the planet.
That is because g modes predominantly resonate in radiative zones.
What this amounts to is a form of feedback between the thermal
and mechanical structures of the planet.
6.1 Dynamical tide
In principle there are two sources of dynamical tides, namely grav-
itational and thermal. We expect that thermal tides do not couple to
MNRAS 469, 1768–1782 (2017)
Tidal heating of hot Jupiters 1775
the g modes considered here. There are two reasons for this. First,
the thermal tide is significant only in the upper layers of the atmo-
sphere where insolation is significant. In particular, the tide damps
as e−κ (Arras & Socrates 2010). The internal radiative zone be-
gins at a comparable column density to the unperturbed radiative–
convective transition. Equation (32) gives κ0t ≈ 5 × 102, so the
damping is of the order of exp (−5 × 102), which suffices to make
this effect negligible. Secondly, the thermal tide relies on time-
scale for redistributing heat being large relative to the orbital time.
We have shown that the temperature distribution becomes spherical
very near the photosphere and well above the convection zone, even
for a tidally locked planet. This means that it will not reach even
the upper convection zone. As a result we restrict our analysis to
gravitational tides.
Due to their frequencies being small relative to the acoustic fre-
quency, g modes are unlikely to substantially compress material in
the planet. As a result we must treat them in the incompressible limit.
This may be done by separating the perturbing tidal potential into
a hydrostatic equilibrium tide and a dynamical tide (Zahn 1975).
The associated radial displacements ξ eq and ξ dyn obey the relations
(Goodman & Dickson 1998)
ξ eq = − δ
d/dr
(92)
and
∂2
∂r2
(r2ξ dyn) + ∂
∂r
(
d ln ρ
dr
r2ξ dyn
)
+ l(l + 1)
(
N2
ω2
− 1
)
ξ dyn
= l(l + 1)ξ eq − ∂
2
∂r2
(
r2ξ eq
)
, (93)
where l is the latitudinal quantum number, ω is the frequency,
 is the unperturbed planetary gravitational potential, δ is the
perturbing tidal potential due to the star and N is the Brunt–Va¨isa¨la¨
frequency, with N2 positive in the radiative zone and negative in the
surrounding convective regions. To analyse this equation we first
solve the homogeneous version, with the right-hand side set to zero,
and then compute the overlap between the resulting modes and the
forcing term given by the right-hand side.
6.2 Mode profile
The homogeneous part of equation (93) is
∂2
∂r2
(r2ξ dyn) + ∂
∂r
(
d ln ρ
dr
r2ξ dyn
)
+ l(l + 1)
(
N2
ω2
− 1
)
ξ dyn = 0. (94)
Defining
ξ ≡ r2ξ dyn, (95)
we find
∂2ξ
∂r2
+ ∂
∂r
(
d ln ρ
dr
ξ
)
+ l(l + 1)
r2
(
N2
ω2
− 1
)
ξ = 0. (96)
We now wish to perform a change of variables that will eliminate
the first-order derivative of ξ . To do this, we note that
∂2
∂r2
=
(
∂y
∂r
)2
∂2
∂y2
+ ∂
2y
∂r2
∂
∂y
. (97)
This may be written as
∂2
∂r2
=
(
∂y
∂r
)2
∂2
∂y2
+ ∂r
∂y
∂2y
∂r2
∂
∂r
. (98)
Using this, we pick
y =
∫
ρ−1 dr, (99)
which gives
∂2
∂r2
= ρ−2 ∂
2
∂y2
− d ln ρ
dr
∂
∂r
. (100)
With this substitution, equation (96) becomes
ρ−2
∂2ξ
∂y2
+ d
2 ln ρ
dr2
ξ + l(l + 1)
r2
(
N2
ω2
− 1
)
ξ = 0. (101)
Qualitatively we expect N to peak near the centre of the radiative
zone and fall to zero at the edges. To fit this, we pick a quadratic
form in our new coordinate y, such that
N2 = N20
(
1 −
(
y − y0
δy
)2)
, (102)
where r0 is the radial coordinate of the centre of the radiative zone
and 2δy is the width of the zone in y. Defining
 ≡ N0
ω
(103)
and
x ≡ y − y0
δy
, (104)
the differential equation (101) becomes
1
ρ2δy2
∂2ξ
∂x2
+ d
2 ln ρ
dr2
ξ + l(l + 1)
r2
(
2−x22−1) ξ = 0. (105)
We now define
q ≡ 1 − d
2 ln ρ
dr2
r2
l(l + 1) , (106)
such that
ρ2
δy2
∂2ξ
∂x2
+ l(l + 1)
r2
(
2 − x22 − q) ξ = 0. (107)
Note that q is positive and large because
− r2 d
2 ln ρ
dr2
≈ r
2
h2

 1. (108)
This follows because ρ has characteristic scale h and because h 	
r except near the core of the planet.
It is now worth noting that the physical width of the zone in r is
lr ≈ ρδy. (109)
This holds because for a thin zone, ρ does not change too much
across it. In thick zones there would be deviations from this which
we neglect. For convenience we now define
β ≡ l(l + 1)
(
lr
r
)2
. (110)
With this, equation (107) becomes
∂2ξ
∂x2
+ β (2 − x22 − q) ξ = 0. (111)
This may also be written as
(
2 − q) ξ =
(
2x2 − β−1 ∂
2
∂x2
)
ξ, (112)
MNRAS 469, 1768–1782 (2017)
1776 A. S. Jermyn, C. A. Tout and G. I. Ogilvie
which is the same as the equation for a quantum harmonic oscillator
with energy 2 − q, mass 2β/2 and zero-point energy /√β. The
eigenvalues are therefore quantized in the form
2 − q = 2√
β
(
1
2
+ n
)
, n = 0, 1, 2, . . . . (113)
The sign of  does not enter into equation (101) so we may take
whichever branch of the solutions to this equation that we choose.
Taking the positive we see that
n = 1 + 2n +
√
1 + 4(βq + n + n2)
2
√
β
, n = 0, 1, 2, . . . . (114)
These correspond to periods and frequencies of
Tn = 2π 1 + 2n +
√
1 + 4(βq + n + n2)
2N0
√
β
, n = 0, 1, 2, . . . (115)
and
ωn = 2N0
√
β
1 + 2n +
√
1 + 4(βq + n + n2) , n = 0, 1, 2, . . . . (116)
The radial profiles of the solutions are the product of an exponential
with a Hermite polynomial. For l = 2, the dominant tidal mode,
these are
ψn,m(r, θ, φ) ≈
√
2n
√
β√
2nn!r3
√
π
e−2n
√
βx2/2
×Hn
(
x
√
2n
√
β
)
Y2m(θ, φ), (117)
where Ylm are the spherical harmonics, n ∈ {0, 1, 2, . . . } and
m ∈{−2, 1, 0, 1, 2}. The modes are normalized so that∫
all space
d3r|ψn,m|2 = 1, (118)
and we take ρ and r as constants throughout the radiative zone to
compute this normalization. This is consistent with approximations
we make elsewhere.
6.3 Overlap integral
In order to compute the tidal forcing Fn,m(ω), we must say some-
thing about the origin of the tidal potential. There are two potential
sources, rotational asynchronization and orbital eccentricity. In the
former, the tidal forcing occurs at a frequency ωrotation − orbit,
while in the latter it occurs at a frequency of orbit. In both cases,
working in the frame corotating with the planet’s orbit,
δ ∝ GMr
2
a3orbit
, (119)
and further involves a sum of l = 2 spherical harmonics. Beyond this
the two cases differ significantly because the eccentricity case has
δ∝ e while the asynchronous case has no such factor. To capture
both cases, we write
δ = GMr
2
a3orbit
∑
m′
Y2m′ (θ, φ)km′ cos(ωt − φm′ ), (120)
where φm′ are phase factors, the factors km′ capture the magnitudes
of the various harmonics and sum in quadrature to unity and  is
a dimensionless factor of order unity in the asynchronous case and
of order e in the eccentric case. From this form and equation (92)
we may write the equilibrium tide as
ξ eq =  Mr
4
ma3orbit
∑
m′
Y2m′ (θ, φ)km′ cos(ωt − φm′ ). (121)
The driving term associated with this equilibrium tide is the right-
hand side of equation (93), given by
d(ξ ) = l(l + 1)ξ eq − ∂
2
∂r2
(r2ξ eq) = −24ξ eq. (122)
In computing the overlap of this with the eigenmodes of the homo-
geneous equation, we may treat factors of r as constant, because the
radiative zone ought to be thin on the scale of the planetary radius.
As a result, the projection is
〈ψ2,m′ |d〉 = −24
∫
ψ2,m′ (r)ξ eq(r) d3r (123)
≈ −24 Mr
6
ma3orbit
km′
∫ lr
−lr
ψ2,m′ (r ′) dr ′ (124)
= −24Mr
6lr
ma3orbit
km′
∫ ∞
−∞
ψ2,m′ (x) dx (125)
= −24 Mr
9/2lr sn
ma3orbit
√
2nn!
√
π
km′
∫ ∞
−∞
e−s
2
nx
2/2Hn(snx) dx (126)
= −24 Mr
9/2lr
ma3orbit
√
2nn!
√
π
km′
∫ ∞
−∞
e−w
2/2Hn(w) dw, (127)
where we have centred the integral on r0, the radial coordinate
corresponding to y0, and defined
sn ≡
√
2n
√
β. (128)
We have also extended the integration bounds to infinity to make the
computation easier because the exponential suppression in x makes
the precise bounds irrelevant. Note that changing variables from
r′ to x is formally quite complicated, though in the approximation
where ρ changes little over the course of the zone it just produces a
pre-factor of lr.
The integral may now be evaluated with the generating function
of the Hermite polynomials,
e2wt−t
2 =
∞∑
n=0
Hn(w) t
n
n!
, (129)
so∫ ∞
−∞
e−w
2/2Hn (w) dw = d
n
dtn
∫ ∞
−∞
e2wt−t
2−w2/2 dw
∣∣∣
t=0
(130)
=
√
2π
dn
dtn
(
et
2
)∣∣∣
t=0
(131)
= n!
√
2π
2
(
1 + n2
) (1 + (−1)n) . (132)
As a result,
〈ψ2,m′ |d〉 ≈ −24 Mr
9/2lrπ
1/4
ma3orbit
(
1 + n2
) km′
√
n!
2n−1
(133)
MNRAS 469, 1768–1782 (2017)
Tidal heating of hot Jupiters 1777
for even n and vanishes for odd n. There is complete degeneracy
in both the dissipation and oscillation over m′, so we may form a
linear combination of spherical harmonics that precisely matches
the forcing term. This amounts to summing the right-hand side of
equation (133) in quadrature over m′ and taking the square root,
which gives
〈ψ2|d〉 ≈ 24 Mr
9/2lrπ
1/4
ma3orbit
(
1 + n2
)
√
n!
2n−1
. (134)
This expression gives the amplitude of the resonance. From this
stage we take it as given that l = 2 and drop the label on ψ .
6.4 Dissipation
The square of the displacement, which is proportional to the dissi-
pation, has maxima at a distance of order ±a from the centre of the
radiative zone so, even if the dampening were uniform, we would
expect the dissipation to be greatest near the edges of the zone. In
practice, convective turbulence increases the dissipation just outside
the zone and this assertion is even stronger. To evaluate the strength
of this effect we turn to various linear dissipation mechanisms. Both
radiative and viscous damping are potentially relevant. For each of
these we may calculate a quality factor Q, giving the number of
undriven cycles required for an e-fold reduction in strength. These
combine as
Q = 11
Qrad
+ 1
Qturb
. (135)
We begin with radiative damping at finite opacity. The quality
factor of mode n is of order
Qn ≈ ωnτn ≈ 34π
(
ωnλn
c
)(
P
aT 4
)
(κρλn), (136)
where c is the speed of light and τ n, ωn and λn are the lifetime,
frequency and wavelength corresponding to mode n (Press 1986).
The wavelength is given by
λn ≈ 2lr
n + 1 , (137)
which just comes from the fact that mode n has n + 1 nodes over
the zone width of 2lr. Thus
Qn ≈ 34π
(
P
aT 4
)
4ωnκρl2r
c(n + 1)2 . (138)
Let
mz ≡ 2lrρ, (139)
an approximate zone mass. We find
Qn ≈ 34π
(
P
aT 4
)
2ωnmzκlr
c(n + 1)2 . (140)
From equation (87) we know that the scale height h is proportional
to T, so(
h
r
)4 (
P
aT 4
)
≈ Pk
4
B
m4pg
4r4a
. (141)
We are interested in regions that are sufficiently shallow so that g is
nearly constant and so
P ≈ g(M − m)
4πr2
(142)
and(
h
r
)4 (
P
aT 4
)
≈
(
k4B
am4p
)(
M − m
4πr6g3
)
(143)
=
(
k4B
4πaG3m4p
)
M − m
m3
≈ 4.6 × 105M2J (M − m)m−3, (144)
so that
Qn ≈ 1.1 × 105
(
1 − m
M
)
M−2
(
M
m
)3 2ωnmzκlr
c(n + 1)2
( r
h
)4
.
(145)
The radiative zone is stably stratified so we expect turbulent
damping to be limited to the evanescent part of the mode that
leaks into the neighbouring convective zones. The Navier–Stokes
equation with a simple viscosity term is
∂v
∂t
+ v · ∇v = g − ∇p
ρ
+ ν∇2v. (146)
We now neglect the non-linear term because we are interested in
understanding the linear growth and decay of modes and expect
the absolute velocities involved to be small. Without the non-linear
term
∂v
∂t
= g − ∇p
ρ
+ ν∇2v. (147)
The balance between gravity and the pressure gradient is what gives
us g modes, so we may write the balance as
∂v
∂t
= ωnv + ν∇2v. (148)
Now let
v′ ≡ eiωntv, (149)
so that
∂v′
∂t
= ν∇2v′. (150)
The kinetic energy density
K = 1
2
v∗ · v = 1
2
v′∗ · v′, (151)
where v∗ is the complex conjugate of v, and evolves as
∂tK = 12νv
′∗ · ∇2v′ + 1
2
νv′ · ∇†2v′∗, (152)
where ∇† is the adjunct gradient operator. When the spatial deriva-
tives are greatest in the radial direction the Laplacian just produces
a factor of (2π/λn)2 so
∂tK =
(
2π
λn
)2
νK. (153)
Integrating this equation over the whole planet we find
d
dt
∫ M
0
K dm =
∫ M
0
∂K
∂t
dm (154)
=
∫ M
0
(
2π
λn
)2
νK dm (155)
=
(
2π
λn
)2 ∫ M
0
νK dm. (156)
MNRAS 469, 1768–1782 (2017)
1778 A. S. Jermyn, C. A. Tout and G. I. Ogilvie
The viscosity ν is only significant in the convection zones on either
side of the radiative zone, where turbulent viscosity dominates. The
kinetic energy density K is only significant inside the radiative zone
and within a few wavelengths of the zone edges on either side. The
damping integral is dominated by the region in which neither is
small, so
d
dt
∫ M
0
K dm ≈
(
2π
λn
)2
λn
lr + λn ν
∫
RadiativeZone
K dm, (157)
where ν is evaluated in the convecting regions and 2(lr + λn) is
roughly the radial extent of the region where the kinetic energy
density is significant. Similarly∫
RadiativeZone
K dm ≈ mzK, (158)
where K on the right-hand side is the average kinetic energy density
in the zone. Then
d
dt
(mzK) ≈
(
2π
λn
)2
λn
lr
νmzK. (159)
Because mz is constant we find that
d lnK
dt
≈
(
2π
λn
)2
λn
lr
ν, (160)
and the damping time-scale is
τturb = λnlr4π2ν (161)
with related quality factor
Qn ≈ ωnτturb = ωnλnlr4π2ν . (162)
The turbulent diffusivity is of order vch when the convective
turnover is on a time-scale shorter than the forcing frequency ω. It is
ω rather than ωn that matters here because the oscillation physically
takes place at the driving frequency, not at the mode period. The
relevant turbulent frequency for motion over length-scale lt is
ωturb(lt) = vc(lt)
lt
. (163)
Taking vc to be given by a Kolmogorov spectrum, we find
ωturb(lt) = vc(h)
h
(
lt
h
)−2/3
. (164)
It follows that the relevant diffusive motions are on a scale
lt
h
= min
[
1,
(
vc(h)
hω
)3/2]
, (165)
and corresponding diffusivity is
ν = vchmin
[
1,
(
vc(h)
hω
)2]
. (166)
This is just the result of Goldreich & Keeley (1977) and yields a
quality factor
Qn ≈ ωnλnlr4π2vch max
[
1,
(
vc(h)
hω
)−2]
. (167)
High-frequency driving leads to a high quality factor and Qn scales
as ω2n because λn ∝ωn. This is a weaker scaling than the radiative
Q, which goes as ω3n ∝ (n + 1)−3. Thus at large n the convective
mechanism dominates.
We are often interested in the lowest n because this mode is the
least suppressed by overlap factors. The convective flux of interest
is generally Fi, which equation (31) shows is of the order of 10−5Fe.
The external flux is typically about 10−2 F, so the relevant con-
vective flux is of the order of 5 × 103 erg cm−2 s−1. For a density of
10−1 g cm−3 the convection speed is then
vc ≈
(
Fi
ρ
)1/3
≈ 30 cm s−1. (168)
So for a scale height of 109 cm, vc/h ≈ 3 × 10−8 Hz. We show later
that we are interested in frequencies on the order of 10−6 Hz. This
means that the factor [hω/vc(h)]2 accounting for the eddy time is of
order 105, so the convective quality factor is Q1 ≈ 300. By compar-
ison, the fiducial radiative quality factor with the same assumptions
is Q1 ≈ 1. Thus we expect radiative damping to dominate by a
reasonable margin unless the fluxes involved are many orders of
magnitude larger or the frequencies are several orders of magnitude
smaller.
6.5 Boundaries
We have shown that, when g-modes dominate the dissipation, ε is
significant primarily near the edges of the radiative zone. It is also
straightforward to show that the dissipation is symmetric because
the squared mode profiles are even. So ε is an even function, the
integral of which is dominated by the regions just outside the zone
boundaries. Suppose that the total luminosity produced by tides is
Lt and that a fraction f of this is produced inside the radiative zone.
In steady-state
L2 − L3 = fLt (169)
and
L1 − L2 = 12 (1 − f )Lt. (170)
Using equation (80) we find
L3 = Li ∇a∇r(Pinject) +
1
2
(1 − f )Lt. (171)
If the heating is large relative to Li and f is small then we expect
L3 ≈ 12 (1 − f )Lt, (172)
so
L1
L2
≈ 2 (173)
and
L2
L3
≈ 1. (174)
Using equations (72)–(75) we find
ln
T1
T2
= ∇a
w
ln 2, (175)
ln
T2
T3
= ∇a
w
ln
Lt
Li,i
, (176)
ln
P1
P2
= 1
w
ln 2 (177)
and
ln
P2
P3
=
(
1
w
+ 1
1 + b
)
ln
Lt
Li,i
. (178)
MNRAS 469, 1768–1782 (2017)
Tidal heating of hot Jupiters 1779
With equations (172) and (173), equation (91) yields
R
RJ
≈ −0.1 ∇a(1 + b) (w)
(
L
L
)1/4 (
M
M
)−1/6 ( τorbit
10 d
)−4/3
×
(
M
MJ
)−1 (
R
RJ
)2 ( 2Li,f
Li,i
)−∇a/w
ln
Li,f
Li,i
. (179)
Recall from equation (116) that the resonant frequencies are
ωn = 2N0
√
β
1 + 2n +
√
1 + 4(βq + n + n2) (180)
= 2N0
√
β
(1 + 2n)
(
1 +
√
1 + 4βq(1+2n)2
) (181)
= 2N0
√
β
(1 + 2n)
(
1 +
√
1 + 4(lr /h)2(1+2n)2
). (182)
The right-hand side has two characteristic regimes, one in which
lr/h is large and one in which it is small or of order unity. In the
former case,ωn ∝ h, while in the latterωn ∝ lr. Both h and lr increase
with P3/P2, so the resonant frequency goes up as the zone width
increases. This means that the resonance shifts up in frequency as
the luminosity increases. That is,
dω0
dLt
> 0. (183)
So an increase in luminosity tends to tune the system towards res-
onance if it is being driven above resonance and pushes it away
from resonance otherwise. The net result is that there is thermo-
mechanical feedback that tends to bias systems towards resonance,
particularly when their resonant frequency is below the driving fre-
quency.
6.6 Power production
Each mode may be treated as a separate damped and forced har-
monic oscillator. Let ξ n be the amplitude for mode n. Then
ω2nξn +
ωn
Qn
˙ξn + ¨ξn = ω2〈ψ |d〉eiωt . (184)
We may solve this differential equation in steady-state and fix the
reference phase to find
ξn = ω
2〈ψ |d〉eiωt
ω2n − ω2 + iωnωQ−1n
. (185)
The power dissipated is
Pn = ρ
(
˙ξ ∗n e
iωtω2〈ψ |d〉) (186)
= ρ (−iξ ∗n eiωt)ω3〈ψ |d〉 (187)
= ρ (ξ ∗n eiωt)ω3〈ψ |d〉 (188)
= ρ
(
1
ω2n − ω2 − iωnωQ−1n
)
ω5|〈ψ |d〉|2 (189)
= ωnωQ
−1
n
(ω2n − ω2)2 + (ωnωQ−1n )2
ρω5|〈ψ |d〉|2. (190)
With equation (134) this becomes
Pn = ωnω
6Q−1n(
ω2n − ω2
)2 + (ωnωQ−1n )2 q (191)
= ω
3
Qnω−1n ω
(
ω2n
ω2
− 1
)2
+ ω−1ωnQ−1n
q, (192)
where
q ≡ 5762ρM
2
 r
9l2rπ
1/2
m2a6orbit
(
n!
2n−1
(
1 + n2
)2
)
. (193)
If the tides are driven by the rotational energy of the planet, then
ω is the planet’s rotation frequency. In many cases however the tides
are driven by either orbital eccentricity, in which case the forcing
frequency is just the orbital frequency (Arras & Socrates 2010). In
this more generally interesting case, the driving frequency is
ω = orbit = 2π
τorbit
≈ 7 × 10−6 Hz
( τorbit
10 d
)−1
. (194)
To compare, the highest resonant frequency occurs when n = 0 and
is
ω0 = 2N0
√
β
1 + √1 + 4βq , (195)
equation (116). Because βq ≈ 1,
ω0 ≈ N0
√
β
1 + √5 ≈
N0lr
R
(196)
for l = 2. Now N0 is the peak Brunt–Va¨isa¨la¨ frequency in the
radiative zone where the temperature gradient is substantially sub-
adiabatic, so
N20 ≈ ∇a
g
h
. (197)
This gives
ω0 ≈ lr
√∇a
R
√
g
h
. (198)
Now lr ≈ h ln (P1/Pt, i) so
ω0 ≈
√
∇a ln P1
Pt,i
√
gh
R2
(199)
= 1√
γ
( cs
R
)
ln
P1
Pt,i
(200)
≈ 2 × 10−5 Hz
(
T
1 × 103 K
)1/2 (
R
RJ
)−1
ln
P1
Pt,i
. (201)
This is quite close to the orbital frequency so we expect that res-
onances are not uncommon. Also note that it is somewhat greater
than the orbital frequency, so there will generally be modes with
frequencies lower than the orbital frequency that are pulled upward
towards it by thermomechanical feedback.
The precise shape and spacing of the resonances depend on our
ansatz for the Brunt–Va¨isa¨la¨ frequency in the radiative zone and so
it is not useful to make predictions that depend on our chosen form.
If we instead average over resonances, we note that the forcing
integral falls exponentially in n while the resonances fall off as
a power law in n so that the n = 0 resonance always dominates
on average. Individual systems far from the n = 0 resonance may
MNRAS 469, 1768–1782 (2017)
1780 A. S. Jermyn, C. A. Tout and G. I. Ogilvie
exhibit significant dissipation by virtue of sitting directly on a higher
resonance but we expect this to be rare. So we assume that ω and
ωn are of the same order and write the net power
P =
∑
i
Pi (202)
≈ P0 (203)
≈ 2 × 103Q−10 2ρ
3orbitM
2
 r
9l2r
m2a6orbit
(
ω0
orbit
)
. (204)
This may be converted to a flux as
Ft = P4πr2 ≈ 2 × 10
2Q−10 
2ρ
3orbitM
2
 r
7l2r
m2a6orbit
(
ω0
orbit
)
. (205)
7 EQU I L I B R I U M R A D I U S
To first order, suppose that lr ≈ h. The tidal flux is then given by
Ft ≈ 2 × 102Q−10 2ρ
3orbitM
2
 r
7h2
m2a6orbit
(
ω0
orbit
)
. (206)
Inserting equation (145) yields
Ft ≈ 2 × 10−32ρ
3
orbitmM2r3h5M2 c
2ω0mzκM3a6orbit
(
ω0
orbit
)(
1 − m
M
)−1
.
(207)
If the radiative zone is shallow but dominates the mass above its
base, then m ≈ M. In addition,
M − m ≈ 4πr2ρlr ≈ 4πr2ρh, (208)
so
Ft ≈ 2 × 10−32ρ r
3h52orbitMπ
1/2c
8πr2ρ2h2κa6orbit
(
M
MJ
)2
(209)
≈ 2 × 1032ρ rh
52orbitMc
8π(ρh)2κa6orbit
(
M
M
)2
. (210)
Noting that
t ≈ hρ, (211)
we find
Ft ≈ 8 × 1012 rh
42orbitMc
tκa
6
orbit
(
M
M
)2
. (212)
Recalling equation (32) and neglecting the logarithmic correction
to t owing to the motion of the zone boundary we write
Ft ≈ 5 × 10−5F2
( τorbit
10 d
)−31/6 ( M
M
)5/12 (
R0
RJ
)
×R−2M8/3
(
L
L
)−5/8
ψ5/2
(
κ
κ0
)−1 (
h
r
)4
. (213)
With the fiducial values and h ≈ 0.2r this flux produces an expansion
at a rate
dR
dt
≈ Ft
ρcpT
≈ F
Pt
≈ 3 × 10−5 cm s−1 2, (214)
which is sufficient to produce expansion of order RJ over Myr time-
scales.
As discussed in Section 5 the expansion eventually increases the
escaping flux to match the generated flux, so the expansion does
not continue forever. The relevant dimensionless parameter for this
equilibrium is the ratio Ft to Fi, i, which is the unperturbed Fi. Recall
from equation (31) that
Fi ≈ 5 × 10−4FeM5/3
(
L
L
)−5/8 (
aorbit
R
)5/4 (
κ
κ0
)−1
ψ5/2
R3 .
(215)
Inserting equation (30) yields
Fi ≈ 1.6 × 10−4 FM5/3
(
L
L
)3/8 (
aorbit
R
)−3/4 (
κ
κ0
)−1
ψ5/2
R3
(216)
≈ 3 × 10−6 FM5/3
(
L
L
)3/8 ( τorbit
10 d
)−1/2
×
(
M
M
)−1/4 (
κ
κ0
)−1
ψ5/2
R3 . (217)
Thus
2Ft
Fi,i
≈ 3 × 1072MR−2
(
L
L
)−1 (
R0
RJ
)
×
( τorbit
10 d
)−14/3 ( M
M
)2/3 (
h
R
)4
. (218)
Inserting equation (90) we get
2Ft
Fi,i
≈ 5 × 1042M−3R2
(
R0
RJ
)5 ( τorbit
10 d
)−10 ( L2i,f
Li,iL3
)−7∇a/w
.
(219)
If the expansion is small, then the luminosity ratios may be re-
placed by flux ratios. If tidal heating is significant, the perturbed
flux escaping from the interior of the planet Fi, f ≈ Ft, so
2Ft
Fi,i
≈ 5 × 1042M−3R2
(
R0
RJ
)5 ( τorbit
10 d
)−10 ( F 2i,f
Fi,iF3
)−7∇a/w
.
(220)
Using equations (173) and (174) we find
2Ft
Fi,i
≈ 5 × 1042M−3R2
(
R0
RJ
)5 ( τorbit
10 d
)−10 ( 2Fi,t
Fi,i
)−7∇a/w
.
(221)
Note that R and R must be the equilibrium radii here, not pre-
expansion radii. Note that the absolute magnitude of the opacity
does not enter our final expression. This cancellation is due to the
assumption of efficient convection.
The quantity
f ≡ −7∇a
w
(222)
is of key importance to the nature of the solution. If f < 1, then the
solution is stable, meaning that a system initially perturbed away
from this equilibrium solution returns to it over time. What this
means physically is that an increase in the generated flux increases
the temperature at the base of the radiative zone enough that the
flux that escapes increases by more, leading to a negative feedback
loop. If f > 1, then the solution is unstable, meaning that an increase
MNRAS 469, 1768–1782 (2017)
Tidal heating of hot Jupiters 1781
in the generated flux increases the temperature in the radiative zone
by less than what is required to allow that additional flux to escape,
leading to a positive feedback loop. If a system has f > 1, then
either the initial perturbation has the radiative zone deep enough
that the runaway process keeps increasing the flux until some of our
assumptions break down or the initial perturbation has the radiative
zone shallow enough that the runaway process prevents it from
migrating inward causing it to stay where it initially forms. If the
pre-factor on the right-hand side of equation (221) exceeds unity,
then the unstable branch is always the relevant one because the initial
perturbation must yield a ratio of at least unity in order to cause a
radiative zone to form. If the pre-factor is less than one, then any
radiative zone formed simply stays where the initial perturbation
produces it.
Despite the uncertainties in the precise numbers involved,
our results are fairly robust because the right-hand side of
equation (221) scales as τ−10orbit , and this scaling only becomes
stronger once f is taken into account. For orbital periods shorter
than of about 10 d, with some uncertainty, the flux generated may
exceed the flux escaping from the centre of the planet. When the
tidal flux dominates the expansion is given by equation (179) and
may be approximated by
R
RJ
≈ 0.1 ∇a(1 + b) (w)
(
L
L
)1/4 (
M
M
)−1/6 ( τorbit
10 d
)−4/3
×
(
M
MJ
)−1 (
R
RJ
)2 ( 2Fi,f
Fi,i
)−∇a/w
ln
Fi,f
Fi,i
. (223)
Neglecting the logarithmic dependence, this may be combined with
equation (221) to give
R
RJ
≈ 0.1 ∇a(1 + b) (w)
(
L
L
)1/4 (
M
M
)−1/6
×
( τorbit
10 d
)−( 43 + 10f7(1−f ) ) ( M
MJ
)−1 (
R
RJ
)2
×
(
5 × 1042M−3R2
(
R0
RJ
)5) f7(1−f )
(224)
Even though the dependence on the flux ratio is small, the ratio itself
can be quite large, particularly at smaller periods. Many cases, such
as a = 0, b = 2 or a = 1, b = 1, have f > 1 and so orbital periods of
order 30 d1/5 suffice to cause expansion of order R0. It is difficult
to say more because many of our approximations break down at this
point. If f < 1, as can be achieved for example with a = 4, b = 2,
R/R ∝ τ−64/3orbit , orbital periods of order 20 d3/16 suffice to pro-
duce unit expansion. This is consistent with most of the known
cases of highly inflated Jupiter-mass planets, assuming  ≈ e ≈
0.1. For low-mass planets, the lower surface gravity makes larger
expansion more feasible. This has been recently observed (Bakos
et al. 2016). Stronger claims are difficult to make analytically be-
cause the dependence of the flux ratio on the specifics of the opacity
are quite severe and the detailed compositions of the atmospheres
of exoplanets at intermediate depths remain largely unknown. Pre-
cision studies of this thermomechanical feedback will likely require
numerical tools in all but the simplest cases.
At sufficiently low masses, large flux ratios become impossible
to attain given the factor of M−3 in equation (221). At this point
further expansion is impossible. Likewise at large enough radii the
central adiabat disappears so that much of this analysis becomes
invalid. We do not expect this to be a limiting factor, however. At
short periods Roche lobe overflow becomes a substantial barrier.
Substantial changes in the opacity may also occur, particularly if
the Kramer regime becomes relevant, and this may invalidate much
of the analysis too. In addition at large radii the neglected factors
of R in converting from luminosities to fluxes become relevant and
these act to limit the expansion.
8 ENERGETI C TI ME-SCALES
If the tides are eccentricity driven, it is important to consider the
time-scale over which the orbit circularizes. It suffices to the level
of accuracy of interest to note that the energy that may be extracted
from an orbit of eccentricity e < 1 is of order e2v2orbitM . The circu-
larization time-scale is therefore
τcirc ≈ e
2v2orbitM
P ≈
( τorbit
10 d
)5/2
2 × 1012 yr, (225)
where we have taken  ≈ e and used our fiducial values for all
parameters other than h/r, which we have taken to be 0.2. From this
it is clear that most systems of interest can be eccentricity driven for
billion-year time-scales, even if they require shorter periods than
the fiducial.
If the tides are driven by the planet’s rotation, they generally have
many orders of magnitude less energy to draw from, and so are not
sustained on the time-scales of interest. They may still produce
bloating, but not for long enough to be easily observable.
9 C O M PA R I S O N
For comparison with other mechanisms it is useful to compute the
Love number associated with these g modes. The Love number is
defined in terms of the power and perturbing tidal potential as
[kml (ω)] =
8πGP
(2l + 1)r|δ|2ω (226)
(Ogilvie 2014). Summing in quadrature over all m and using l = 2
and equation (213) we find
(k) ≈ 3 × 10−6
( τorbit
10 d
)−17/6 ( M
M
)−11/12
×R−3M8/3
(
L
L
)−5/8
ψ5/2
(
κ
κ0
)−1 (
h
0.2r
)4
. (227)
For comparison, inertial waves result in a frequency-averaged value
of
(k) ≈ 5R
32
GM
(
Rc
R
)5
(228)
≈ 8 × 10−9
(
M
MJ
)−1 (
R
RJ
)3 ( τorbit
10 d
)−2 ( Rc
0.1R
)5
(229)
(Ogilvie 2013), where Rc is the core radius. Likewise viscoelastic
dissipation in a solid core potentially yields (k) ≈ 10−6, scaling
with the elastic properties of the core (Remus et al. 2012; Guenel,
Mathis & Remus 2014) and once more averaging over frequency.
These mechanisms are therefore of comparable order, depending
on precisely which fiducial value is chosen. The reason the g-mode
mechanism inflates planets more readily despite having comparable
or somewhat less power dissipation is that it heats primarily near
the surface where the radius is more easily perturbed.
MNRAS 469, 1768–1782 (2017)
1782 A. S. Jermyn, C. A. Tout and G. I. Ogilvie
1 0 C O N C L U S I O N S
We have characterized the response of heavily insolated Jupiter-like
planets to tidal heating for a wide range of tidal heating models. A
necessary condition for significant bloating of these planets is deep
heating. We find that tidal heating, either directly through tide–core
interactions or indirectly through resonance-sensitive migration of
radiative zones induced by g modes, is of the right order of mag-
nitude to induce the observed bloating if it is sufficiently deep. We
have further shown that nearly every tidal heating model results in
deep heating so long as the atmosphere is sufficiently irradiated. This
explains the observation of substantially bloated hot Jupiters with
a physically reasonable orbital period cut-off for such effects. The
migration of interior radiative zones provides a natural explanation
for the matching of tidal frequency with orbital frequency despite
the observed wide range of orbital frequencies of hot Jupiters.
This entire analysis hinges on there being a luminosity perturba-
tion to start. This luminosity then produces a self-sustaining interior
radiative zone that dissipates substantially more heat. The initial
perturbation may come from non-linear instabilities and so may
provide an indirect probe of these effects. It may also come from
planetary migration. If a planet migrates inward and if opacity falls
as a result, the incident flux can temporarily force the creation of a
radiative zone while the convection zone adjusts to the reduced flux
it must carry. The g-mode hysteresis described in this paper may
then prevent the zone from disappearing, even if its location shifts
to better match resonance. Inflated planets may therefore carry a
record of their migration histories.
Finally, the thermomechanical feedback mechanism we propose
highlights the importance of considering dynamical effects across
many time-scales. Feedback is possible both from short time-scales
to long, as in tidal heating, and from long time-scales to short, as
in the dynamical tuning of g modes. By their very nature couplings
across so many scales are difficult to track down and so there may
be many more that have yet to be discovered.
AC K N OW L E D G E M E N T S
We thank Sterl Phinney for productive conversations on planetary at-
mospheres and heated planets. ASJ acknowledges support from the
Goldwater scholarship and the Marshall scholarship. CAT thanks
Churchill College for his fellowship.
R E F E R E N C E S
Arras P., Socrates A., 2010, ApJ, 714, 1
Bakos G. ´A. et al., 2016, AAS, preprint (arXiv:1606.04556)
Batygin K., Stevenson D. J., Bodenheimer P. H., 2011, ApJ, 738, 1
Galperin B., Sukoriansky S., Anderson P. S., 2007, Atmos. Sci. Lett., 8, 65
Ginzburg S., Sari R., 2015, ApJ, 803, 111
Goldreich P., Keeley D. A., 1977, ApJ, 211, 934
Goodman J., Dickson E. S., 1998, ApJ, 507, 938
Guenel M., Mathis S., Remus F., 2014, A&A, 566, L9
Hartman J. D. et al., 2012, AJ, 144, 139
Hartman J. D. et al., 2016, AJ, 152, 182
Hellier C. et al., 2012, MNRAS, 426, 739
Kippenhahn R., Weigert A., 1990, Stellar Structure and Evolution. Springer-
Verlag, Berlin
Lopez E. D., Fortney J. J., 2016, ApJ, 818, 4
Mathis S., Palacios A., Zahn J.-P., 2004, A&A, 425, 243
Miller N., Fortney J. J., Jackson B., 2009, ApJ, 702, 1413
Ogilvie G. I., 2013, MNRAS, 429, 613
Ogilvie G. I., 2014, ARA&A, 52, 171
Press W. H., 1986, in Sturrock P. A., Holzer T. E., Mihalas D. M., Ulrich R.
K., eds, Physics of the Sun. Reidel, Dordrecht, Vol. 1, chapt. 4, p. 77
Rein H., 2012, preprint (arXiv:1211.7121)
Remus F., Mathis S., Zahn J.-P., Lainey V., 2012, A&A, 541, A165
Socrates A., 2013, preprint (arXiv:1304.4121)
Spiegel D. S., Burrows A., 2013, ApJ, 772, 76
Stevenson D. J., 1991, ARA&A, 29, 163
Weiss L. M. et al., 2013, ApJ, 768, 14
Youdin A. N., Mitchell J. L., 2010, ApJ, 721, 1113
Zahn J.-P., 1975, A&A, 41, 329
This paper has been typeset from a TEX/LATEX file prepared by the author.
MNRAS 469, 1768–1782 (2017)