Orbital Decay of Hot Jupiters due to Weakly Nonlinear Tidal Dissipation

We study tidal dissipation in hot Jupiter host stars due to the nonlinear damping of tidally driven g-modes, extending the calculations of Essick & Weinberg to a wide variety of stellar host types. This process causes the planet’s orbit to decay and has potentially important consequences for the evolution and fate of hot Jupiters. Previous studies either only accounted for linear dissipation processes or assumed that the resonantly excited primary mode becomes strongly nonlinear and breaks as it approaches the stellar center. However, the great majority of hot Jupiter systems are in the weakly nonlinear regime in which the primary mode does not break but instead excites a sea of secondary modes via three-mode interactions. We simulate these nonlinear interactions and calculate the net mode dissipation for stars that range in mass from 0.5M ⊙ ≤ M ⋆ ≤ 2.0M ⊙ and in age from the early main sequence to the subgiant phase. We find that the nonlinearly excited secondary modes can enhance the tidal dissipation by orders of magnitude compared to linear dissipation processes. For the stars with M ⋆ ≲ 1.0M ⊙ of nearly any age, we find that the orbital decay time is ≲100 Myr for orbital periods P orb ≲ 1 day. For M ⋆ ≳ 1.2M ⊙, the orbital decay time only becomes short on the subgiant branch, where it can be ≲10 Myr for P orb ≲ 2 days and result in significant transit time shifts. We discuss these results in the context of known hot Jupiter systems and examine the prospects for detecting their orbital decay with transit timing measurements.


Introduction
The orbits of hot Jupiters decay over time due to the tide-induced transfer of energy and angular momentum from the orbit to the host star.The orbital decay rate depends on the efficiency of tidal dissipation within the star and is sensitive to its structure and evolutionary state.The rate can therefore be a strong function of not just orbital period, but also stellar mass and age (for a review of tidal dissipation processes in stars and giant planets, see Ogilvie 2014).
The most direct observational evidence of hot Jupiter orbital decay comes from the measured transit time shifts of WASP-12b (Maciejewski et al. 2016;Patra et al. 2017;Yee et al. 2020) and Kepler-1658b (Vissapragada Corresponding author: Nevin N. Weinberg nevin@uta.eduet al. 2022).A number of studies also find evidence from the statistical analysis of hot Jupiter populations (Jackson et al. 2009;Teitler & Königl 2014;Penev et al. 2018;Hamer & Schlaufman 2019).For example, Jackson et al. (2009) find that older planets tend to be farther from their hosts than younger planets, which they argue is evidence for the ongoing destruction of planets by tides.Using the measured Galactic velocity dispersion, Hamer & Schlaufman (2019) show that hot Jupiter host stars are preferentially younger than a matched sample of field stars, which can be explained if the planets are destroyed while the hosts are on the main sequence.McQuillan et al. (2013) find a dearth of close-in planets orbiting rapidly rotating stars, which Teitler & Königl (2014) attribute to tidal ingestion of giant planets.The recently reported infrared transient ZTF SLRN-2020 appears to capture the moments of a planet's ingestion by a main sequence or early subgiant branch star with mass around 0.8 − 1.5M , and could be the culmination of tide-induced orbital decay (De et al. 2023).
The dominant source of dissipation in most hot Jupiter host stars is damping of the resonantly excited internal gravity waves that comprise the dynamical tide (Barker & Ogilvie 2010;Ivanov et al. 2013;Essick & Weinberg 2016;Barker 2020).In stars with thick outer convective envelopes, an internal gravity wave is excited near the radiative-convective boundary and propagates inwards towards the core.As the wave approaches the stellar center, its amplitudes grows due to geometric focusing and it can become nonlinear.If the amplitude of the wave is not too large, it reflects at an inner turning point and forms a standing wave, i.e., a g-mode.
Most studies either ignore nonlinearities and treat the wave as linear throughout the star or they take the other extreme and assume it becomes strongly nonlinear and undergoes wave breaking in the core.However, as we show (see also Barker 2020), most hot Jupiter systems are in an intermediate regime where the wave is weakly nonlinear.In this regime, the wave excites a sea of secondary waves through nonlinear wave-wave interactions.Since the secondary waves have much shorter wavelengths than the primary wave, they have much larger damping rates.Treating the primary wave as linear is therefore not only invalid, it can greatly underestimate the efficiency of tidal dissipation.On the other hand, treating it as strongly nonlinear overestimates the efficiency because it assumes the wave transfers all of its energy and angular momentum on its first journey into the stellar center.In the weakly nonlinear regime, by contrast, the primary wave deposits only a fraction of its energy and angular momentum each journey.The value of the deposited fraction depends on the detailed interaction between the primary and secondary waves and determines the rate of tidal dissipation.
The first calculation of hot Jupiter orbital decay in the weakly nonlinear regime was by Essick & Weinberg (2016;hereafter EW16).They assumed a solartype host star and considered a range of planet masses M p and orbital periods P orb .For a solar-type host, the tide is in the weakly nonlinear regime for M p 3.6M J (P orb /day) −0.1 ; above this planet mass the waves are strongly nonlinear (Barker & Ogilvie 2011).By solving the dynamics of large networks of nonlinearly interacting waves, EW16 calculated the weakly nonlinear tidal dissipation and found a stellar tidal quality factor Q 3 × 10 5 (M p /M J ) 0.5 (P orb /day) 2.4 .The dissipation can thus be highly efficient, causing hot Jupiters with solar-type hosts and P orb 2 days to decay on timescales that are small compared to the main-sequence lifetime of their host star.They found that the decrease in P orb associated with these short decay times could produce detectable transit timing variations.EW16 only considered hot Jupiters orbiting a solartype star.However, as we show in Figure 1, there are many hot Jupiters orbiting non-solar type stars.The results of EW16 do not inform these other systems since the efficiency of tidal dissipation is often sensitive to stellar structure and therefore stellar mass M and age.For example, on the main sequence the mass of the convective envelope rapidly decreases with increasing M , and a convective core appears above M ≈ 1.2M .While the former effect may lead to only modest changes in dissipation rate for M 1.2M , the appearance of a convective core may drastically decrease the degree of wave nonlinearity, causing a sudden decrease in tidal dissipation relative to slightly lower mass stars.As we show, the efficiency of tidal dissipation can also be sensitive to stellar age, especially as stars begin to evolve up the subgiant branch.
The paper is organized as follows.In Section 2, we describe the formalism we use to study weakly nonlinear tidal dissipation, including the equation of motion for our mode decomposition.In Section 3, we present our stellar models and their linear and nonlinear mode properties.In Section 4, we describe our method for building and integrating networks of nonlinearly interacting modes.In Section 5, we present the main results of our calculations, including how the orbital decay rate and transit time shift depend on M , stellar age, M p , and P orb .In Section 6, we discuss the implications of In hot Jupiter systems, the tide raised by the planet excites high-order g-modes in the host star through linear and, we assume, weakly nonlinear forces.Since the planet's orbital period will in general be much shorter than the rotational period of the star, the g-modes are not strongly modified by Coriolis forces and we therefore neglect the star's rotation. 1We also assume that the planet's orbit is circular, as is the case for most of the observed hot Jupiters. 2The planet's rotation should then be synchronous with the orbit (Storch & Lai 2014;Barker & Lithwick 2014) and thus there is no tidal dissipation within the planet.We study this problem using the formalism developed in Weinberg et al. (2012; hereafter WAQB; see also Schenk et al. 2001;Van Hoolst 1994).Here we summarize the method and refer the reader to WAQB for a more detailed discussion.

