# Tidal dissipation in rotating and evolving giant planets with application to
exoplanet systems

Yaroslav A. Lazovik,1,2 Adrian J. Barker,3 Nils B. de Vries3 and Aurélie
Astoul3
1Lomonosov Moscow State University, Faculty of Physics, 1 Leninskie Gory,
bldg.2, Moscow, 119991, Russia
2Sternberg Astronomical Institute, Lomonosov Moscow State University,
Universitetsky pr. 13, Moscow, 119234, Russia
3 School of Mathematics, University of Leeds, Leeds LS2 9JT, UK E-mail:
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
(Accepted 2023 November 26. Received 2023 November 24; in original form 2023
August 07)

###### Abstract

We study tidal dissipation in models of rotating giant planets with masses in
the range $0.1-10M_{\mathrm{J}}$ throughout their evolution. Our models
incorporate a frequency-dependent turbulent effective viscosity acting on
equilibrium tides (including its modification by rapid rotation consistent
with hydrodynamical simulations) and inertial waves in convection zones, and
internal gravity waves in the thin radiative atmospheres. We consider a range
of planetary evolutionary models for various masses and strengths of stellar
instellation. Dissipation of inertial waves is computed using a frequency-
averaged formalism fully accounting for planetary structures. Dissipation of
gravity waves in the radiation zone is computed assuming these waves are
launched adiabatically and are subsequently fully damped (by wave
breaking/radiative damping). We compute modified tidal quality factors
$Q^{\prime}$ and evolutionary timescales for these planets as a function of
their ages. We find inertial waves to be the dominant mechanism of tidal
dissipation in giant planets whenever they are excited. Their excitation
requires the tidal period ($P_{\mathrm{tide}}$) to be longer than half the
planetary rotation ($P_{\mathrm{rot}}/2$), and we predict inertial waves to
provide a typical $Q^{\prime}\sim 10^{3}(P_{\mathrm{rot}}/1\mathrm{d})^{2}$,
with values between $10^{5}$ and $10^{6}$ for a 10-day period. We show
correlations of observed exoplanet eccentricities with tidal circularisation
timescale predictions, highlighting the key role of planetary tides. A major
uncertainty in planetary models is the role of stably-stratified layers
resulting from compositional gradients, which we do not account for here, but
which could modify predictions for tidal dissipation rates.

###### keywords:

planet-star interactions – planetary systems – planets and satellites:
interiors – planets and satellites: physical evolution

††pubyear: 2023††pagerange: Tidal dissipation in rotating and evolving giant
planets with application to exoplanet systems–Tidal dissipation in rotating
and evolving giant planets with application to exoplanet systems

##  1 Introduction

Tidal interactions play a major role in the dynamics of star-planet and
stellar binary systems, leading to planetary orbital migration (e.g. Jackson
et al., 2008, 2016; Villaver & Livio, 2009; Bolmont & Mathis, 2016; Gallet et
al., 2017; Bolmont et al., 2017; Ahuir et al., 2021; Lazovik, 2021), orbital
circularization (e.g. Witte & Savonije, 2002; Nagasawa et al., 2008; Beaugé &
Nesvorný, 2012; Barker, 2022), spin-orbit re-alignment (e.g. Barker & Ogilvie,
2009; Winn et al., 2010; Lai, 2012; Hamer & Schlaufman, 2022), and rotational
evolution (e.g. Bolmont et al., 2012; Penev et al., 2018; Gallet et al., 2018;
Gallet & Delorme, 2019; Barker, 2022). Indeed, tidal interactions alter the
architectures of exoplanetary systems, and the rotations of stars and planets.
However, theoretical predictions for tidal evolution are currently uncertain,
and they highly depend on the specific prescriptions employed in the
aforementioned papers. This motivates us to work towards developing more
realistic treatments of tidal flows in stars and planets.

Some early steps towards a theory of tides were made long before the discovery
of the first exoplanet from studying the Earth-Moon system (Darwin, 1879).
More than a hundred years of research has substantially refined our knowledge
of tidal interactions in fluid bodies such as stars and giant planets.
Motivated by binary stars, some crucial steps were made by Zahn (1975, 1977,
1989), who separated the contributions of equilibrium (non-wavelike) and
dynamical (wavelike) tides in linear theory. Dynamical tides can be further
divided into inertial waves (hereafter IWs, or magneto-inertial) and internal
gravity waves (hereafter GWs, or magneto-gravito-inertial), propagating in the
convective and radiative regions, respectively (e.g. Terquem et al., 1998;
Goodman & Dickson, 1998; Ogilvie & Lin, 2004; Wu, 2005b; Ogilvie & Lin, 2007;
Goodman & Lackner, 2009; Weinberg et al., 2012; Ivanov et al., 2013; Lin &
Ogilvie, 2018). Recently, Barker (2020, hereafter B20) applied the latest
tidal theory to compute modified tidal quality factors $Q^{\prime}$
(essentially the ratio of the maximum tidal energy stored to the energy
dissipated in one tidal period – a quantity essential for tidal modelling) and
tidal evolutionary timescales in a range of stellar models, developing
prescriptions for the dissipation of tides of various types (i.e. equilibrium
tides, IWs and GWs). These were implemented in Lazovik (2021, 2023) to explore
the secular evolution of hot Jupiter systems over a wide parameter space, and
employed to provide an explanation for close solar-type binary circularization
in Barker (2022).

Even after the first exoplanet detection, attention has primarily focused on
computing tidal dissipation inside stars, while dissipation in planetary
interiors has typically been restricted to Solar System objects. A
comprehensive application of tidal theory to rotating planetary models was
performed in Ogilvie & Lin (2004) (see also Wu 2005b and Ivanov & Papaloizou
2007), which opened the doors to a new direction of research. Giant planets
share many similarities in structure with stars. They are both primarily
multi-layer fluid bodies containing both convective and radiative regions,
albeit planets may also have solid cores. But there are important differences:
planets tend to rotate faster (compared with both their dynamical and
convective frequencies) such that IWs are almost always likely to be important
and hot Jupiters also have larger tidal amplitudes than planet-hosting stars,
thereby requiring consideration of nonlinear tidal mechanisms like the
elliptical instability (that excites IWs in convection zones, e.g. Barker,
2016; de Vries et al., 2023) or other nonlinear IW interactions (e.g. Favier
et al., 2014; Astoul & Barker, 2022, 2023). The interaction between
equilibrium tides and turbulent convection is also expected to be far into the
regime of fast tides (where tidal frequencies exceed convective turnover
frequencies) because convection is typically much slower in planets than in
stars (Goldreich & Nicholson, 1977; de Vries et al., 2023). Hence, the
reduction in the turbulent viscosity for fast tides must be accounted for
(e.g. Goldreich & Nicholson, 1977; Ogilvie & Lesur, 2012; Vidal & Barker,
2020a; Duguid et al., 2020b). Furthermore, convection is likely to be
influenced by rapid rotation, which can modify this interaction (e.g.
Stevenson, 1979; Barker et al., 2014; Mathis et al., 2016; Currie et al.,
2020; Dandoy et al., 2023; de Vries et al., 2023).

The role of IWs for planetary tidal dissipation has been explored in prior
work (e.g. Ogilvie & Lin, 2004; Wu, 2005a, b; Ogilvie, 2009; Ivanov &
Papaloizou, 2007; Goodman & Lackner, 2009; Ivanov & Papaloizou, 2010;
Papaloizou & Ivanov, 2010; Ogilvie, 2013; de Vries et al., 2023). However, a
detailed study of the evolution of tidal dissipation rates and tidal quality
factors $Q^{\prime}$ from equilibrium and dynamical tides following planetary
evolution has never been performed previously to our knowledge, though Terquem
& Martin (2021) performed computations for just the equilibrium tide, assuming
a different mechanism dissipates the tidal flow to what we consider. We use
new models of giant planet interiors with masses in the range from $0.1$ to
$10M_{\mathrm{J}}$ computed with the MESA code (Paxton et al., 2011, 2013,
2015, 2018, 2019), with various strengths of stellar instellation so as to
model both hot and cold planets, to theoretically calculate tidal dissipation
rates (thereby extending Barker, 2020, for low-mass stars). In Sec. 2, we
describe our model. Our results are presented in Sec. 3 and applied to
exoplanet eccentricities in Sec. 4.

##  2 Methods

###  2.1 Tidal dissipation mechanisms

We consider a giant planet of mass $M_{\mathrm{pl}}$ and radius
$R_{\mathrm{pl}}$ and employ spherical coordinates centred on the body with
radial coordinate $r$. The intensity of tidal dissipation is often quantified
by the (modified) tidal quality factor $Q^{\prime}$, which is proportional to
the ratio of the maximum energy stored in the tide to the amount dissipated in
one period (e.g. Goldreich, 1963; Ogilvie, 2014). More effective dissipation
corresponds to lower $Q^{\prime}$. Here we consider three mechanisms of tidal
dissipation: equilibrium (non-wavelike) tides damped by their interaction with
turbulent (rotating) convection, IWs, and GWs. Each mechanism is characterized
by its corresponding tidal quality factor ($Q^{\prime}_{\mathrm{eq}}$,
$Q^{\prime}_{\mathrm{iw}}$, and $Q^{\prime}_{\mathrm{gw}}$ for equilibrium
tides, IWs, and GWs, respectively), and these are calculated within the
formalism of B20 (building upon many prior works) with a few modifications
that we will describe below. We focus on tides with spherical harmonic degree
$l=2$ and azimuthal wavenumber $m=2$, which is usually the dominant component
in systems with low obliquities. This is likely to be the dominant component
of tidal forcing in asynchronously rotating bodies, as well as for
eccentricity tides in weakly eccentric but spin-synchronised bodies – see e.g.
Eqs. 4 and 5 in Ogilvie & Lin (2007) or Table 1 in Ogilvie (2014).

