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

We study tidal dissipation in models of rotating giant planets with masses in the range $0.1 - 10 M_\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'$ 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'\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.

Some early steps towards a theory of tides were made long before ★ E-mail: yaroslav.lazovik@gmail.com† E-mail: A.J.Barker@leeds.ac.uk ‡ E-mail: mmnbdv@leeds.ac.uk § E-mail: A.A.V.Astoul@leeds.ac.uk 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 (1975Zahn ( , 1977Zahn ( , 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  ′ (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 (2021Lazovik ( , 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 andIvanov &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  ′ 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 10 J computed with the MESA code (Paxton et al. 2011(Paxton et al. , 2013(Paxton et al. , 2015(Paxton et al. , 2018(Paxton et al. , 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.

Tidal dissipation mechanisms
We consider a giant planet of mass  pl and radius  pl and employ spherical coordinates centred on the body with radial coordinate .The intensity of tidal dissipation is often quantified by the (modified) tidal quality factor  ′ , 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  ′ .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 ( ′ eq ,  ′ iw , and  ′ 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  = 2 and azimuthal wavenumber  = 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: where  is the gravitational constant,  v is the rate of viscous dissipation of the equilibrium tide, and  is the amplitude of the tidal potential component (this will not be specified further as it cancels because  v ∝ | | 2 in linear theory).The tidal forcing frequency is  t = 2/ tide (where  tide is the tidal period), which is  t = 2(−Ω pl ) for a circular, aligned orbit, where  and Ω pl are the orbital mean motion and planetary spin, respectively.For an eccentric orbit with synchronised1 (and aligned) spin we instead have The viscous dissipation rate  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  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  MLT ∝  c  c , with  c the convective velocity and  c =   the mixing-length ( is mixing-length parameter and   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),  E is reduced for fast tides (where  t >   , with the convective frequency  c =  c / c ) in a frequency-dependent manner.This can be accounted for using a piece-wise continuous correction factor depending on the ratio  t / 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) andde Vries et al. (2023) have shown that such rapid rotation reduces  E even further (though interestingly the regime with  t ≫  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  c and  c at each radius by Ro 3/5 and Ro 1/5 , respectively, where Ro is the convective Rossby number (Ro =  c /Ω 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   -corresponding to tidal antidissipation -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 oversimplification).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: where the parameters  l ,  l−1 , and  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  l and  l±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 2020Barker , 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  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): ( The quantity G depends on the planetary conditions at the radiative/convective interface: Subscript b denotes the base of the radiative envelope, so  b and  b are the corresponding radius and density, respectively, N is the Brunt-Väisälä frequency, and the parameter  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 gravitoinertial) 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.

Planetary model
We compute planetary models using the MESA code (version r11701; Paxton et al. 2011Paxton et al. , 2013Paxton et al. , 2015Paxton et al. , 2018Paxton et al. , 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.19dex, 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  J , where  J is the mass of Jupiter and we fix the initial radius to 2 J (increasing to 4 J for "hot-start" models with higher initial entropy does produce substantial differences in  ′ ).The core mass is 10  ⊕ , where the subscript ⊕ refers to Earth units, and its density is 5 g cm −3 , giving a fixed core radius of approximately 0.2 J .Note that the core radius, when normalised by  pl , varies because the planetary radius (rather than the core size) evolves in time.The fiducial value of the incident flux is 1000  ⊕ and we use the canonical mixing length value  = 2, though we have explored smaller  values and found minimal differences in our results.The column depth for irradiation is fixed at 330 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).

Evolution of tidal quality factors
We now turn to present our results for  ′ computed in planetary models.We begin by showing the dependence of the tidal quality factors for each mechanism on the tidal forcing period  tide in Fig. 1 for a range of planetary masses, ages, and instellations.Here, the top, middle, and bottom rows show planets with  pl = 0.3, 1.0, and 10  J , respectively.The left column corresponds to models of young Jupiters ( = 10 Myr), and the right column represents models of old Jupiters ( = 3 Gyr).The fiducial hot planets ( = 1000  ⊕ ) are shown with solid lines and cold planets ( =  ⊕ ) 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  ′ eq ∼ 10 9 .At low tidal periods, the rotationally-modified tidal quality factor  ′ eq,RMLT,FIT , shown in blue, is the same as the one obtained based on non-rotating convection  ′ eq,FIT (black, where subscript FIT refers to a fit from numerical simulations).This is because, in the high-frequency regime, , 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  ′ eq ∝  −1 tide because | t |/ ′ eq ∝  v ∝  2 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  rot can be obtained simply by re-scaling since  ′ iw ∝  2 rot , so if  rot = 1 day,  ′ iw should be a factor of 5.76 larger), inertial wave dissipation is the dominant tidal mechanism ( ′ iw is smallest) over the full range of tidal periods considered, except for the 'old' model of 0.3  J planet, for which gravity waves begin to prevail at low  tide (i.e.,  ′ gw <  ′ iw ).Note that  ′ iw only strictly operates if the tidal period  tide >  rot /2, otherwise inertial waves are not (linearly) excited and we should not employ  ′ 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),  ′ gw should be used instead according to Fig. 1.
The prediction for  ′ 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  = 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 , where  = 1 and  = 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  = 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, highlyirradiated 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  = 1 polytrope and less sensitive to age and instellation.We, therefore, conclude that most of our models can be approximated by polytropic solutions with  = 1 or  = 1.5 with sufficient accuracy.
Adopting a polytropic model with index  = 1 (1.5), we find  ′ iw = 230.22 2 dyn /Ω 2 pl (or 130.83  2 dyn /Ω 2 pl ), where  2 dyn =   pl / 3 pl 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  ′ iw (for a fixed  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  ′ 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  ′ gw ∝  8/3 tide under the assumptions of our model, and thus shorter tidal periods imply more efficient dissipation.
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  tide = 1 day and  rot = 10 hr, respec- tively, to focus on the age and mass dependence here.The choice of  rot = 10 hr is made for comparison with Jupiter, but we note that for inertial waves, we predict  ′ iw ∝  2 rot 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  ′ eq,RMLT,FIT ), insufficient to cause a significant change in orbital or spin parameters.The relevant tidal quality factor  ′ 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  ′ 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  ′ iw (for this  rot ) ranging between 10 2 and 2×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  rot because  ′ iw ∝  2 dyn /Ω 2 pl ∝  −3 pl .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  ′ 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  tide and planetary mass.In addition to its strong dependence on  tide ( ′ gw ∝  8/3 tide ), 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  ′ gw ) in lower-mass objects, with values as small as  ′ gw ≈ 10 4 , whereas the most massive objects we consider are much less dissipative.These values are sensitive to the radius  b and density  b at the launching region, which is manifested through the factor G being proportional to  5 b (Eq.4).
In Fig. 4 we explore in more detail the reasons for the substantial variation of  ′ gw with planetary mass in our models.Here, we illustrate the evolution of the quantities involved in the expression for  ′ gw , given by Eqs. ( 3) and ( 4), for the models with  pl = 0.1, 1, and 10  J , depicted in red, green, and blue, respectively.The top left and bottom left panels display d N 2 /d ln  = b and  b , while the top right and bottom right panels display  b and  b as a function of age, respectively.The value of  ′ gw obtained for planets with  pl = 10 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  b and  b , but there is an order of magnitude difference attributable to d N 2 /d ln  = 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  b .The remaining factor of 10 2 results from the presence of  2 pl in Eq. (3).At the same time, decreasing the plan-etary mass from 1 to 0.1  J reduces  ′ gw by a factor of ∼ 50 − 10 2 .A factor of ∼ 4 arises due to differences in d N 2 /d ln  = b .One can see that  b is substantially smaller in the case of a lower-mass planet, leading to a factor of ∼ 3 in  ′ gw .An additional factor of ∼ 4 (at t ∼ 1 Gyr) comes from differences in  b , and one of ∼ 1.5 arises from the combination  5 b  pl / 2 pl .Combining these (crude) factors results in a total reduction of ∼ 50 − 10 2 from 1 to 0.1  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  ′ gw in our models, rather than any one of these parameters.

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 ( =  ⊕ ), and the red triangles represent a typical 'hot' Jupiter irradiation ( = 1000  ⊕ ).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,  ′ 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  ′ 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  ′ 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  ′ gw do not differ substantially.It is possible that the emer-gence 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  pl = 1 and 10  J reveal substantial variation in gravity wave dissipation rate with the amount of irradiation, with hot Jupiters being more dissipative.For 0.3  J planets, however, the changes in  ′ 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  ′ 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  ′ 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  tide >  rot /2) in almost all models, except perhaps for the latest ages for low masses when  ′ iw >  ′ gw .

APPLICATION TO STAR-PLANET AND PLANET-MOON SYSTEMS
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  J <  pl < 10  J ,  orb < 20 days;  orb is the orbital period).In addition, we filter out systems containing another planet with  orb < 100 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 ( = 0) when  calculating the relative number of eccentric planets (i.e., the planets with  > 0.1).
For every system in our sample, we calculate a corresponding circularization timescale due to planetary tides,  e,pl , following the equation (e.g.Goldreich & Soter 1966): (5) Here, the tidal quality factor  ′ pl is set equal to  ′ 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. ′ iw may be represented as the product of two components, namely the structural tidal quality factor  ′ iw,s and the parameter . The first component,  ′ 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  Ω is inferred using observational data, assuming spin-orbit synchronization (Ω pl = , 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 ( 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  > 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 higheccentricity 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  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  > 0.1 ("eccentric planets") and  < 0.1 ("non-eccentric planets").We select bin sizes along the -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  > 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  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  e,pl ∼ 1 Gyr, i.e., when  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. (2009Lainey et al. ( , 2017) constrained tidal dissipation rates in our Solar System's giants.They inferred  2 / = (1.102± 0.203) × 10 −5 for Jupiter and  2 / = (1.59 ± 0.54) × 10 −4 for Saturn, giving approximately  ′ = (1.59 ± 0.25) × 10 5 and  ′ = (9.43 ± 4.39) × 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 ≲  ′ iw ≲ 10 4 for  rot ≈ 10 hr.Hence, iner-tial waves are sufficiently dissipative -according to the frequencyaveraged measure we have computed -to explain observations.Given that the actual dissipation rate, and hence  ′ 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.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 ( rot = 0.44 days) and the orbital period of Enceladus ( orb = 1.37 days), our Saturn model predicts  ′ gw ∼ 3×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  rot = 0.41 and  orb = 1.77 days (Jupiter's rotation period and Io's orbital period, respectively) yields  ′ gw ∼ 4 × 10 4 for a hot Jupiter model and  ′ gw ∼ 2 × 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. 2017André et al. , 2019;;Pontin 2022;Pontin et al. 2023b,a;Lin 2023;Dewberry 2023;Dhouib et al. 2023).

CONCLUSIONS
We have studied theoretically the evolution of tidal dissipation rates and modified quality factors  ′ in rotating giant planets following their evolution using MESA interior models with masses ranging from 0.1 to 10 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 excited2 -i.e., when the tidal period  tide >  rot /2 -and are capable of providing  ′ iw ∼ 10 3 ( rot /10hr) 2 .This implies  ′ iw ∼ 10 5 − 10 6 for orbital periods of order 10 days (assuming spinorbit 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.1 J ), our models also predict efficient dissipation of gravity waves in the radiative envelope with  ′ gw ∼ 10 4 ( tide /1 d) 8/3 (see also Lubow et al. 1997;Ogilvie & Lin 2004).This indicates efficient dissipation via this mechanism is also possible, though  ′ 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  ′ we have obtained can be compared with the latest statistical inferences from modelling exoplanetary eccentricity damping in Mahmud et al. (2023).They found  ′ = 10 5±0.5 for hot Jupiters with  tide ∈ [0.8, 7] days, with no strong evidence of any tidal period dependence.For eccentricity tides, assuming spin-orbit synchronism (hence  tide =  orb =  rot ), we predict  ′ iw ≈ 10 3.5 for  tide = 0.8 days, and  ′ iw ≈ 10 5.5 for  tide = 7 days, with a value of  ′ iw ≈ 10 4.5 for  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 frequencyaveraged 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.

Figure 1 .
Figure1.Tidal quality factors  ′ for each mechanism as a function of tidal period  tide .Solid and dashed lines correspond to the hot ( = 1000  ⊕ ) and cold ( =  ⊕ ) models, respectively.The planetary mass, rotation period (fixed here to  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 ( tide =  rot /2).

FFigure 2 .
Figure2.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 (dashdotted) 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  = 1 and  = 1.5, respectively.

Figure 4 .
Figure 4. Evolution of various quantities that are important in computing tidal quality factor due to gravity waves  ′ gw according to (3) and (4) in hot planetary models with  tide = 1 day.

Figure 5 .
Figure 5. Evolution of  ′ due to inertial waves (top panel) and gravity waves (middle panel) for a Jupiter-mass planet irradiated by low ( =  ⊕ ) and high ( = 1000  ⊕ ) 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.

Figure 6 .
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  > 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  e>0.1 of eccentric planets to the number  of planets within each bin.