Mode amplitude equation
The position x of a fluid element in the perturbed star at time t is related to its position x in the unperturbed star by x = x + ξ(x, t), where ξ(x, t) is the Lagrangian displacement vector.The equation of motion for ξ to lowest nonlinear order is where ρ is the background density, f 1 and f 2 are the linear and leading-order nonlinear restoring forces, is the tidal acceleration, and U is the tidal potential.We include only the dominant l = 2 tidal harmonic and since we assume that the orbit is circular, in a spherical coordinate system (r, θ, φ) centered on the star.Here = (M p /M )(R /a)3 is the tidal 1 Hot stars on the main sequence with M 1.2M and effective temperatures T eff ≥ 6250 K (i.e., above the Kraft break; Kraft 1967) can rotate with spin periods of a few days, comparable to P orb of hot Jupiters.However, we do not include such stars in our study since their tidal dissipation rate is negligible owing to the presence of a convective core, as explained in Section 3.1.The stars with M ≥ 1.2M we consider are all on the subgiant branch and are observed to have spin periods 10 days (Avallone et al. 2022). 2 https://exoplanetarchive.ipac.caltech.edu/strength, ω = (GM /R 3 ) 1/2 is the dynamical frequency of a star with mass M and radius R , a and Ω = G(M + M p )/a 3 are the orbital semi-major axis and frequency, and W 20 = −(π/5) 1/2 , W 2±2 = (3π/10) 1/2 , W 2±1 = 0. We use the linear eigenmodes to expand the six-dimensional phase-space vector as (Schenk et al. 2001) where each eigenmode is specified by its amplitude q a , frequency ω a , and eigenfunction ξ a .The sum over a runs over all mode quantum numbers (radial order n a , angular degree l a , and azimuthal order m a ) and frequency signs to allow both a mode and its complex conjugate.We normalize the eigenmodes as so that the energy of a mode is E a (t) = q * a (t)q a (t)E . 3Substituting Equation (4) into Equation (1), adding a linear damping term, and using the orthogonality of the eigenmodes leads to a coupled, nonlinear amplitude equation for each mode where The coefficient γ a is the linear damping rate of the mode, U a and U ab are the linear and nonlinear tidal coefficients, and κ abc is the three-mode coupling coefficient.
The nonlinear tide cancels significantly with threemode coupling to the equilibrium tide such that U ab + 2 c κ abc U c 0 (WAQB).By treating the cancellation as perfect, we can write Equation (6) as We solve this coupled set of equations to find the amplitude evolution q a (t) of each mode in our network and from this determine the total tidal dissipation rate (EW16) In Section 3, we present the values of γ a , U a , and κ abc for our stellar models.

Parametric instability
As discussed in WAQB and EW16, in the absence of nonlinear mode coupling the linear tide drives a mode to an energy where ∆ a = ω a − m a Ω is the linear detuning.This 'parent' mode is unstable to the parametric instability if there exists a pair of 'daughter' modes b and c such that E lin E thr , where the threshold energy is and ∆ abc = ω b + ω c + m a Ω is the nonlinear detuning.
If the parent energy E a E thr , the daughters grow exponentially at a rate while their energy E b , E c E a .Eventually, the system reaches a nonlinear equilibrium as described in Appendix B of EW16.
In the host star of a hot Jupiter, the linearly resonant parents have E lin E thr and they can therefore excite many daughter pairs (WAQB).For example, in a solar model even a 0.1 M J companion in a 3 day orbit has ∼ 10 3 daughter pairs for which E lin > E thr (see Figure 1 in EW16).Many of these daughters are nonlinearly driven to such large amplitudes that they in turn excite many granddaughters to large amplitudes, and the granddaughters excite great-granddaughters, and so on.The total number of unstable modes is therefore very large and the system's nonlinear equilibrium is much more complicated than the simple, analytic three-mode equilibrium given in Appendix B of EW16.In Section 4, we describe our procedure for building and integrating these large networks of nonlinearly excited modes.

Orbital decay rate
Dissipation of the tidally excited modes removes energy from the orbit and as a result the planet inspirals.As in EW16, we assume that linear damping of the waves excited within the star is the only dissipation in the system.Although the rotational energies of the star and the synchronized planet increase as the orbit decays, these changes are small compared to the corresponding change in orbital energy.Similarly, the energy in the excited stellar modes themselves and in their interaction energy may change with orbital period, but these too are small effects.
Since | Ė/E orb | Ω, where E orb = −GM M p /2a is the orbital energy, we model the back-reaction on the orbit as a steady decrease in E orb of quasi-Keplerian circular orbits.At each orbital period P orb = 2π/Ω, the timescale of orbital energy decay is then given by We can compute a corresponding time-averaged decay time where Ė is the time-averaged energy dissipation rate.
In our calculations, we average over timescales of ∼ 10 6 P orb , which is much longer than the time it takes for the system to reach a nonlinear equilibrium but much less than the orbital decay timescales.It is common to parameterize τ in terms of the star's tidal quality factor (Goldreich & Soter 1966; see also Jackson et al. 2008) where the expression assumes a circular orbit.

Transit time shift
The decrease in P orb due to tide-induced orbital decay will cause a planet's observed transits to arrive early.Over a duration T dur , the transit time will shift by an amount (Birkby et al. 2014) Based on the typical uncertainties in transit timing measurements of hot Jupiter systems (see, e.g., Patra et al.

Stellar Models and their Mode Properties
Section 3.1 presents the set of stellar models we use in our calculations.Section 3.2 discusses the dependence of mode displacement ξ r on radius and stellar structure and the conditions under which wave breaking occurs.Section 3.3 describes our calculations of the various coefficients in the amplitude equation and presents their values as a function of the stellar and mode properties.