The equilibrium tide is a quasi-static fluid response of a perturbed body that
is thought to be dissipated through the action of “turbulent viscosity" in
convective zones. The corresponding tidal quality factor is obtained via the
expression:

$\frac{1}{Q^{\prime}_{\mathrm{eq}}}=\frac{16\pi
G}{3(2l+1)R_{\mathrm{pl}}^{2l+1}|A|^{2}}\frac{D_{\mathrm{v}}}{|\omega_{\mathrm{t}}|},$
(1)

where $G$ is the gravitational constant, $D_{\mathrm{v}}$ is the rate of
viscous dissipation of the equilibrium tide, and $A$ is the amplitude of the
tidal potential component (this will not be specified further as it cancels
because $D_{\mathrm{v}}\propto|A|^{2}$ in linear theory). The tidal forcing
frequency is $\omega_{\mathrm{t}}=2\pi/P_{\mathrm{tide}}$ (where
$P_{\mathrm{tide}}$ is the tidal period), which is
$\omega_{\mathrm{t}}=2(n-\Omega_{\mathrm{pl}})$ for a circular, aligned orbit,
where $n$ and $\Omega_{\mathrm{pl}}$ are the orbital mean motion and planetary
spin, respectively. For an eccentric orbit with synchronised111Perhaps the
spin should instead be pseudo-synchronised for an eccentric orbit, but it is
not clear that the classical formula of Hut (1981) is valid since it is
derived by assuming equilibrium tide damping with a constant time-lag (see
also Ivanov & Papaloizou, 2004). (and aligned) spin we instead have
$|\omega_{\mathrm{t}}|=n$.

The viscous dissipation rate $D_{\mathrm{v}}$ is computed using Eqs. (20) to
(22) in B20. This quantity depends on the equilibrium tidal displacement
vector (defined in section 2 of B20) and the turbulent effective viscosity
$\nu_{\mathrm{E}}$ at each radius. The latter is assumed to be a function of
radius and to act like an isotropic shear viscosity, linked to the crude
mixing-length theory expectation $\nu_{\mathrm{MLT}}\propto
u_{\mathrm{c}}l_{\mathrm{c}}$, with $u_{\mathrm{c}}$ the convective velocity
and $l_{\mathrm{c}}=\alpha H_{p}$ the mixing-length ($\alpha$ is mixing-length
parameter and $H_{p}$ is pressure scale height). As demonstrated in
hydrodynamical simulations (e.g. Ogilvie & Lesur, 2012; Duguid et al., 2020b;
Vidal & Barker, 2020a) and as previously hypothesised using phenomenological
arguments (Zahn, 1966; Goldreich & Nicholson, 1977, even if the arguments of
both are not supported by simulations, despite agreement regarding the final
result with the latter), $\nu_{\mathrm{E}}$ is reduced for fast tides (where
$\omega_{\mathrm{t}}>\omega_{c}$, with the convective frequency
$\omega_{\mathrm{c}}=u_{\mathrm{c}}/l_{\mathrm{c}}$) in a frequency-dependent
manner. This can be accounted for using a piece-wise continuous correction
factor depending on the ratio $\omega_{\mathrm{t}}/\omega_{\mathrm{c}}$ (for
which we employ Eq. (27) of B20, inferred from detailed numerical simulations)
at each radius in the planet. Moreover, the rapid rotation expected for giant
planets stabilizes convection on large length scales, and a steeper
temperature gradient is required to sustain a given heat flux (e.g. Stevenson,
1979; Barker et al., 2014; Currie et al., 2020). Mathis et al. (2016) and de
Vries et al. (2023) have shown that such rapid rotation reduces
$\nu_{\mathrm{E}}$ even further (though interestingly the regime with
$\omega_{\mathrm{t}}\gg\omega_{\mathrm{c}}$ is not modified by rotation).
Following these works, and as confirmed by their simulations, we take into
account rapid rotation according to rotating mixing-length theory (RMLT), by
multiplying $l_{\mathrm{c}}$ and $u_{\mathrm{c}}$ at each radius by $\rm
Ro^{3/5}$ and $\rm Ro^{1/5}$, respectively, where $\rm Ro$ is the convective
Rossby number ($\mathrm{Ro}=\omega_{\mathrm{c}}/\Omega_{\mathrm{pl}}$, based
on the non-rotating convective frequency). We demonstrate in Sec. 3 that the
dissipation of equilibrium tides is insufficient to provide significant
orbital or spin evolution compared with wavelike tides.

Motivated by results from the (albeit very idealised), numerical simulations
described above, we assume equilibrium tides are damped by their interaction
with convection in a way that can be modelled as a local (frequency-dependent)
effective viscosity that is positive at each radial location in the planet
(and is isotropic for simplicity). Negative values for $\nu_{E}$ –
corresponding to tidal anti-dissipation – have been found to occur,
particularly at very high tidal frequencies (Ogilvie & Lesur, 2012; Duguid et
al., 2020a; Vidal & Barker, 2020a), but these are typically negligibly small
in magnitude so we neglect their contribution here. On the other hand, we
ignore possible contributions to the tidal energy transfer from Reynolds
stresses involving tide-tide correlations and gradients of the convective flow
(as proposed to be the only important term for fast tides by Terquem, 2021,
even if this interpretation is a drastic over-simplification). This is because
we believe that it is not currently possible to estimate contributions from
this term without detailed numerical simulations (e.g. Barker & Astoul, 2021,
suggesting our overall conclusions regarding the ineffectiveness of
equilibrium tides are likely to hold).

Turning to wavelike tides, the tidal quality factor representing inertial wave
dissipation is computed following the (low frequency) frequency-averaged
formalism of Ogilvie (2013), and is calculated according to:

$\frac{1}{Q^{\prime}_{\mathrm{iw}}}=\frac{32\pi^{2}G}{3(2l+1)R_{\mathrm{pl}}^{2l+1}|A|^{2}}(E_{\mathrm{l}}+E_{\mathrm{l-1}}+E_{\mathrm{l+1}}),$
(2)

where the parameters $E_{\mathrm{l}}$, $E_{\mathrm{l-1}}$, and
$E_{\mathrm{l+1}}$, are specified by Eqs. (31)–(33) in B20, fully accounting
for the planetary structure. These coefficients are proportional to the
squared spin rate, implying that inertial wave dissipation is more efficient
in rapidly rotating bodies. Note that $E_{\mathrm{l}}$ and $E_{\mathrm{l\pm
1}}$ involve radial integrals that depend to some extent on the assumed core
size (inner boundary of the fluid envelope). This dependence on core size is
found to be weak in our realistic models though (which is consistent with Fig.
10 of Ogilvie, 2013, notably showing compressible polytropic models with
indices between 1 and 1.5), unlike results obtained for incompressible models.
We adopt impenetrable boundary conditions (with vanishing radial velocity) for
inertial waves (not the total tide) at the core-envelope boundary and
planetary surface, which is appropriate at low frequencies, and follows
Ogilvie 2013 and Barker 2020. We do neglect the possibility of very large
stably-stratified cores in this work though due to the uncertainties in such
planetary models.

This is a simple measure to represent the typical level of dissipation due to
inertial waves over the full range of propagation of these waves, which is
straightforward to compute in a given planetary or stellar model (e.g. Mathis,
2015). Note that this quantity is independent of the specific damping
mechanism, and is computed in a model assuming an impulsive encounter to
excite all inertial waves, which are then assumed to be subsequently fully
dissipated. Modelling tidal evolution of nearly circular or aligned orbits
using this quantity involves making assumptions, since this is not rigorously
valid (despite its wide usage in this context e.g. Bolmont & Mathis, 2016),
but it is believed to be a representative value for inertial wave dissipation.
We choose to adopt this approach (following Mathis, 2015; Bolmont & Mathis,
2016; Barker, 2020, 2022, and many others) because this measure is both
simpler and much faster to compute (hence amenable to evolutionary studies),
and it is also much more robust to the incorporation of additional (or
variation in) model physics than the direct linear (or nonlinear) response at
a particular frequency.

It should be remembered however that the actual dissipation due to inertial
waves at a given $\omega_{\mathrm{t}}$ could differ substantially from this
value (e.g. Ogilvie, 2013; Astoul & Barker, 2022, 2023). In particular,
predictions for inertial wave dissipation find substantial deviations
(potentially by orders of magnitude, either larger or smaller) at a particular
frequency – and between the frequency-averaged measure and the response at a
particular frequency – depending on the degree of density stratification
(Ogilvie, 2013), the presence of magnetic fields (Lin & Ogilvie, 2018; Wei,
2018), differential rotation (Baruteau & Rieutord, 2013; Guenel et al., 2016;
Astoul & Barker, 2022, 2023), nonlinearity (Favier et al., 2014; Barker, 2016;
Astoul & Barker, 2022, 2023), convection, varying the microscopic
diffusivities (Ogilvie & Lin, 2004; Ogilvie, 2009), and the presence of
stably-stratified (or different density) inner fluid layers (as opposed to a
rigid core) (Ogilvie, 2013; Pontin, 2022; Lin, 2023; Dewberry, 2023; Pontin et
al., 2023a). However, the frequency-averaged measure has been found to be much
more robust regarding the incorporation of magnetic fields (Lin & Ogilvie,
2018), nonlinearity, and to a limited extent differential rotation (Astoul &
Barker, 2023). It is an open question how reliable this approach will be at
modelling a population of individual systems that are each forced at a
particular tidal frequency (or range of these) at a given epoch.

We assume that upward propagating gravity waves are excited (adiabatically) at
the base of the radiative envelope and are fully damped (e.g. by radiative
diffusion or wave breaking) before propagating back to their launching sites
(e.g. Zahn, 1975; Lubow et al., 1997; Goodman & Dickson, 1998). The
corresponding tidal quality factor is then given by (B20):

$\frac{1}{Q^{\prime}_{\mathrm{gw}}}=\frac{2\left[\Gamma\left(\frac{1}{3}\right)\right]^{2}}{3^{\frac{1}{3}}(2l+1)(l(l+1))^{\frac{4}{3}}}\frac{R_{\mathrm{pl}}}{GM_{\mathrm{pl}}^{2}}\mathcal{G}|\omega_{\mathrm{t}}|^{\frac{8}{3}}.$
(3)

The quantity $\mathcal{G}$ depends on the planetary conditions at the
radiative/convective interface:

$\mathcal{G}=\sigma_{\mathrm{b}}^{2}\rho_{\mathrm{b}}r_{\mathrm{b}}^{5}\bigg{|}\frac{\operatorname{d}\\!{\,}\mathcal{N}^{2}}{\operatorname{d}\\!{\,}\ln\,r}\bigg{|}_{r=r_{\mathrm{b}}}^{-\frac{1}{3}}.$
(4)

Subscript $\mathrm{b}$ denotes the base of the radiative envelope, so
$r_{\mathrm{b}}$ and $\rho_{\mathrm{b}}$ are the corresponding radius and
density, respectively, $\mathcal{N}$ is the Brunt-Väisälä frequency, and the
parameter $\sigma_{\mathrm{b}}$ is determined numerically by the derivative of
the dynamical tide radial displacement (see Eq.(43) in B20). This is the
simplest measure of gravity wave dissipation that applies if the waves are
fully damped – regardless of the specific damping mechanism. Whether or not
this is valid is an open question. This estimate is the simplest one to
estimate the effects of gravity waves and was the one adopted in many prior
works (e.g. Lubow et al., 1997; Ogilvie & Lin, 2004) but future work should
explore in detail the validity of this assumption in thin envelopes. We omit
the influence of Coriolis forces here, partly for simplicity, and partly
because the typical level of dissipation due to gravity (or gravito-inertial)
waves is unlikely to differ substantially in this fully damped regime (e.g.
Ogilvie & Lin, 2004; Ivanov et al., 2013, though will likely differ to a
greater extent for certain tidal frequencies involving resonances with
inertial modes in neighbouring convective regions).

When a planet possesses multiple radiative and convective envelopes, the total
dissipation rates (not tidal quality factors) are derived by summing up the
contribution from each layer where the dissipation takes place.

###  2.2 Planetary model

We compute planetary models using the MESA code (version r11701; Paxton et
al., 2011, 2013, 2015, 2018, 2019). Most parameters in our inlist files are
adopted from the make planets test suit. We set initial Y and initial Z to
0.2804 and 0.02131, respectively, to reproduce the high average metallicity of
hot Jupiter hosts (<[Fe/H]> = +0.19 dex, see Petigura et al. 2018). According
to our exploration of parameter space, the tidal quality factors obtained with
metallicities in the range between -0.5 and +0.5 dex are similar within an
order of magnitude for any given age (consistent with findings for the lowest
mass stars considered in Bolmont et al., 2017). Given that abundances do not
seem to play a major role in any of our results, the effects of chemical
composition will not be reported further in this paper. The planetary mass is
varied between 0.1 and 10 $M_{\mathrm{J}}$, where $M_{\mathrm{J}}$ is the mass
of Jupiter and we fix the initial radius to $2R_{\mathrm{J}}$ (increasing to
$4R_{\mathrm{J}}$ for “hot-start" models with higher initial entropy does
produce substantial differences in $Q^{\prime}$). The core mass is 10
$M_{\mathrm{\oplus}}$, where the subscript $\oplus$ refers to Earth units, and
its density is 5 $\mathrm{g\;cm^{-3}}$, giving a fixed core radius of
approximately $0.2R_{\mathrm{J}}$. Note that the core radius, when normalised
by $R_{\mathrm{pl}}$, varies because the planetary radius (rather than the
core size) evolves in time. The fiducial value of the incident flux is 1000
$F_{\oplus}$ and we use the canonical mixing length value $\alpha=2$, though
we have explored smaller $\alpha$ values and found minimal differences in our
results. The column depth for irradiation is fixed at 330 $\rm g\;cm^{-2}$ to
reproduce the mean opacity from Guillot (2010). The control use dedt form of
energy eqn is set to .false. to avoid convergence issues at late ages. We
found that, for lower planetary masses, the default spatial resolution is too
low to provide accurate solutions for the tidal response. Therefore, we choose
a higher resolution by setting max dq = 1d-3 if there are no convergence
issues at the beginning of the MESA run. Otherwise, the parameter max dq is
gradually increased until convergence issues are avoided.

Our planetary models generate a neutrally (adiabatically) stratified interior,
as we might expect in convective regions, with only the surface layers being
stably stratified (radiative). However, observational inferences from the
Solar System’s gas giants suggest this assumption may not be valid due to
interior compositional gradients, and Jupiter or Saturn could possess extended
dilute stably stratified fluid cores (e.g. Mankovich & Fuller, 2021; Howard et
al., 2023). The consequences of interior stably stratified layers are outside
the scope of the present paper, and are currently a major uncertainty in
planetary models (but see e.g. Pontin et al., 2023b; Lin, 2023; Dewberry,
2023; Dhouib et al., 2023; Pontin et al., 2023a).

##  3 Results

###  3.1 Evolution of tidal quality factors

Figure 1: Tidal quality factors $Q^{\prime}$ for each mechanism as a function
of tidal period $P_{\mathrm{tide}}$. Solid and dashed lines correspond to the
hot ($F$ = 1000 $F_{\oplus}$) and cold ($F$ = $F_{\oplus}$) models,
respectively. The planetary mass, rotation period (fixed here to
$P_{\mathrm{rot}}=10$ hr), and age are shown on the top of each plot. Black
and blue lines represent tidal quality factors due to equilibrium tide damping
without and with rotational modification, respectively. Red and magenta lines
display tidal quality factors due to inertial and gravity waves, respectively.
Grey dotted line illustrates the minimal tidal period required for inertial
wave excitation ($P_{\mathrm{tide}}=P_{\mathrm{rot}}/2$).

We now turn to present our results for $Q^{\prime}$ computed in planetary
models. We begin by showing the dependence of the tidal quality factors for
each mechanism on the tidal forcing period $P_{\mathrm{tide}}$ in Fig. 1 for a
range of planetary masses, ages, and instellations. Here, the top, middle, and
bottom rows show planets with $M_{\mathrm{pl}}$ = 0.3, 1.0, and 10
$M_{\mathrm{J}}$, respectively. The left column corresponds to models of young
Jupiters ($t$ = 10 Myr), and the right column represents models of old
Jupiters ($t$ = 3 Gyr). The fiducial hot planets ($F$ = 1000 $F_{\oplus}$) are
shown with solid lines and cold planets ($F$ = $F_{\oplus}$) are shown with
dashed lines. These panels show a wide range of tidal frequencies, and hence
our results can be applied to model spin-orbit synchronisation, the dominant
component driving orbital circularisation, and aspects of tidal obliquity
evolution.

As reported previously, equilibrium tide dissipation is very weak in all cases
displayed according to our assumptions outlined in § 2, with the minimum
$Q^{\prime}_{\mathrm{eq}}\sim 10^{9}$. At low tidal periods, the rotationally-
modified tidal quality factor $Q^{\prime}_{\mathrm{eq,RMLT,FIT}}$, shown in
blue, is the same as the one obtained based on non-rotating convection
$Q^{\prime}_{\mathrm{eq,FIT}}$ (black, where subscript FIT refers to a fit
from numerical simulations). This is because, in the high-frequency regime,
$\nu_{\mathrm{E}}\propto\nu_{\mathrm{MLT}}\left(\frac{\omega_{\mathrm{c}}}{\omega_{\mathrm{t}}}\right)^{2}\propto\frac{u_{\mathrm{c}}^{3}}{l_{\mathrm{c}}}\frac{1}{\omega_{\mathrm{t}}^{2}}$,
and the adopted rotationally-induced scalings for the convective velocity and
length scale counteract each other (de Vries et al., 2023). In the fast tides
regime (with or without rotational inhibition of convection) – which is
typically the most relevant one in giant planets (see e.g. de Vries et al.,
2023) – we thus have $Q^{\prime}_{\mathrm{eq}}\propto P_{\mathrm{tide}}^{-1}$
because $|\omega_{\mathrm{t}}|/Q^{\prime}_{\mathrm{eq}}\propto
D_{\mathrm{v}}\propto\omega_{\mathrm{t}}^{2}\nu_{E}\propto\omega_{t}^{0}$.
Thus, the slowness of convective flows relative to tides leads to substantial
reductions in equilibrium tide damping.