Stellar models
As we showed in Figure 1, most of the hot Jupiter host stars are in the mass range 0.5 M /M 2.0 and span a range of ages (the age is often quite uncertain).Motivated by these observations, we use the MESA stellar evolution code (Paxton et al. 2011(Paxton et al. , 2013(Paxton et al. , 2015(Paxton et al. , 2018(Paxton et al. , 2019;;Jermyn et al. 2022) to construct stellar models with M and age as listed in Table 1.The key parameters of the MESA inlist files we use to build the models are provided in Appendix A.
Our models range in age from the early main sequence through the subgiant branch, and all have a radiative interior and a convective envelope.For M ≥ 1.2M , we only consider stars on the subgiant branch; nonlinear mode coupling of g-modes is negligible in the interior of pre-main sequence stars or main sequence stars with M 1.2M .This is because they have convective cores and their g-modes therefore do not steepen sufficiently, as explained below, before reflecting at the radiativeconvective boundary.

Wave steepening and breaking
High-order g-modes are restored by buoyancy and propagate between inner and outer turning points determined by the locations at which ω a N (r), where N is the Brunt-Väisälä buoyancy frequency (Aerts et al. 2010).Their inner turning points, r inner , are very close to the stellar center and their outer turning points are near the radiative-convective interface r c .
The mode displacements ξ r steepen towards the stellar center due to geometric focusing and conservation of wave flux (Goodman & Dickson 1998;Barker & Ogilvie 2010).Their maximum displacement peaks at r inner and scales as ξ r,max ∝ r −2 inner .As a result, nearly all the nonlinear mode coupling takes place near the inner turning points of the modes (see Appendix A in WAQB).
If the initial amplitude of the linearly resonant g-mode is sufficiently large, it can become strongly nonlinear (k r ξ r 1, where k r is the radial wavenumber) before reaching r inner .This will cause it to overturn the stratification and break rather than reflect.It should then be treated as a traveling wave rather than a global standing wave.For a given stellar model, the critical initial amplitude ξ r for wave breaking depends linearly on planet mass and only weakly on orbital period (P 1/6 orb ).Barker (2020) calculates the critical planet mass M p,crit above which the g-mode will break for a range of stellar models.By comparing with his Figure 9, we see that all of our main sequence models have M p,crit > 3.0M J and thus their g-modes are comfortably in the weakly nonlinear regime (the closest is the solar model which has M p,crit ≈ 3.3M J ; the rest have much larger M p,crit ).
The value of M p,crit decreases dramatically as a star evolves off the main sequence (Barker 2020).This is because as the star evolves its core contracts and the core's nearly constant gradient of N Cr subsequently increases.Since r inner of the resonantly excited g-mode is located where 2Ω N , i.e., r inner 2Ω/C, as the core contracts and C increases, the mode propagates to smaller r inner .As a result, k r ξ r,max increases, all else being equal.
This effect can be important for our models with M ≥ 1.2M , which, as explained in Section 3.1, are all on the subgiant branch.There are two models for each mass in this range, a less evolved and a slightly more evolved subgiant star.In the latter, the linearly resonant g-modes are close to the wave breaking limit for typical hot Jupiter masses (Sun et al. 2018;Barker 2020).These stars thus represent the transition between the weakly and strongly nonlinear regimes.We do not consider even more evolved subgiant stars or stars on the red giant branch since the efficiency of their tidal dissipation is instead found by analyzing the dynamical tide in the traveling wave regime (see, e.g., Barker & Ogilvie 2010;Barker 2011;Weinberg et al. 2017;Sun et al. 2018;Barker 2020).

Values of the coefficients in the amplitude equation
Since the g-modes we consider are all high-order (n a 30), their properties are well-approximated by asymptotic WKB expressions (Aerts et al. 2010).For example, their eigenfrequencies are approximately given by where Λ a = l a (l a + 1) and We illustrate the accuracy of this scaling in the top panel of Figure 2 for three of the models: (M /M , Age/Gyr) = (0.5, 4.97), (1.0, 4.60), and (1.5, 2.03).In Table 1, we give the value of ∆P for all our models.
As we now show, the various coefficients in Equation ( 8) are also all well-fit by power-law expressions.We use these expressions to build our mode networks since they enable a much faster search for mode triplets with low parametric threshold energies E th (see Section 4).(M/M , Age/Gyr) = (0.5, 4.97) in red, (1.0, 4.60) in green, and (1.5, 2.03) in cyan.The solid black lines show the fitted expressions given by Equations ( 17), ( 20), ( 24), and (25).

Linear damping rate γa
The dominant linear damping mechanism of the highorder g-modes in all our models is radiative damping, similar to the solar models shown previously (Terquem et al. 1998;Goodman & Dickson 1998).We use GYRE's solution of the non-adiabatic oscillation equations to determine the radiative damping rate γ (rad) a . We also calculate the damping due to turbulent viscosity in convection zones, which we find by computing where F (r) is given by the expression derived in Higgins & Kopal (1968) and depends on ξ a (x).The turbulent effective viscosity ν turb is a function of the ratio of the convective turnover frequency (provided by MESA) to the mode frequency and decreases as this ratio decreases.We calculate ν turb using the power-law expression given in Duguid et al. (2020) from a fit to their numerical simulations.We find that for the dynamically relevant modes, turbulent damping is always much smaller than radiative damping and γ a γ (rad) a .As in the case of a solar model (WAQB), the linear damping rate in all our models is well-fit by the expression This dependence stems from the fact that the radiative diffusion is proportional to the second derivative of the temperature fluctuation, and thus the damping rate of short-wavelength perturbations scales as the square of the radial wavenumber k 2 a ∝ (Λ a /ω a ) 2 .We illustrate the accuracy of this scaling in Figure 2 for three of the models (second panel from the top) and in Table 1 we give the value of γ 0 for all our models.We find that γ 0 increases dramatically with increasing M and stellar age.
Damping attenuates a mode's amplitude by exp(−α) during a round-trip between the turning points, where and travel time of a mode with group velocity v group = ∂ω a /∂k a .The second equality in Equation ( 21) follows from the WKB dispersion relation for g-modes (Equation ( 17)) and the scaling relation for γ a (Equation ( 20)).
If α 1, a significant fraction of the mode's energy is lost during a round trip and the mode is effectively a traveling wave (see, e.g., Goodman & Dickson 1998;Burkart et al. 2013;Sun et al. 2018).In that case, our treatment of nonlinear mode coupling, which assumes the modes are global standing waves, is no longer valid.By Equation ( 21), we see that α > 1 if the mode period P a = 2π/ω a exceeds a critical value In Table 1 we give the values of P crit for each of our models.For all the M ≤ 1.0M models and for the early subgiant models with M ≥ 1.2M , the l a = 2 values are in the range P crit [3,40] days.Thus, for these models, the linearly resonant parents (P a P orb /2) are comfortably in the standing wave regime.For the older subgiant models with M ≥ 1.2M , this is only true for hot Jupiter systems with P orb 1.5 days; at longer periods, the parent modes are in the traveling wave regime.However, we will see that even if we continue to treat them as standing waves, the mode network calculation yields a tidal dissipation rate that is close to that found with a traveling wave treatment.
The daughter and granddaughter modes in our networks have periods that are up to a few times larger than the parent modes.Moreover, they can have l a > 2. These modes are therefore more likely to be in the traveling wave regime, especially for the older and larger M models.For such modes, a proper treatment of their nonlinear interactions may require a hybrid formulation in which the parent mode is a standing wave and the secondary modes are traveling waves.Such a formalism has not been developed, as far as we know, and is left to future work.We suspect that since the nonlinear mode coupling happens very close to the stellar center (near the inner turning point of the modes) whereas the bulk of the linear dissipation occurs at considerably larger radii, the tidal dissipation rate due to the nonlinear excitation of secondary traveling waves may not be very different from that found here.

Linear tidal coefficient Ua
By plugging Equation (3) into Equation (7a), we can express the linear tidal coefficient U a in terms of the dimensionless linear overlap integral We evaluate I a using the method given by Equation (B23) in Burkart et al. (2013), which offers good numerical stability even for high-order modes.As in the solar model (see Figure 11 of WAQB), we find that the l a = 2 overlap integrals of our models are generally well-fit by We illustrate the accuracy of this scaling in Figure 2 for three of the models (third panel from the top) and in Table 1 we give the value of I 0 for all our models.We find that I 0 decreases with increasing M and stellar age (except for the most evolved and massive models, where it increases with age up the subgiant branch).For our high M subgiant models, I a can be larger than the scaling at short periods (see, e.g., the 1.5M model in Figure 2).This is because in those models, the mode wavelength becomes larger than the size of the convective envelope at short mode periods and the ω 11/6 a scaling only applies in the opposite limit.For those cases, we use the numerically computed I a rather than the scaling relation.
In our mode networks, we assume that only linearly resonant modes (parents) have non-zero linear tidal forcing U a .This is justified because U a is much smaller for the daughter modes, granddaughter modes, etc., and their driving is far off resonance.It therefore has a negligible secular effect compared to the resonant threemode interactions.Ignoring such forcing allows us to adopt the convenient change of coordinates described in EW16 (see Section 4.4 below) and significantly speeds up the integration of the amplitude equations.

Three-mode coupling coefficient
We calculate the three-mode coupling coefficient κ abc , which is symmetric under the interchange of mode indices, using Equations A55 through A62 in WAQB.Angular momentum conservation leads to the following angular selection rules for the three modes: (i) l a + l b + l c must be even, (ii) m a + m b + m c = 0, and (iii) the triangle inequality, |l a − l b | ≤ l c ≤ l a + l b .We focus on the parametric instability involving three-mode interactions between a high-order parent g-mode and a pair of high-order daughter g-modes whose summed frequency nearly equals the parent's frequency (and similarly, a daughter can couple to pairs of granddaughters, etc.).For such a triplet, the coupling is strongest in the stellar core, where the displacements of the modes peak (Section 3.2).As in the solar model (Appendix A in WAQB), we find that the coupling coefficients of our models are well-fit by where P a is the period of the parent mode and T abc ≈ 0.1 − 1 is an angular integral that depends on each mode's l and m and is easily evaluated in terms of Wigner 3-j symbols.The coupling occurs mostly near the parent's inner turning point r a,inner and scales as P 2 a because the parent's displacement there varies as ξ r,a ∼ r −2 a,inner ∼ P 2 a .
We illustrate the accuracy of this scaling in  12 in WAQB).In Table 1, we give the value of κ 0 for all our models.We find that it increases somewhat with increasing M and stellar age.

Building and Integrating the Mode Networks
As in the solar model calculations of EW16, we find that the parent can drive many parametrically unstable daughters to large amplitudes.These daughters can then drive parametrically unstable granddaughters to large amplitudes, and so on.The total number of potentially unstable modes and the number of couplings is larger than the number we can realistically simulate, especially in a survey that explores the parameter space of M , stellar age, P orb , and M p .Fortunately, as EW16 found in the case of a solar-type host and as we illustrate in Section 5 for non-solar type hosts, the total tidal dissipation rate can be reliably computed with a mode network that contains only a subset of the potentially unstable modes.This requires adopting EW16's approach to constructing mode networks, which we now summarize.

Selecting daughter modes
The approach in EW16 relies on the fact that despite the many daughter pairs with E thr < E lin (see Equations (11) and ( 10)), those with the lowest E thr overwhelmingly dominate the dynamics in large multi-mode systems.By adding pairs with progressively higher E thr to their networks, EW16 found that the system converges to a dissipation rate Ė that does not change significantly as even more modes are added.This requires including a sufficient number of generations (at least parents, daughters, and granddaughters) in order to obtain convergent results.Following EW16, we thus build our mode networks by systematically searching the mode parameter space and constructing, for each generation, a complete list of pairs ranked by E thr .
In order to carry out our search, we use the scaling relations for ω a , γ a , and κ abc (Equations ( 17 Angular degree l and radial order n for each mode in our P orb = 2.0 day networks for models (M /M , Age/Gyr) = (0.5, 4.97) in red, (1.0, 4.60) in green, and (1.5, 2.03) in cyan (the same models and colors as Figure 2).The networks consist of one parent mode (solid circles), ten pairs of daughter modes (crosses), and fifty pairs of granddaughter modes per daughter (small triangles).
angular integral T abc ).There is a tradeoff between finding daughters with smaller ∆ abc , which favors higher n and thus higher l for a given ω, and daughters with smaller γ, which favors smaller l since γ ∼ l 2 .The regions of small E thr therefore tend to occur where these two countering effects are balanced.After finding the local minima, we expand our search around the minima and find pairs with progressively higher E thr .Because γ ∼ l 2 , at high enough l the damping dominates detuning, and E thr increases with increasing l.We truncate our search upon reaching an l max such that E thr > E lin (i.e., a stable triplet).

Selecting granddaughter modes
The search for the lowest E thr daughtergranddaughter triplets is carried out in a similar way.Note that their E thr is generally smaller than that of the most unstable parent-daughter triplets.This is because the mode frequencies decrease with each generation and thus κ abc ∝ ω −2 a is larger (Equation ( 25)) while ∆ abc /ω a ∝ ω 2 a is smaller (Wu & Goldreich 2001).For a fuller discussion of this point, see Appendix F of EW16.Physically, the minimum E thr tends to decrease with each generation for two reasons.First, lower frequency modes penetrate deeper into the core.Their peak displacement ξ r,max are therefore greater (for a given mode amplitude) and thus their κ abc is larger.Second, such modes are more densely spaced in fre-quency and therefore can have smaller detunings.As a result, each generation is more susceptible than its antecedents to three-mode instabilities.

Reference networks
We find that for parent-daughter (daughtergranddaughter) coupling, the dissipation Ė is overwhelmingly dominated by the ≈ 10 (≈ 50) lowest E thr pairs per parent (daughter), similar to what EW16 found.Thus, as in EW16, our reference networks consist of one parent mode, its 10 lowest E thr daughter pairs, and the 50 lowest E thr granddaughter pairs per daughter.For most networks, this corresponds to 20 daughters and 2000 granddaughters.The exception is the more evolved subgiant models with M ≥ 1.2M and systems with larger P orb .In those cases, the total number of granddaughters in the network can be several times smaller since many of them couple to more than one daughter (the number of pairs is still 50).This is because in those models the modes have much higher damping rates (see γ 0 values in Table 1), which strongly favors low-l modes and narrows the pool of distinct granddaughters pairs with low E thr .
When integrating our networks, we choose the orbital period such that the parent is located midway between linear resonances, i.e., at a resonance trough.Thus, we might set the orbital period to, say, P orb = 1.01 day rather than P orb = 1.00 day.We also carry out a restricted set of runs in which the parent is located at a resonance peak.We will see that for some cases, τ is nearly the same regardless of whether the parent is at a resonance trough or peak; for example, this is true of a solar model, as also found by EW16 (see their Figure 2).In other cases, τ is significantly smaller for a parent at a resonance peak.Such systems will move quickly through resonances and, as a result, the long term decay rate will be close to the off-resonance τ value.
In Figure 3 we show the network structure for three of our stellar models at P orb = 2.0 day; our other networks have qualitatively similar structure to these.We see that the daughters and granddaughters with the smallest E th are typically l 10 modes.However, for the models with larger ∆P (i.e., those with smaller M and age), their l values can sometimes be even larger.

Integration method
We integrate each mode's coupled amplitude equation (Equation ( 8)) using a method similar to EW16's (see their Section 3.3 for a discussion).In particular, we change coordinates to x a = q a e iωat and assume that only the linearly resonant parent modes have a non-zero linear tidal forcing U a .This is a good approximation since the daughters and granddaughters are far from being linearly resonant and they also have much smaller U a ∝ ω 11/6 a (Equation ( 24)).This allows us to greatly speed up the integrations by canceling out the relatively high frequency ω a term in the amplitude equation and replacing it with a much lower frequency term proportional to ∆ abc .
We use the CVODES Adams solver from the SUNDIALS package (Hindmarsh et al. 2005) and parallelize the computations across multiple CPUs using MPI.Each simulation typically takes one to two days to run on a node of a local cluster.
A given simulation run consists of a particular stellar model, orbital period, and planetary mass.For each of our sixteen stellar models (see Table 1), we consider P orb = {0.5, 1, 2, 3, 4} days and M p = {0.5, 1, 2, 3} M J for a total of 300 runs.

Results
In this section we present the results of integrating the amplitude equation for large sets of nonlinearly coupled modes (Equation ( 8)) using the models and mode parameters presented in Section 3 and the mode networks built according to the procedure described in Section 4. In Section 5.1, we show examples of the time-dependent mode dynamics for two representative networks.In Sections 5.2 and 5.3 we show the full set of results, expressed in terms of the average orbital decay time τ , the tidal quality factor Q , and the transit time shift T shift (Equations ( 14), (15), and (16), respectively) for the runs from the lower mass stellar models (0.5 ≤ M /M ≤ 1.0) and the higher mass stellar models (1.2 ≤ M /M ≤ 2.0), respectively.Table 4 in Appendix B also lists the values of τ , Q , and T shift from each of the runs.

Example mode dynamics
In Figure 4 we show the time-dependent energy E a (t) and energy dissipation rate Ėa (t) of oscillation modes excited in a star due to its tidal interaction with a Jupiter-mass planet on a one day orbit.The left panels show results for a lower mass main sequence star and the right panels show results for a higher mass early subgiant star.In both cases, we see in the top panels that the energy of the linearly resonant parent mode (black line) is similar to the energies of individual daughter modes (red points) and granddaughter modes (blue points).This is consistent with the notion that upon reaching nonlinear equilibrium, a large system of coupled modes will be in approximate energy equipartition.
The solid red and blue lines show the summed energy a E a (t) of all the daughters and granddaughters, respectively.We see that the summed energy of all the modes (purple line) is dominated by the contribution of the granddaughters and, to a slightly lesser extent, the daughters and that their contribution is ∼ 10 − 100 times larger than the parent's (i.e., while the energy of individual modes tend to be similar across generations, there are many daughters and granddaughters and their summed contribution dominates over the parent's).Thus, at equilibrium, the energy being pumped into the parent by linear tidal driving is efficiently transferred to the secondary modes.
The bottom panels of Figure 4 show that unlike the energy E a (t), the energy dissipation rate Ėa (t) of individual daughter and granddaughter modes is often much larger than that of the parent.This is because their wavelengths are often much shorter.They thus have higher linear damping rates γ a and at the same energy are much more dissipative.The solid lines show that the total energy dissipation rate is dominated by the contribution of the daughters and granddaughters and is ∼ 10 3 − 10 4 times larger than the parent's contribution.We also see that the system reaches a nonlinear equilibrium within ∼ 10 5 P orb and, unlike the energy in individual modes, the total dissipation has a fairly small relative standard deviation (1.2 and 0.2 in the left and right panels, respectively).
As described in Section 4, our reference network consists of one parent mode, its 10 lowest E thr daughter pairs, and the 50 lowest E thr granddaughter pairs per daughter; this typically corresponds to ≈ 2000 total modes.EW16 found that such a network yields total Ė comparable to (to within a factor of ≈ 2) much larger networks that contain more modes and generations.We carried out a few numerical experiments with larger networks and likewise find that they yield Ė values comparable to that of our reference network.We illustrate this with the green line in the bottom right panel of Figure 4, which shows the total Ė for the same planetary and stellar parameters as the purple line but for a network that consists of one parent mode, its 10 lowest E thr daughter pairs, the 50 lowest E thr granddaughter pairs per daughter, and the 100 lowest E thr great-granddaughter pairs per granddaughter.The network contains about 2.3 × 10 4 modes and is thus 10 times larger than the reference network (and comparable in size to the largest networks EW16 ran; see their Figure 3).Since such a large network runs much slower, we stopped the integration after about 2 × 10 5 P orb (the wall time was about ten days).As the green line shows, this is long enough to see that Ė is within a factor of 2 of the reference network.Similar experiments with a few other networks yield similar results, suggesting  The blue points show the ≈ 2000 granddaughter modes and the blue lines show their cumulative contribution; for clarity, we only plot individual points for one out of every 20 granddaughter modes (many of them also lie below the range of the plotted vertical scale).The purple lines show the total energy and energy dissipation rate found by summing over all the modes in the network.The green line in the bottom right panel shows the total energy dissipation rate for a network that includes 100 great-granddaughter pairs per granddaughter.
that, like in EW16, the reference network is sufficiently large to yield convergent results.
5.2.Lower mass stars: 0.5 ≤ M /M ≤ 1.0 In Figure 5, we show the orbital decay time τ (solid circles) as a function of P orb and M p for the runs with the M = {0.5, 0.8, 1.0}M stellar models.We also show the corresponding stellar tidal quality factor Q (open squares) for the M p = 1.0MJ case; in Table 4 of Appendix B, we list the values for all the cases.When computing the time-average, we neglect the early portion of a run (typically the first 10 6 orbits) to allow transients from the initial conditions to die away and to ensure the system has reached an equilibrium.In Figure 6, we show the values of T shift for these same models assuming P orb = 1.0 day and M p = 1.0MJ or M p = 3.0M J (the full set of values are given in Table 4).
As expected, the value of τ is most sensitive to P orb ; in many cases it has a nearly power-law dependence with τ ∝ P 6.5 orb , approximately.It is also sensitive to the stellar model, with τ decreasing as M and stellar age Figure 5. Orbital decay time τ (solid circles) as a function of orbital period P orb for the M = {0.5, 0.8, 1.0}M models at different ages, as indicated by the label in each panel.The green, blue, red, and black solid circles show τ for a parent that is off-resonance (i.e., between resonance peaks) and Mp = {0.5, 1.0, 2.0, 3.0}MJ, respectively.The cyan squares show the corresponding stellar tidal quality factor Q for Mp = 1.0MJ (the full set of results is given in Table 4).The open blue circles show τ for a parent on a resonance peak for the case Mp = 1.0 MJ, P orb 1.0 day; in four of the panels, the open blue circle blends in with the solid circles since the off-resonance and on-resonance τ are similar.The blue triangles show τ for a parent treated as a traveling wave rather than a standing wave for the case Mp = 1.0 MJ.In the upper left panel, we subtract 5 dex and 3 dex from the off-resonance τ and Q , respectively, so that they appear on the plotted scale.In the bottom right panel, the x's show τ for the 11.0 Gyr model using the same color scheme as the solid circles.
increase.For all of the M = 0.8M and 1.0M models -from the early main sequence to the late main sequence or early subgiant phase -we find τ 10 8 yr for M p ≥ 1.0M J and P orb ≤ 1.0 day.The corresponding quality factors and transit time shifts are Q ≈ 10 5 −10 6 and T shift ≈ 10 − 100 s.For the M = 0.5M model, τ is negligible on the early main sequence (this is largely because κ 0 is signficantly smaller for this model; see Table 1); however, for the older M = 0.5M models at P orb 1.0 day and M p 1.0M J , we find that τ is less than the age of the system, Q 10 7 , and T shift 10 s.The value of τ tends to be fairly insensitive to planet mass as long as M p 1.0M J ; at smaller planet mass, the nonlinearities start to become less significant and τ can become very large, especially for the younger, lower M models.
Our standard calculation assumes the parent is between linear resonance peaks (see Section 4).To examine how sensitive τ is to the parent's degree of resonance, for each stellar model we also did a run with a  4.
parent located exactly on a resonance peak.The open blue circles in Figure 5 show the result for the particular case M p = 1.0MJ and P orb = 1.0 day.For the three M = 0.5M models and for the young (1.0 Gyr) M = 0.8M and 1.0M models, τ is significantly smaller for a parent on a resonance peak.The system will therefore evolve relatively quickly through the resonance peak and ultimately spend most of its time between peaks.For the other M = 0.8M and 1.0M models, τ is nearly the same on a resonance peak as it is between peaks (consistent with what EW16 found for the solar model).Whether on a resonance peak or between, we conclude that the solid circles should be a reasonably accurate estimate of τ over long timescales.A phenomenological model for estimating nonlinear tidal dissipation that can account for differences between onand off-resonance driving was presented by Yu et al. (2020) in the context of white dwarf binaries.In future work, we plan to investigate whether this model can be applied to hot Jupiter systems.The blue triangles in Figure 5 show τ for a parent treated as a traveling wave rather than a standing wave assuming M p = 1.0MJ .The results are found by interpolating across the Q IGW values given in Figure 8 of Barker (2020), applying the P 8/3 orb scaling (assuming the tidal forcing frequency is twice the orbital frequency), and using Equation ( 15) to convert the results to a decay time τ .The parent would be a traveling wave if its amplitude is large enough to undergo wave breaking or if its linear damping is high enough.Although most hot Jupiters are not in this regime (see Sections 3.2 and 3.3.1),it is nonetheless interesting to compare the τ predictions for the two regimes.We see that for an off-resonance parent, the traveling wave τ is generally smaller then the standing wave τ , although for the higher M models at short P orb they can be similar.By contrast, for an on-resonance parent (see the open blue circles in Figure 5), the traveling wave τ can sometimes be larger then the standing wave τ . 4In such cases, this suggests that the parent will be a traveling wave while the system is passing through a resonance and that the correct τ during those times is given by the traveling wave value.5.3.Higher mass stars: In Figure 7 we show τ and Q for the higher mass stellar models: M = {1.2,1.5, 2.0}M .The corresponding T shift values are shown in Figure 6 and Table 4.These models are all in the subgiant phase since, as explained in Section 3.1, main sequence stars with M 1.2M have convective cores and nonlinear coupling of g-modes is negligible in their interior.
For each M , we find that τ decreases significantly as the star evolves from the early to later subgiant phase.We find that at P orb 1.0 day, the dissipation can be especially efficient in the later subgiant models, with τ 10 5 yr, Q 10 5 and T shift 10 4 (T dur /10 yr) 2 s.As Figure 6 and Table 4 show, the T shift of these models are orders of magnitude larger than that of the M ≤ 1.0M main sequence models.Since T shift ≈ 100 s should be detectable, such systems should produce a detectable transit time shift after only T dur ≈ 1 yr.
The orbital decay is much faster for more evolved subgiants because they have significantly larger γ 0 and I 0 , and to a lesser extent, larger κ 0 (see Table 1).The damping rate is larger because as the stars evolve up the subgiant branch, the core contracts and C = N/r increases.Since k r ∝ C/ω, at a given P orb the resonant modes will have shorter wavelengths and thus larger damping rates.Even at P orb 4 day, the longest orbital period we consider, τ ≈ 10 8 yr for the older subgiant models.This is comparable to the subgiant evolutionary timescale, and thus at these larger P orb the changing stellar structure can dictate the planet's decay  Similar to Figure 5 but for the M = {1.2,1.5, 2.0}M models at different ages, as indicated by the label in each panel.The green, blue, red, and black solid circles show τ for an off resonance parent (i.e., between resonance peaks) and Mp = {0.5, 1.0, 2.0, 3.0}MJ, respectively.The cyan squares show the corresponding stellar tidal quality factor Q for Mp = 1.0MJ (the full set of results is given in Table 4).The blue triangles in the 1.2M and 1.5M panels show τ for a parent treated as a traveling wave rather than a standing wave for the case Mp = 1.0 MJ. trajectory (see also Sun et al. 2018).Interestingly, in those models τ is nearly independent of M p .Note too that τ is also close to the traveling wave value (blue triangles); indeed, for these older stellar models, the critical planet mass for wave breaking is around 1.0M J (Sun et al. 2018;Barker 2020) and thus the parent is close to, if not in, the traveling wave regime.

Implications for Known Hot Jupiters
We now consider the implications of our orbital decay calculations for known hot Jupiter systems.In Section 6.1, we discuss the two systems, WASP-12 and Kepler-1658, with observational evidence of tide-induced orbital decay and place them in the context of our results.In Section 6.2, we discuss the systems that our calculations suggest are especially likely to be undergoing rapid orbital decay and we recommend be included in campaigns that search for tide-induced transit timing variations.
The orbital decay of WASP-12b can be understood if the host star is a subgiant with M 1.2M ; the linearly excited primary mode (the parent) is then in the wave breaking regime and the expected Ṗorb agrees well with the observed value (Weinberg et al. 2017).However, Bailey & Goodman (2019) find that subgiant models are in greater tension with current observational constraints than main-sequence models with M 1.3 − 1.4M .Since such a main-sequence star has a convective core, the excited primary mode should be a standing wave, i.e., a g-mode, and stable to the parametric instability.It is then not clear how to explain the observed decay since the linear and nonlinear dissipation of the g-mode is too small.
Kepler-1658 has a mass M 1.5M and is more definitively a subgiant star.Vissapragada et al. (2022) explain that the measured decay is in good agreement with theoretical predictions for inertial wave dissipation (Barker 2020).Such dissipation is possible because the tidal frequency 2(Ω − Ω spin ) is less than twice the star's spin frequency Ω spin (this is because of the relatively large P orb = 3.85 days and the star's relatively short spin period P spin = 2π/Ω spin = 5.66 days; WASP-12b orbits too quickly and the host star spins too slowly for the tide to excite inertial waves).Although the linearly excited gravity wave breaks as it approaches the stellar center, we find that the wave luminosity is too small to explain the observed orbital decay (here we are ignoring rotation; the wave luminosity scales approximately as   2016

Systems likely to be undergoing rapid orbital decay
In Tables 2 and 3 we list hot Jupiter systems that our calculations suggest are favorable targets for the search for orbital decay, separated into systems with M ≤ 1.1M and M > 1.1M , respectively.Given the strong dependence of τ on P orb , these tables essentially list the systems with especially small P orb .However, we exclude short P orb systems whose host stars are likely to have convective cores (main sequence stars with M 1.2M ) since we find that such stars are unlikely to be undergoing significant orbital decay; examples of such systems include WASP-18 (P orb = 0.94 days, M = 1.3M , Age = 1.6 +1.4 −0.9 Gyr; Cortés-Zuleta et al. 2020), HATS-52 (P orb = 1.37 days, M = 1.1M , Age = 1.2 +1.5 −1.1 Gyr; Henning et al. 2018), and CoRoT-1 (P orb = 1.51 days, M = 1.2M , Age = 1.6 ± 0.5 Gyr; Bonomo et al. 2017).That said, the systems listed in Table 3 have sufficiently large age uncertainties that many are likely still on the main sequence, in which case their τ would be large.
For the systems listed in Table 2 with P orb 1.0 day, our results suggest τ 100 Myr and Q 10 6 (given the uncertainties in the measured parameters, especially the stellar age, we do not attempt to make system specific predictions).If a system listed in Table 3 has a host star on the subgiant branch, then it too should have τ 100 Myr and Q 10 6 .By Equation ( 16), we expect systems with such small τ to have T shift 10 s after T dur = 10 yr.Thus, the transit timing shifts of some of these systems, many of which are part of existing transit timing measurement campaigns (e.g., Patra et al. 2020;Ivshina & Winn 2022;Maciejewski et al. 2022;Mannaday et al. 2022;Rosário et al. 2022;Harre et al. 2023), may soon be detectable.

Summary and Conclusions
We studied the dissipation of the dynamical tide due to the excitation and damping of weakly nonlinear gmodes in hot Jupiter host stars.This is the dominant source of tidal dissipation in the great majority of hot Jupiter systems.By integrating the amplitude equations for large networks of nonlinearly interacting g-modes consisting of a linearly driven parent exciting a sea of secondary modes (daughters and granddaughters), we calculated the tidal dissipation rate as a function of stellar mass M , age, orbital period P orb , and planetary mass M p .We expressed our results in terms of the orbital decay time τ , stellar tidal quality factor Q , and transit time shift T shift .In order to span the range of observed hot Jupiter systems, our analysis considered 16 stellar models with masses between M = 0.5 − 2.0M and ages between the early main sequence and the subgiant branch.For each stellar model we considered orbital periods between P orb = 0.5 − 4 days and planetary masses between M p = 0.5 − 3.0M J .Our main results are shown in Figures 5-7 and in Table 4.
For hosts with M 1.0M , we found that the tidal dissipation rate increases with M and age, and can be significant throughout the main sequence.The dissipation is sensitive to P orb and in many cases follows a power-law τ ∝ P α orb with α 6−7.For P orb 1.0 day and M p 0.5M J we in general found that τ 100 Myr, Q 10 6 , and T shift 10(T dur /10 yr) 2 s, where T dur is the duration of a system's observation time.
For hosts with M 1.2M , we found that the tidal dissipation rate is negligible on the main sequence but becomes highly significant as these stars ascend the subgiant branch (owing primarily to the respective presence and absence of a convective core).Even at P orb ≈ 3 days, the dissipation rate on the subgiant branch is rapid enough to produce a detectable transit time shift (T shift 10 s) within 10 years.For P orb ≈ 1.0 day, a detectable transit time shift takes only ≈ 1 yr, which is ∼ 100 times faster than for main sequence stars with M 1.0M (see Figure 6).We compared our results to known hot Jupiters and identified a number of systems that could be undergoing rapid orbital decay (see Tables 2 and 3).However, since our results are sensitive to the age of the star (es-pecially for M 1.2M ), which is usually the least well-measured parameter of the hosts, there can be considerable uncertainty in the expected τ and T shift of individual systems.
Our analysis assumed that the excited modes are all global standing waves.We found that this is a good approximation for the linearly excited parent, except perhaps in our most evolved subgiant models.In those models, the parent's radiative damping can exceed its group travel time through the star or its displacement amplitude ξ r can approach the wave breaking limit (k r ξ r 1) near the stellar center.If either is true, the parent should instead be treated as a traveling wave since it dissipates all of its energy and angular momentum before it can reflect near the stellar center and form a standing wave.Several studies have considered the tidal dissipation rate in hot Jupiter systems under the assumption that the parent is a traveling wave and thus maximally dissipative (Barker & Ogilvie 2010;Barker 2011;Chernov et al. 2017;Weinberg et al. 2017;Sun et al. 2018;Lazovik 2021).
Treating all the daughters and granddaughters as standing waves may not always be a good approximation, however, even for some of our main sequence models.This is because their wavelengths can be much shorter than the parent that excites them.They therefore have higher radiative damping rates and at a given mode energy are closer to the wave breaking limit than the parent.This suggests that the interaction between the parent and the secondary modes may in some cases involve a standing wave (the parent) nonlinearly exciting traveling waves (the daughters, granddaughters, etc.).The computational methods for studying weakly nonlinear interactions between standing waves and traveling waves have not been developed, as far as we know, and is left for future work.Since the secondary modes dissipate their energy by radiative damping far from the stellar center where they are nonlinearly excited, it is possible the total dissipation will be insensitive to whether they are treated as standing waves or traveling waves.
Another caveat of our calculation is that we do not account for possible changes to the stellar structure due to the transfer of energy and angular momentum from the modes to the background star.This type of interaction was investigated recently by Guo et al. (2023), who performed 2D hydrodynamical simulations of tidally excited nonlinear gravity waves in the cores of solar-type stars.They found that linear damping of the waves gradually spins up the core and that subsequent incoming (parent) waves are absorbed in an expanding critical layer.Importantly, this process was found to occur even when the parents are below the wave breaking threshold, suggesting that such parents are effectively in the traveling wave regime and thus maximally dissipative.Due to computational limitations, they say that the secondary waves generated by the parametric instability, which are the focus of our analysis, would be very difficult to observe in their simulations.Properly resolving such secondary waves might be important, however.For example, they could propagate outwards from where they are excited and transfer their energy and angular momentum at larger radii rather than locally.Their dissipation might also modify the stratification of the central region, moving the parent's inner turning point outwards and causing it to reflect at larger radii.Such effects can potentially impact the formation of a critical layer and the subsequent absorption of the incoming parent wave.This work was supported by NSF grant No. AST-2054353.Software: MESA (Paxton et al. 2011(Paxton et al. , 2013(Paxton et al. , 2015(Paxton et al. , 2018(Paxton et al. , 2019;;Jermyn et al. 2022, http://mesa.sourceforge.net), GYRE (Townsend & Teitler 2013;Townsend et al. 2018, https://bitbucket.org/rhdtownsend/gyre/wiki/Home), SUNDIAL (Hindmarsh et al. 2005, https: //computing.llnl.gov/projects/sundials).The stellar models are built using version 15140 of MESA (Paxton et al. 2011(Paxton et al. , 2013(Paxton et al. , 2015(Paxton et al. , 2018(Paxton et al. , 2019;;Jermyn et al. 2022).The key parameters of the inlist files we use are provided below.We change initial_mass to build different mass models and we change max_center_cell_dq, R_function_weight, and R_function_param if finer grid resolution is needed in the stellar center (e.g., to calculate κ abc ).     4. Orbital decay time τ , stellar tidal quality factor Q , and transit time shift T shift from the mode network integrations.Each column is for a particular stellar model (M and Age) and each triplet of rows is, from top to bottom, for an orbital period P orb = 0.5, 1.0, 2.0, 3.0, 4.0 days.For a given P orb , the first row gives log 10 ( τ /yr), the second row gives log 10 Q , and the third row gives log 10 (T shift /s) assuming T dur = 10 yr.The four comma-separated entries in each row are, from left to right, for planet mass Mp = 0.5, 1.0, 2.0, 3.0MJ.The top of the table illustrates the format.− 11.9, 11.5, 11.4, 11.4 8.8, 8.4, 8.4 − 7.6, 7.8, 8.0, 8.1 4.5, 4.6, 4.8, 4.9 − − − 1.9, 2.0, 2.1, 2.2 5.9, 6.1, 6.2, 6.3 1.0 − − 9.5, 9.2, 9.0, 8.9 5.7, 5.6, 5.5, 5.4 − 8.4, 8.4, 8.4, 8.5 5.5, 5.6, 5.8, 5 − 11.4, 11.4, 11.4, 11.4 6.4, 6.7, 7.

Figure 1 .
Figure 1.Distribution of hot Jupiters in the plane of stellar mass M and orbital period P orb for planets with mass Mp sin i > 0.5MJ , where i is the orbital inclination.The color scale shows the best-fit stellar ages in Gyrs.Data is from the NASA Exoplanet Archive and only stars with reported ages are shown.
Figure 2 for three of the models (bottom panel).In the figure, the triplets consist of l a = l b = l c = 2 modes and nearly resonant daughters with detuning |ω a +ω b +ω c | < |10 −3 ω a | and |n b − n c | < |0.8n a |.The magnitude of κ abc varies by factors of order unity depending on the particular daughter pair and is negligible if |n b − n c | |n a | (see Figure Figure 3.Angular degree l and radial order n for each mode in our P orb = 2.0 day networks for models (M /M , Age/Gyr) = (0.5, 4.97) in red, (1.0, 4.60) in green, and (1.5, 2.03) in cyan (the same models and colors as Figure2).The networks consist of one parent mode (solid circles), ten pairs of daughter modes (crosses), and fifty pairs of granddaughter modes per daughter (small triangles).

Figure 4 .
Figure 4. Mode energy (top panels) and energy dissipation rate (bottom panels) as a function of time for the planetary parameters (Mp, P orb ) = (1.0MJ,1.0 day).The left panels are for the stellar model (M , Age) = (0.8M , 4.87 Gyr) and the right panels are for (1.2M , 4.00 Gyr).The black lines show the parent mode.The red points show the twenty daughter modes and the red lines show their cumulative contribution ( Ea and Ėa, where the sums run over only the daughter modes).

Figure 6 .
Figure 6.Transit time shift T shift as a function of stellar mass M for P orb = 1.0 day with Mp = 1.0MJ (blue points) and Mp = 3.0MJ (red points) assuming T dur = 10 yr.The ages of the stellar models increase from bottom to top.In three cases, T shift is below the scale of the plot.The full set of T shift results are given in Table4.
Figure 7.Similar to Figure5but for the M = {1.2,1.5, 2.0}M models at different ages, as indicated by the label in each panel.The green, blue, red, and black solid circles show τ for an off resonance parent (i.e., between resonance peaks) and Mp = {0.5, 1.0, 2.0, 3.0}MJ, respectively.The cyan squares show the corresponding stellar tidal quality factor Q for Mp = 1.0MJ (the full set of results is given in Table4).The blue triangles in the 1.2M and 1.5M panels show τ for a parent treated as a traveling wave rather than a standing wave for the case Mp = 1.0 MJ.

Facility:
Exoplanet Archive This research has made use of the NASA Exoplanet Archive, which is operated by the California Institute of Technology, under contract with the National Aeronautics and Space Administration under the Exoplanet Exploration Program.APPENDIX A. PARAMETERS OF THE MESA INLISTS , and 7).

Table 1 .
Properties of the stellar models and their modes.

Table 2 .
Favorable targets for the search for tidal orbital decay for host stars with M ≤ 1.1M .The age column includes the uncertainty, which is symmetric about the best fit if enclosed in parentheses.The fractional uncertainties are generally 0.001% for P orb and 10% for M and Mp.References.(1) McCormac et al. (2020); (2) Bonomo et al.

Table 3 .
Favorable targets for the search for tidal orbital decay for host stars with M > 1.1M .We exclude WASP-12 and Kepler-1658 from this list since their orbital decay has been detected (see Section 6.1).The age column includes the uncertainty, which is symmetric about the best fit if enclosed in parentheses.The fractional uncertainties are generally 0.001% for P orb and 10% for M and Mp.
TABLE OF ORBITAL DECAY RESULTS FROM THE MODE NETWORK INTEGRATIONS In Table 4 we list the values of τ , Q , and T shift from each of our mode network integrations (see also Figures 5, 6