With our chosen rotation period of 10 hr adopted for illustration (results for
different $P_{\mathrm{rot}}$ can be obtained simply by re-scaling since
$Q^{\prime}_{\mathrm{iw}}\propto P_{\mathrm{rot}}^{2}$, so if
$P_{\mathrm{rot}}=1$ day, $Q^{\prime}_{\mathrm{iw}}$ should be a factor of
5.76 larger), inertial wave dissipation is the dominant tidal mechanism
($Q^{\prime}_{\mathrm{iw}}$ is smallest) over the full range of tidal periods
considered, except for the ‘old’ model of 0.3 $M_{\mathrm{J}}$ planet, for
which gravity waves begin to prevail at low $P_{\mathrm{tide}}$ (i.e.,
$Q^{\prime}_{\mathrm{gw}}<Q^{\prime}_{\mathrm{iw}}$). Note that
$Q^{\prime}_{\mathrm{iw}}$ only strictly operates if the tidal period
$P_{\mathrm{tide}}>P_{\mathrm{rot}}/2$, otherwise inertial waves are not
(linearly) excited and we should not employ $Q^{\prime}_{\mathrm{iw}}$ to
model tidal evolution. This is independent of frequency because we have
adopted the frequency-averaged measure here – in reality, inertial wave
dissipation is expected to be strongly frequency-dependent, though this value
is thought to be a representative one for tidal modelling of planetary
populations as discussed in § 2. In this short tidal period regime (left
region compared to the dotted line), $Q^{\prime}_{\mathrm{gw}}$ should be used
instead according to Fig. 1.

The prediction for $Q^{\prime}_{\mathrm{iw}}$ is similar in all models since
they have the same rotation period and a similar internal structure. Indeed,
each model has a structure for all ages that is very similar to a polytrope
with a polytropic index ranging from $n=1$ (commonly thought to be appropriate
for Jupiter) to $1.5$ (thought to apply to fully convective low-mass stars).
As shown in Fig. 2, where we plot density normalised by the mean density and
radius normalised by the planetary radius, our models are well described by
Lane-Emden polytropes with such a range of $n$, where $n=1$ and $n=1.5$ are
represented by black dotted and dashed lines, respectively. The only exception
is the models of young low-mass planets displayed in red in the top panel.
Nevertheless, as these planets cool down with age, they approach the profile
of the $n=1$ polytrope. Accordingly, the models at 3 Gyr, depicted in green,
yield flatter density profiles. In contrast to the ‘cold’ models, shown with
solid lines, highly-irradiated planets, represented by dash-dotted lines, are
characterized by a steeper density gradient. At the same time, one can see
that the internal structure of more massive planets, displayed in the bottom
panel, is closer to the $n=1$ polytrope and less sensitive to age and
instellation. We, therefore, conclude that most of our models can be
approximated by polytropic solutions with $n=1$ or $n=1.5$ with sufficient
accuracy.

Adopting a polytropic model with index $n=1$ ($1.5$), we find
$Q^{\prime}_{\mathrm{iw}}=230.22\,\omega_{\mathrm{dyn}}^{2}/\Omega_{\mathrm{pl}}^{2}$
(or $130.83\,\omega_{\mathrm{dyn}}^{2}/\Omega_{\mathrm{pl}}^{2}$), where
$\omega_{\mathrm{dyn}}^{2}=GM_{\mathrm{pl}}/R_{\mathrm{pl}}^{3}$ is the
squared dynamical frequency, implying a value approximately 2558 (1454) for a
Jupiter-like model/rotation, similar to the values shown in Fig. 1. Hence
$Q^{\prime}_{\mathrm{iw}}$ (for a fixed $P_{\mathrm{rot}}$) varies only
modestly with planetary mass, age, and instellation within the ranges we
consider, ultimately because planetary structures always remain very similar
(and similar to polytropes) in our models.

Internal gravity waves become more dissipative (smaller
$Q^{\prime}_{\mathrm{gw}}$) in planets with thicker radiative envelopes,
typically corresponding to higher stellar instellations (solid purple lines),
lower planetary masses, and older ages (where planets have had time to develop
thicker envelopes). In all cases $Q^{\prime}_{\mathrm{gw}}\propto
P_{\mathrm{tide}}^{8/3}$ under the assumptions of our model, and thus shorter
tidal periods imply more efficient dissipation.

Figure 2: Density profiles as a function of radius for the planetary gaseous
envelopes of the models displayed in Fig. 1. Density is normalised by the mean
density, and radius is normalised by the planetary radius. Solid (dash-dotted)
lines correspond to the ‘cold’ (‘hot’) models; the red (green) color refers to
the age of 10 Myr (3 Gyr). Black dotted and dashed lines correspond to the
polytropic model with indices $n=1$ and $n=1.5$, respectively. Figure 3:
Evolution of tidal quality factors with planetary age in a range of models due
to equilibrium tides (first panel), inertial waves (second panel), and gravity
waves (third panel). In the fourth panel, the evolution of the planetary
radius is shown. Here, we fix tidal period $P_{\mathrm{tide}}$ = 1 day,
rotation period $P_{\mathrm{rot}}$ = 10 hr, and incident flux $F$ = 1000
$F_{\oplus}$ (i.e. these are ‘hot’ planets). Figure 4: Evolution of various
quantities that are important in computing tidal quality factor due to gravity
waves $Q^{\prime}_{\mathrm{gw}}$ according to (3) and (4) in hot planetary
models with $P_{\mathrm{tide}}$ = 1 day.

In Fig. 3, we show the evolution with planetary age of the tidal quality
factors corresponding to equilibrium tides (first panel), inertial waves
(second panel), and gravity waves (third panel) for a set of ‘hot’ planetary
models with different masses. We now fix both the tidal and rotation periods
at $P_{\mathrm{tide}}=1$ day and $P_{\mathrm{rot}}=10$ hr, respectively, to
focus on the age and mass dependence here. The choice of $P_{\mathrm{rot}}=10$
hr is made for comparison with Jupiter, but we note that for inertial waves,
we predict $Q^{\prime}_{\mathrm{iw}}\propto P_{\mathrm{rot}}^{2}$ in general,
so our results can be easily scaled for different rotation rates.

According to the prescriptions we have adopted, the equilibrium tide is
characterized by negligible damping inside the convective envelope (using
$Q^{\prime}_{\mathrm{eq,RMLT,FIT}}$), insufficient to cause a significant
change in orbital or spin parameters. The relevant tidal quality factor
$Q^{\prime}_{\mathrm{eq}}>10^{10}$ throughout the evolution of each of these
planets. This is similar to the conclusions in B20 regarding the inefficiency
of equilibrium tide dissipation in stellar interiors. Note that
$Q^{\prime}_{\mathrm{eq}}$ increases with age because the planet cools, thus
convection slows down as the planet evolves, and because the planet shrinks.

In contrast, dissipation of inertial waves appears to be the most important
mechanism in almost all cases when they are excited, with
$Q^{\prime}_{\mathrm{iw}}$ (for this $P_{\mathrm{rot}}$) ranging between
$10^{2}$ and $2\times 10^{4}$, with higher values corresponding to higher-mass
objects. As shown in the second panel of Fig. 3, planets gradually become less
dissipative with age. We have explored the reason for this, and found that it
is primarily not related to structural changes (consistent with Fig. 4 of B20
for low mass fully convective objects), but is instead explained by the
shrinking radius as the planet cools for a fixed $P_{\mathrm{rot}}$ because
$Q^{\prime}_{\mathrm{iw}}\propto\omega_{\mathrm{dyn}}^{2}/\Omega_{\mathrm{pl}}^{2}\propto
R_{\mathrm{pl}}^{-3}$. The evolution of planetary radius is displayed in the
bottom panel. Furthermore, given that planetary spin-down is typically the
natural outcome of long-term tidal and planetary evolution (unless, e.g., the
planet is spiralling into its star while remaining tidally locked), we expect
inertial wave damping to become less efficient at later epochs.

On the other hand, the tidal quality factor due to gravity waves
$Q^{\prime}_{\mathrm{gw}}$ does not exhibit substantial evolution during the
planetary lifetime, and its variation for each planetary model is within
approximately an order of magnitude for a fixed $P_{\mathrm{tide}}$ and
planetary mass. In addition to its strong dependence on $P_{\mathrm{tide}}$
($Q^{\prime}_{\mathrm{gw}}\propto P_{\mathrm{tide}}^{8/3}$), gravity wave
damping also strongly depends on the planetary mass, spanning over six orders
of magnitude for our mass range characteristic of gas giants. This is
illustrated in the third panel of Fig. 3. Similar to inertial waves, gravity
waves dissipate more efficiently (smaller $Q^{\prime}_{\mathrm{gw}}$) in
lower-mass objects, with values as small as $Q^{\prime}_{\mathrm{gw}}\approx
10^{4}$, whereas the most massive objects we consider are much less
dissipative. These values are sensitive to the radius $r_{\mathrm{b}}$ and
density $\rho_{\mathrm{b}}$ at the launching region, which is manifested
through the factor $\mathcal{G}$ being proportional to $r_{\mathrm{b}}^{5}$
(Eq. 4).

In Fig. 4 we explore in more detail the reasons for the substantial variation
of $Q^{\prime}_{\mathrm{gw}}$ with planetary mass in our models. Here, we
illustrate the evolution of the quantities involved in the expression for
$Q^{\prime}_{\mathrm{gw}}$, given by Eqs. (3) and (4), for the models with
$M_{\mathrm{pl}}$ = 0.1, 1, and 10 $M_{\mathrm{J}}$, depicted in red, green,
and blue, respectively. The top left and bottom left panels display
$\big{|}\operatorname{d}\\!{\,}\mathcal{N}^{2}/\operatorname{d}\\!{\,}\ln\,r\big{|}_{r=r_{\mathrm{b}}}$
and $\sigma_{\mathrm{b}}$, while the top right and bottom right panels display
$r_{\mathrm{b}}$ and $\rho_{\mathrm{b}}$ as a function of age, respectively.
The value of $Q^{\prime}_{\mathrm{gw}}$ obtained for planets with
$M_{\mathrm{pl}}=10M_{\mathrm{J}}$ is, on average, 3.5 orders of a magnitude
higher than for Jupiter-mass planets. These two planets are characterized by
similar values of $\sigma_{\mathrm{b}}$ and $r_{\mathrm{b}}$, but there is an
order of magnitude difference attributable to
$\big{|}\operatorname{d}\\!{\,}\mathcal{N}^{2}/\operatorname{d}\\!{\,}\ln\,r\big{|}_{r=r_{\mathrm{b}}}$
(note that this quantity is raised to the power -1/3 in Eq.(3)), and half an
order of magnitude from the variation in $\rho_{\mathrm{b}}$. The remaining
factor of $10^{2}$ results from the presence of $M_{\mathrm{pl}}^{2}$ in Eq.
(3). At the same time, decreasing the planetary mass from 1 to 0.1
$M_{\mathrm{J}}$ reduces $Q^{\prime}_{\mathrm{gw}}$ by a factor of $\sim
50-10^{2}$. A factor of $\sim 4$ arises due to differences in
$\big{|}\operatorname{d}\\!{\,}\mathcal{N}^{2}/\operatorname{d}\\!{\,}\ln\,r\big{|}_{r=r_{\mathrm{b}}}$.
One can see that $\sigma_{\mathrm{b}}$ is substantially smaller in the case of
a lower-mass planet, leading to a factor of $\sim 3$ in
$Q^{\prime}_{\mathrm{gw}}$. An additional factor of $\sim 4$ (at t $\sim 1$
Gyr) comes from differences in $\rho_{\mathrm{b}}$, and one of $\sim 1.5$
arises from the combination
$r_{\mathrm{b}}^{5}R_{\mathrm{pl}}/M_{\mathrm{pl}}^{2}$. Combining these
(crude) factors results in a total reduction of $\sim 50-10^{2}$ from 1 to 0.1
$M_{\mathrm{J}}$, in agreement with the overall differences outlined above.
Therefore, we conclude that differences in several parameters come into play
to produce the strong dependence on planetary mass exhibited by
$Q^{\prime}_{\mathrm{gw}}$ in our models, rather than any one of these
parameters.

###  3.2 Impact of the incident flux

The strength of external irradiation affects the locations of the interfaces
between convective and radiative regions near the surface, thereby altering
the conditions inside the layers where gravity waves are excited and
propagate. This is demonstrated with the example of a Jupiter-mass planet in
Fig. 5. The top and middle panels display the evolution of the tidal quality
factors of inertial and gravity waves, respectively. In the following plots,
blue circles correspond to the low incident flux, characteristic of a ‘cold’
Jupiter ($F$ = $F_{\oplus}$), and the red triangles represent a typical ‘hot’
Jupiter irradiation ($F$ = 1000 $F_{\oplus}$). One can see that the
enhancement of stellar flux received by a planet slightly increases the
inertial wave dissipation rate, which is most noticeable at early ages, which
is primarily because these planets have slightly inflated radii (see
discussion above). Nonetheless, differences between the two models are always
within an order of magnitude. On the contrary, the incident flux plays a
crucial role in modifying gravity wave damping for most of the planetary
lifetime. According to our ‘cold’ model, $Q^{\prime}_{\mathrm{gw}}$ rises by
almost two orders of magnitude up to $10^{9}$ after 30 Myrs. Eventually, the
dissipation rates are amplified to the values obtained for hot Jupiters after
6 Gyrs.

As we show in the bottom panel of Fig. 5, these dramatic changes in
$Q^{\prime}_{\mathrm{gw}}$ are directly linked to the location and number of
radiative zones near the surface that arise in these models. Here, we display
the separation of the base of each radiative zone from the planetary surface.
Note that the lower the point is, the closer the corresponding interface is to
the planetary surface. The prominent feature in the cold planetary model is
the emergence of a second radiative region after 3 Myr. Thus, the bottom base
is depicted by the blue circles, while the top base, when present, is
illustrated by the black circles. As a result of the above, both envelopes
could contribute to the dissipation of gravity waves. The overall tidal
quality factor due to gravity waves is computed (crudely) by summing up the
dissipation rates associated with each radiative region. As we mentioned
previously, the efficiency of gravity wave damping is determined by the
conditions at their launching sites, including the local density,
characterized by a steep gradient in the vicinity of the planetary surface.
Thereby, the occurrence of the outer radiative layer does not significantly
alter the evolution of the tidal quality factor as long as the bottom layer
exists. However, at the planetary age of 30 Myrs, two convective envelopes
merge into one, leaving the top radiative region as the only contributor to
gravity wave damping, which immediately manifests in a sharp increase in
$Q^{\prime}_{\mathrm{gw}}$, as shown in the middle panel. Finally, after 6
Gyrs, the radiative envelope becomes divided into two separate regions again,
allowing the dissipation rates obtained for cold and hot Jupiters to converge.
We have found that the appearance of multiple radiative layers is insensitive
to initial-Y and initial-Z in the ranges $[0.25,0.28]$ and $[0.004,0.03]$,
respectively, and the resulting values of $Q^{\prime}_{\mathrm{gw}}$ do not
differ substantially. It is possible that the emergence of these layers would
differ with different equations of state to those used in our version of MESA.

The comparison between models with different incident fluxes was also provided
in Fig. 1, where solid and dashed lines represent hot and cold Jupiters,
respectively. Contrary to the earlier epoch, at late ages, planets with
$M_{\mathrm{pl}}$ = 1 and 10 $M_{\mathrm{J}}$ reveal substantial variation in
gravity wave dissipation rate with the amount of irradiation, with hot
Jupiters being more dissipative. For 0.3 $M_{\mathrm{J}}$ planets, however,
the changes in $Q^{\prime}_{\mathrm{gw}}$ between cold and hot models are
manifested at early times. In addition, for lower-mass planets, tidal quality
factors due to equilibrium tides are also sensitive to the incident flux. In
contrast to gravity waves, the dissipation of equilibrium tides is somewhat
more effective in such low mass cold planets, which have smaller
$Q^{\prime}_{\mathrm{eq}}$.

We have assumed gravity waves are launched in each layer and are then fully
damped before returning to their launching sites, even when there multiple
radiative regions. It is uncertain whether this is justified, as is the
emergence of a second radiative region, which may either be the natural
outcome of planetary evolution for low incident fluxes, or it could be an
artifact caused by uncertainties in the equation of state – or various
neglected physics – in current versions of MESA. We have discovered that this
feature can be eliminated with the introduction of additional interior
heating, which may arise due to tidal heating or Ohmic dissipation. The fully
damped assumption can potentially be justified by linear radiative damping,
though this might only be effective in the outer zone, but there are
alternative possibilities (including differential rotation, non-linearity, and
magnetic fields). These should be explored further in future work, but for now
we caution that there are potentially large uncertainties in
$Q^{\prime}_{\mathrm{gw}}$, particularly in our cold models.

To summarise, we find gravity wave damping can be effective for highly-
irradiated planets with extended stable layers near their surfaces that are
deeper than one percent of their radii. Otherwise, inertial waves are
predicted to be the most important tidal mechanism when they are excited (i.e.
for $P_{\mathrm{tide}}>P_{\mathrm{rot}}/2$) in almost all models, except
perhaps for the latest ages for low masses when
$Q^{\prime}_{\mathrm{iw}}>Q^{\prime}_{\mathrm{gw}}$.

Figure 5: Evolution of $Q^{\prime}$ due to inertial waves (top panel) and
gravity waves (middle panel) for a Jupiter-mass planet irradiated by low ($F$
= $F_{\oplus}$) and high ($F$ = 1000 $F_{\oplus}$) incident fluxes (blue and
red colors, respectively). Bottom panel: radius of the base of each radiative
envelope with respect to the planetary surface. Blue and black circles
represent the base of the bottom and top envelopes of a cold Jupiter,
respectively. Red triangles represent the only radiative envelope of a hot
Jupiter.

##  4 Application to star-planet and planet-moon systems

Figure 6: Top panel: eccentricity distribution for observed planets as a
function of predicted tidal circularization timescale from planetary tides due
to dissipation of inertial waves. Systems with the star above (below) the
Kraft break are displayed in black (light blue). Blue dashed and black dotted
lines illustrate the mean eccentricity among the planets with $e>0.1$ (i.e.,
eccentric planets) orbiting stars below and above the Kraft break,
respectively. Histogram in the middle panel shows the average eccentricity of
the eccentric planets. Histogram in the bottom panel shows the ratio of the
number $N_{\mathrm{e>0.1}}$ of eccentric planets to the number $N$ of planets
within each bin.

The dissipation of planetary tides can potentially explain an important aspect
of the eccentricity distribution in star-planet systems, which is that hot
Jupiters tend to have smaller eccentricities (and a stronger preference for
circularity) than warm and cold Jupiters. To explore this scenario, we collect
data from the NASA Exoplanet Archive
(https://exoplanetarchive.ipac.caltech.edu/) representing massive close-in
planets ($0.1\;M_{\mathrm{J}}<M_{\mathrm{pl}}<10\;M_{\mathrm{J}}$,
$P_{\mathrm{orb}}<20\;{\rm days}$; $P_{\mathrm{orb}}$ is the orbital period).
In addition, we filter out systems containing another planet with
$P_{\mathrm{orb}}<100\;{\rm days}$ to avoid scenarios where eccentricity
excitation due to planet-planet interactions may be competing with tidal
eccentricity damping. Our main sample consists of 162 systems with a known
eccentricity, stellar effective temperature, stellar and planetary mass, and
planetary radius. This sample has been further extended by eight systems,
namely HAT-P-2, HAT-P-4, HD 118203, HD 149026, HD 189733, HD 209458,
Kepler-91, and WASP-8, with no accessible data on the planetary radii. The
radii of the corresponding planets have been derived using the mass-radius-
flux parametrization from Lazovik (2023) given by his Eqs. (16)—(18), which
encompasses the relations from Valsecchi et al. (2014) and Thorngren et al.
(2021). We also consider 136 additional planets with an upper bound on the
eccentricity below 0.1. We assume that these planets are on circular orbits
($e=0$) when calculating the relative number of eccentric planets (i.e., the
planets with $e>0.1$).

For every system in our sample, we calculate a corresponding circularization
timescale due to planetary tides, $\tau_{\mathrm{e,pl}}$, following the
equation (e.g. Goldreich & Soter, 1966):

$\tau_{\mathrm{e,pl}}=\frac{4}{63}\frac{Q^{\prime}_{\mathrm{pl}}}{n}\frac{M_{\mathrm{pl}}}{M_{\mathrm{*}}}\left(\frac{a}{R_{\mathrm{pl}}}\right)^{5}.$
(5)

Here, the tidal quality factor $Q^{\prime}_{\mathrm{pl}}$ is set equal to
$Q^{\prime}_{\mathrm{iw}}$ since, as shown in Sec. 3, inertial waves provide
the main contribution to the overall tidal dissipation inside the planets in
almost all cases in our models. $Q^{\prime}_{\mathrm{iw}}$ may be represented
as the product of two components, namely the structural tidal quality factor
$Q^{\prime}_{\mathrm{iw,s}}$ and the parameter
$\epsilon_{\mathrm{\Omega}}^{-2}\equiv\left(\Omega_{\mathrm{pl}}/\sqrt{GM_{\mathrm{pl}}/R_{\mathrm{pl}}^{3}}\right)^{-2}$.
The first component, $Q^{\prime}_{\mathrm{iw,s}}$, is computed via linear
interpolation between our models of hot planets with the adjacent masses and
ages plotted in Fig. 3 (note the large observational uncertainties in ages
does not lead to large differences according to this figure), while
$\epsilon_{\mathrm{\Omega}}$ is inferred using observational data, assuming
spin-orbit synchronization ($\Omega_{\mathrm{pl}}=n$, since this is expected
to occur more rapidly than circularization).

Our sample is illustrated in the top panel of Fig. 6. One can see that
eccentricities tend to increase with the predicted circularization timescale.
This is especially true for planets orbiting stars above the Kraft break
($T_{\mathrm{eff}}$ > 6250 K, Kraft 1967), as shown in black. Hot stars have
thin convective envelopes, leading to weaker tidal dissipation in stellar
interiors (see B20), suggesting planetary tides may be even more important for
eccentricity evolution in these systems. This weaker stellar tidal dissipation
may have several observational manifestations. In particular, star-planet
systems with a hot star typically sustain higher obliquities, as shown in
Spalding & Winn (2022); Attia et al. (2023). Here, we draw similar conclusions
concerning the eccentricity distribution, which reveals the same trend, albeit
with caution due to the low numbers involved. Indeed, among the systems with
$e>0.1$, planets orbiting stars above (below) the Kraft break have an average
eccentricity of 0.33 (0.24). Apart from this trend, a correlation between
eccentricity and eccentricity damping timescale appears to be more pronounced
in hot stars, which may imply that stellar tides also contribute to the
orbital circularization of hot Jupiters, or it could be related to the shorter
main-sequence ages of hotter stars.

As reviewed in Dawson & Johnson (2018), hot Jupiters might have formed via two
channels, namely disc migration and high-eccentricity migration (e.g.
triggered by planet-planet scattering or secular/Kozai migration). In contrast
to the former channel, the latter allows the formation of highly eccentric hot
planets. It is still unknown which channel (if any) dominates within the
overall hot Jupiter population. The presence of almost circular systems with
$\tau_{\mathrm{e,pl}}$ a few orders of magnitude higher than the age of the
universe, seen in Fig. 6, suggests that disc migration/low eccentricity
formation is likely to be a favorable scenario for some fraction of our
sample. To avoid the planets which might have formed with low initial
eccentricities, we separate the planets with $e>0.1$ (“eccentric planets") and
$e<0.1$ (“non-eccentric planets"). We select bin sizes along the $x$-axis to
have roughly equal numbers of eccentric planets in each one. For every bin, we
calculate the average eccentricity of the eccentric planets and plot it in the
middle panel of Fig. 6. In addition, we derive the fraction of planets per bin
with $e>0.1$, and we display this in the bottom panel. Both quantities are
found to increase with our predicted tidal circularization timescale. There is
only a handful of eccentric planets with $\tau_{\mathrm{e,pl}}<10^{8}$ yrs.
This might be expected if planetary tides have acted here given that the
average age of the observed systems is on the order of a few Gyr. Another
prominent detail is that the mean eccentricity of the eccentric sub-sample
increases when $\tau_{\mathrm{e,pl}}\sim 1$ Gyr, i.e., when
$\tau_{\mathrm{e,pl}}$ becomes comparable with the systems’ mean age. The
above features strongly suggest that tidal dissipation due to inertial waves
can play an important role in shaping the orbital architectures of star-planet
systems containing giant planets. On the other hand, under the assumptions of
our models, neither equilibrium tides nor gravity waves can explain tidal
circularisation timescales consistent with observations (that are shorter than
or comparable to the ages of the systems). Gravity waves could circularise
only the very closest planets under the assumptions we have made to model them
(neglecting the possibility of resonance locking, e.g. Fuller et al., 2016).

From the detailed analysis of astrometric observations of Jupiter’s and
Saturn’s satellites, Lainey et al. (2009, 2017) constrained tidal dissipation
rates in our Solar System’s giants. They inferred $k_{2}/Q=(1.102\pm
0.203)\times 10^{-5}$ for Jupiter and $k_{2}/Q=(1.59\pm 0.54)\times 10^{-4}$
for Saturn, giving approximately $Q^{\prime}=(1.59\pm 0.25)\times 10^{5}$ and
$Q^{\prime}=(9.43\pm 4.39)\times 10^{3}$ for these planets, respectively.
According to our models for evolved planets at the Solar System’s age (Figs. 1
and 3), we find $10^{3}\lesssim Q^{\prime}_{\mathrm{iw}}\lesssim 10^{4}$ for
$P_{\mathrm{rot}}\approx 10$ hr. Hence, inertial waves are sufficiently
dissipative – according to the frequency-averaged measure we have computed –
to explain observations. Given that the actual dissipation rate, and hence
$Q^{\prime}$ value, due to inertial waves at a given tidal frequency can vary
by orders of magnitude from this “typical value" represented by the frequency-
average (e.g. Ogilvie, 2013; Astoul & Barker, 2022, in both linear and
nonlinear calculations), this suggests that the orbital evolution of Jupiter’s
and Saturn’s moons can be explained by inertial waves (see also Pontin, 2022;
Lin, 2023; Dewberry, 2023; de Vries et al., 2023; Pontin et al., 2023a; Dhouib
et al., 2023). However further work is required to explore this scenario in
more detail, and to determine the validity of the frequency-averaged formalism
in modelling tidal evolution.

According to Fig. 3, the above constraints on Jupiter and Saturn can also be
obtained via gravity wave damping in the envelope. For the rotation period of
Saturn ($P_{\mathrm{rot}}=0.44$ days) and the orbital period of Enceladus
($P_{\mathrm{orb}}=1.37$ days), our Saturn model predicts
$Q^{\prime}_{\mathrm{gw}}\sim 3\times 10^{3}$ (not strongly depending on
instellation for the relevant mass and age). In turn, the present-epoch
gravity wave dissipation rate calculated for the Jupiter model strongly
depends on the incident flux. Adopting $P_{\mathrm{rot}}=0.41$ and
$P_{\mathrm{orb}}=1.77$ days (Jupiter’s rotation period and Io’s orbital
period, respectively) yields $Q^{\prime}_{\mathrm{gw}}\sim 4\times 10^{4}$ for
a hot Jupiter model and $Q^{\prime}_{\mathrm{gw}}\sim 2\times 10^{6}$ for a
cold model. Thus, our crude gravity wave dissipation estimates are also
(surprisingly) in reasonable agreement with the observations. It would be
interesting to explore the role of interior stably-stratified layers in future
work, which we have neglected in our models due to the large uncertainties
involved (i.e., to build upon Fuller et al., 2016; André et al., 2017, 2019;
Pontin, 2022; Pontin et al., 2023b, a; Lin, 2023; Dewberry, 2023; Dhouib et
al., 2023).

##  5 Conclusions

We have studied theoretically the evolution of tidal dissipation rates and
modified quality factors $Q^{\prime}$ in rotating giant planets following
their evolution using MESA interior models with masses ranging from $0.1$ to
$10M_{\mathrm{J}}$, for various incident stellar fluxes. We compute the
dissipation of equilibrium tides by rotating turbulent convection (assuming an
effective viscosity consistent with hydrodynamical simulations), dissipation
of gravity waves in the thin radiative envelope, and inertial waves in the
convective interior.

Our models indicate that inertial waves are almost always likely to be the
dominant mechanism of tidal dissipation in giant planets whenever they are
excited222These waves can also be excited “nonlinearly” by the elliptical
instability for the hottest (very shortest period) planets, which we do not
study here (e.g. de Vries et al., 2023). – i.e., when the tidal period
$P_{\mathrm{tide}}>P_{\mathrm{rot}}/2$ – and are capable of providing
$Q^{\prime}_{\mathrm{iw}}\sim 10^{3}(P_{\mathrm{rot}}/10\mathrm{hr})^{2}$.
This implies $Q^{\prime}_{\mathrm{iw}}\sim 10^{5}-10^{6}$ for orbital periods
of order $10$ days (assuming spin-orbit synchronism). Note that the frequency-
averaged measure we have adopted can differ by orders of magnitude from the
predictions at a specific tidal frequency according to linear and nonlinear
calculations (e.g. Ogilvie, 2013; Lin & Ogilvie, 2018; Astoul & Barker, 2022,
2023, and can depend on magnetic fields and differential rotation), but is
likely to represent a rather robust “typical value" of dissipation due to
these waves.

In hot low-mass planets (approx $0.1M_{\mathrm{J}}$), our models also predict
efficient dissipation of gravity waves in the radiative envelope with
$Q^{\prime}_{\mathrm{gw}}\sim 10^{4}(P_{\mathrm{tide}}/1\,\mathrm{d})^{8/3}$
(see also Lubow et al., 1997; Ogilvie & Lin, 2004). This indicates efficient
dissipation via this mechanism is also possible, though
$Q^{\prime}_{\mathrm{gw}}$ values ranging up to six orders of magnitude larger
are found in the more massive planets we modelled.

We have shown that our predicted circularization timescales from the
dissipation of inertial waves correlate well with observed planetary
eccentricities. This provides evidence that inertial wave dissipation may have
played an important role in planetary tidal evolution.

The values of $Q^{\prime}$ we have obtained can be compared with the latest
statistical inferences from modelling exoplanetary eccentricity damping in
Mahmud et al. (2023). They found $Q^{\prime}=10^{5\pm 0.5}$ for hot Jupiters
with $P_{\mathrm{tide}}\in[0.8,7]$ days, with no strong evidence of any tidal
period dependence. For eccentricity tides, assuming spin-orbit synchronism
(hence $P_{\mathrm{tide}}=P_{\mathrm{orb}}=P_{\mathrm{rot}}$), we predict
$Q^{\prime}_{\mathrm{iw}}\approx 10^{3.5}$ for $P_{\mathrm{tide}}=0.8$ days,
and $Q^{\prime}_{\mathrm{iw}}\approx 10^{5.5}$ for $P_{\mathrm{tide}}=7$ days,
with a value of $Q^{\prime}_{\mathrm{iw}}\approx 10^{4.5}$ for
$P_{\mathrm{tide}}=2.4$ days. Our results are therefore consistent with the
range they obtained for tidal periods longer than about 2.4 days, and thus we
argue that inertial waves are likely to be able to explain their results in
these cases. For the shortest tidal periods, we find more effective
dissipation than they do, though this may be mitigated when considering the
frequency-dependent tidal dissipation. It is also unclear whether their
assumption of a constant planetary radius affects their results.

Convective damping of equilibrium tides is estimated to be negligible in giant
planets compared with wavelike tides because of the strong frequency-reduction
of the effective turbulent viscosity due to the slow convection (relative to
the tide) in these bodies (Goldreich & Nicholson, 1977; Ogilvie & Lesur, 2012;
Duguid et al., 2020b; Vidal & Barker, 2020b), – though see Terquem (2021) for
an alternative viewpoint. Rapid rotation minimally affects the resulting
convective turbulent viscosities in the fast tides regime (e.g. de Vries et
al., 2023), though it reduces the effective viscosity even further for slow
tides. Further work should explore the interaction between tidal flows and
convection in more realistic numerical models.

It is essential in future work to study whether the frequency-averaged
formalism for inertial wave dissipation is appropriate to model tidal
interactions when global inertial modes are excited in realistic density-
stratified models of planets, and whether it faithfully reproduces overall
trends resulting from the dynamical evolution in a population of planetary
systems. The role of interior stably-stratified layers such as inferred for
the dilute cores of Jupiter & Saturn should also be explored further, as
should the effects of differential rotation and magnetic fields.

## Acknowledgements

We would like to thank the initial referee, Caroline Terquem, for reading two
versions of our manuscript and for providing critical comments each time, even
if we disagree with many of them, and our second referee, Jim Fuller, for a
constructive report that helped us to improve the paper. AJB was supported by
STFC grants ST/S000275/1 and ST/W000873/1. NBV was supported by EPSRC
studentship 2528559. AA was supported by a Leverhulme Trust Early Career
Fellowship (ECF-2022-362).

## Data Availability

The data underlying this article will be shared on reasonable request to the
corresponding author.

## References

  * Ahuir et al. (2021) Ahuir J., Strugarek A., Brun A.-S., Mathis S., 2021, arXiv e-prints, p. arXiv:2104.01004
  * André et al. (2017) André Q., Barker A. J., Mathis S., 2017, A&A, 605, A117
  * André et al. (2019) André Q., Mathis S., Barker A. J., 2019, A&A, 626, A82
  * Astoul & Barker (2022) Astoul A., Barker A. J., 2022, MNRAS, 516, 2913
  * Astoul & Barker (2023) Astoul A., Barker A. J., 2023, arXiv e-prints, p. arXiv:2309.02520
  * Attia et al. (2023) Attia O., Bourrier V., Delisle J. B., Eggenberger P., 2023, A&A, 674, A120
  * Barker (2016) Barker A. J., 2016, MNRAS, 459, 939
  * Barker (2020) Barker A. J., 2020, MNRAS, 498, 2270
  * Barker (2022) Barker A. J., 2022, ApJ, 927, L36
  * Barker & Astoul (2021) Barker A. J., Astoul A. A. V., 2021, MNRAS, 506, L69
  * Barker & Ogilvie (2009) Barker A. J., Ogilvie G. I., 2009, MNRAS, 395, 2268
  * Barker et al. (2014) Barker A. J., Dempsey A. M., Lithwick Y., 2014, ApJ, 791, 13
  * Baruteau & Rieutord (2013) Baruteau C., Rieutord M., 2013, Journal of Fluid Mechanics, 719, 47
  * Beaugé & Nesvorný (2012) Beaugé C., Nesvorný D., 2012, ApJ, 751, 119
  * Bolmont & Mathis (2016) Bolmont E., Mathis S., 2016, Celestial Mechanics and Dynamical Astronomy, 126, 275
  * Bolmont et al. (2012) Bolmont E., Raymond S. N., Leconte J., Matt S. P., 2012, A&A, 544, A124
  * Bolmont et al. (2017) Bolmont E., Gallet F., Mathis S., Charbonnel C., Amard L., Alibert Y., 2017, A&A, 604, A113
  * Currie et al. (2020) Currie L. K., Barker A. J., Lithwick Y., Browning M. K., 2020, MNRAS, 493, 5233
  * Dandoy et al. (2023) Dandoy V., Park J., Augustson K., Astoul A., Mathis S., 2023, A&A, 673, A6
  * Darwin (1879) Darwin G. H., 1879, Philosophical Transactions of the Royal Society of London Series I, 170, 1
  * Dawson & Johnson (2018) Dawson R. I., Johnson J. A., 2018, ARA&A, 56, 175
  * Dewberry (2023) Dewberry J. W., 2023, MNRAS, 521, 5991
  * Dhouib et al. (2023) Dhouib H., Baruteau C., Mathis S., Debras F., Astoul A., Rieutord M., 2023, arXiv e-prints, p. arXiv:2311.03288
  * Duguid et al. (2020a) Duguid C. D., Barker A. J., Jones C. A., 2020a, MNRAS, 491, 923
  * Duguid et al. (2020b) Duguid C. D., Barker A. J., Jones C. A., 2020b, MNRAS, 497, 3400
  * Favier et al. (2014) Favier B., Barker A. J., Baruteau C., Ogilvie G. I., 2014, MNRAS, 439, 845
  * Fuller et al. (2016) Fuller J., Luan J., Quataert E., 2016, MNRAS, 458, 3867
  * Gallet & Delorme (2019) Gallet F., Delorme P., 2019, A&A, 626, A120
  * Gallet et al. (2017) Gallet F., Bolmont E., Mathis S., Charbonnel C., Amard L., 2017, A&A, 604, A112
  * Gallet et al. (2018) Gallet F., Bolmont E., Bouvier J., Mathis S., Charbonnel C., 2018, A&A, 619, A80
  * Goldreich (1963) Goldreich P., 1963, MNRAS, 126, 257
  * Goldreich & Nicholson (1977) Goldreich P., Nicholson P. D., 1977, Icarus, 30, 301
  * Goldreich & Soter (1966) Goldreich P., Soter S., 1966, Icarus, 5, 375
  * Goodman & Dickson (1998) Goodman J., Dickson E. S., 1998, ApJ, 507, 938
  * Goodman & Lackner (2009) Goodman J., Lackner C., 2009, ApJ, 696, 2054
  * Guenel et al. (2016) Guenel M., Baruteau C., Mathis S., Rieutord M., 2016, A&A, 589, A22
  * Guillot (2010) Guillot T., 2010, A&A, 520, A27
  * Hamer & Schlaufman (2022) Hamer J. H., Schlaufman K. C., 2022, AJ, 164, 26
  * Howard et al. (2023) Howard S., et al., 2023, A&A, 672, A33
  * Hut (1981) Hut P., 1981, A&A, 99, 126
  * Ivanov & Papaloizou (2004) Ivanov P. B., Papaloizou J. C. B., 2004, MNRAS, 353, 1161
  * Ivanov & Papaloizou (2007) Ivanov P. B., Papaloizou J. C. B., 2007, MNRAS, 376, 682
  * Ivanov & Papaloizou (2010) Ivanov P. B., Papaloizou J. C. B., 2010, MNRAS, 407, 1609
  * Ivanov et al. (2013) Ivanov P. B., Papaloizou J. C. B., Chernov S. V., 2013, MNRAS, 432, 2339
  * Jackson et al. (2008) Jackson B., Greenberg R., Barnes R., 2008, ApJ, 678, 1396
  * Jackson et al. (2016) Jackson B., Jensen E., Peacock S., Arras P., Penev K., 2016, Celestial Mechanics and Dynamical Astronomy, 126, 227
  * Kraft (1967) Kraft R. P., 1967, ApJ, 150, 551
  * Lai (2012) Lai D., 2012, MNRAS, 423, 486
  * Lainey et al. (2009) Lainey V., Arlot J.-E., Karatekin Ö., van Hoolst T., 2009, Nature, 459, 957
  * Lainey et al. (2017) Lainey V., et al., 2017, Icarus, 281, 286
  * Lazovik (2021) Lazovik Y. A., 2021, MNRAS, 508, 3408
  * Lazovik (2023) Lazovik Y. A., 2023, MNRAS, 520, 3749
  * Lin (2023) Lin Y., 2023, A&A, 671, A37
  * Lin & Ogilvie (2018) Lin Y., Ogilvie G. I., 2018, MNRAS, 474, 1644
  * Lubow et al. (1997) Lubow S. H., Tout C. A., Livio M., 1997, ApJ, 484, 866
  * Mahmud et al. (2023) Mahmud M. M., Penev K. M., Schussler J. A., 2023, Monthly Notices of the Royal Astronomical Society, p. stad2298
  * Mankovich & Fuller (2021) Mankovich C. R., Fuller J., 2021, Nature Astronomy, 5, 1103
  * Mathis (2015) Mathis S., 2015, A&A, 580, L3
  * Mathis et al. (2016) Mathis S., Auclair-Desrotour P., Guenel M., Gallet F., Le Poncin-Lafitte C., 2016, A&A, 592, A33
  * Nagasawa et al. (2008) Nagasawa M., Ida S., Bessho T., 2008, ApJ, 678, 498
  * Ogilvie (2009) Ogilvie G. I., 2009, MNRAS, 396, 794
  * Ogilvie (2013) Ogilvie G. I., 2013, MNRAS, 429, 613
  * Ogilvie (2014) Ogilvie G. I., 2014, ARA&A, 52, 171
  * Ogilvie & Lesur (2012) Ogilvie G. I., Lesur G., 2012, MNRAS, 422, 1975
  * Ogilvie & Lin (2004) Ogilvie G. I., Lin D. N. C., 2004, ApJ, 610, 477
  * Ogilvie & Lin (2007) Ogilvie G. I., Lin D. N. C., 2007, ApJ, 661, 1180
  * Papaloizou & Ivanov (2010) Papaloizou J. C. B., Ivanov P. B., 2010, MNRAS, 407, 1631
  * Paxton et al. (2011) Paxton B., Bildsten L., Dotter A., Herwig F., Lesaffre P., Timmes F., 2011, ApJS, 192, 3
  * Paxton et al. (2013) Paxton B., et al., 2013, ApJS, 208, 4
  * Paxton et al. (2015) Paxton B., et al., 2015, ApJS, 220, 15
  * Paxton et al. (2018) Paxton B., et al., 2018, ApJS, 234, 34
  * Paxton et al. (2019) Paxton B., et al., 2019, ApJS, 243, 10
  * Penev et al. (2018) Penev K., Bouma L. G., Winn J. N., Hartman J. D., 2018, AJ, 155, 165
  * Petigura et al. (2018) Petigura E. A., et al., 2018, AJ, 155, 89
  * Pontin (2022) Pontin C. M., 2022, Wave propagation and tidal dissipation in giant planets containing regions of stable stratification, https://etheses.whiterose.ac.uk/31119/
  * Pontin et al. (2023a) Pontin C. M., Barker A. J., Hollerbach R., 2023a, arXiv e-prints, p. arXiv:2311.03273
  * Pontin et al. (2023b) Pontin C. M., Barker A. J., Hollerbach R., 2023b, ApJ, 950, 176
  * Spalding & Winn (2022) Spalding C., Winn J. N., 2022, ApJ, 927, 22
  * Stevenson (1979) Stevenson D. J., 1979, Geophysical and Astrophysical Fluid Dynamics, 12, 139
  * Terquem (2021) Terquem C., 2021, MNRAS, 503, 5789
  * Terquem & Martin (2021) Terquem C., Martin S., 2021, MNRAS, 507, 4165
  * Terquem et al. (1998) Terquem C., Papaloizou J. C. B., Nelson R. P., Lin D. N. C., 1998, ApJ, 502, 788
  * Thorngren et al. (2021) Thorngren D. P., Fortney J. J., Lopez E. D., Berger T. A., Huber D., 2021, ApJ, 909, L16
  * Valsecchi et al. (2014) Valsecchi F., Rasio F. A., Steffen J. H., 2014, ApJ, 793, L3
  * Vidal & Barker (2020a) Vidal J., Barker A. J., 2020a, MNRAS, 497, 4472
  * Vidal & Barker (2020b) Vidal J., Barker A. J., 2020b, ApJ, 888, L31
  * Villaver & Livio (2009) Villaver E., Livio M., 2009, ApJ, 705, L81
  * Wei (2018) Wei X., 2018, ApJ, 854, 34
  * Weinberg et al. (2012) Weinberg N. N., Arras P., Quataert E., Burkart J., 2012, ApJ, 751, 136
  * Winn et al. (2010) Winn J. N., Fabrycky D., Albrecht S., Johnson J. A., 2010, ApJ, 718, L145
  * Witte & Savonije (2002) Witte M. G., Savonije G. J., 2002, A&A, 386, 222
  * Wu (2005a) Wu Y., 2005a, ApJ, 635, 674
  * Wu (2005b) Wu Y., 2005b, ApJ, 635, 688
  * Zahn (1966) Zahn J. P., 1966, Annales d’Astrophysique, 29, 313
  * Zahn (1975) Zahn J. P., 1975, A&A, 41, 329
  * Zahn (1977) Zahn J. P., 1977, A&A, 500, 121
  * Zahn (1989) Zahn J. P., 1989, A&A, 220, 112
  * de Vries et al. (2023) de Vries N. B., Barker A. J., Hollerbach R., 2023, MNRAS, 524, 2661