Tides in the high-eccentricity migration of hot Jupiters: Triggering diffusive growth by nonlinear mode interactions

High eccentricity migration is a possible formation channel for hot Jupiters. However, in order for it to be consistent with the observed population of planets, tides must circularize the orbits in less than $\approx$ a Myr. A potential mechanism for such rapid circularization is the diffusive growth of the tidally driven planetary f-mode. Such growth occurs if the f-mode's phase at pericenter varies chaotically from one pericenter passage to the next. Previous studies focused on the variation of the orbital period due to tidal back-reaction on the orbit as the source of chaos. Here we show that nonlinear mode interactions can also be an important source. Specifically, we show that nonlinear interactions between a parent f-mode and daughter f-/p-modes induce an energy-dependent shift in the oscillation frequency of the parent. This frequency shift varies randomly from orbit to orbit because the parent's energy varies. As a result, the parent's phase at pericenter varies randomly, which we find can trigger it to grow diffusively. We show that the phase shift induced by nonlinear mode interactions in fact dominates the shift induced by tidal back-reaction and significantly lowers the one-kick energy threshold for diffusive growth by about a factor of 5 compared to the linear theory's prediction. Nonlinear interactions could thus enhance the formation rate of hot Jupiters through the high-eccentricity migration channel and potentially mitigate the discrepancy between the observed and predicted occurrence rates for close-in gas giants as compared to those further from the star.


INTRODUCTION
More than 25 years after the first detection of a hot Jupiter (Mayor & Queloz 1995), we still do not know their dominant formation channel. Possibilities include in situ formation, gas disk migration, and higheccentricity tidal migration (see Dawson & Johnson 2018 for a review). In the latter scenario, the planet is born beyond the snow line at 1 AU and is driven to high eccentricity through planet-planet scattering (e.g., Rasio & Ford 1996;Chatterjee et al. 2008) or secular interactions with another planet or star (e.g., Wu & Mur-ray 2003;Nagasawa et al. 2008;Wu & Lithwick 2011;Hamers et al. 2017;Teyssandier et al. 2019). Strong tidal interactions during close pericenter passages subsequently damp the eccentricity and shrink the semi-major axis, culminating in a planet that resides in a days-long circular orbit.
An outstanding problem with this formation channel is the lack of very high eccentricity systems (e > 0.9) among the observed population of hot Jupiters. In order to sufficiently speed up the circularization and thereby reduce the likelihood of catching a planet in the high-e state, the tidal dissipation must be at least ten times more efficient than our own Solar System's Jupiter (Socrates et al. 2012;Dawson et al. 2015). Such an enhanced efficiency is not necessarily inconceivable, however, as tidal dissipation in hot Jupiter systems can be sensitive to the strength and frequency of the tidal forcing and the structure of the components (see, e.g., Ogilvie & Lin 2004;Barker 2011;Essick & Weinberg 2016). One should therefore consider the problem from first-principles rather than rely on parametrized extrapolations.
Indeed, Wu (2018) showed that the diffusive growth of the the planet's l = 2 f -mode during high-eccentricity migration could lead to very rapid orbit circularization (see also Vick & Lai 2018;Vick et al. 2019). This process was first considered in the hot Jupiter context by Ivanov & Papaloizou (2004) and has also been considered in a number of other high eccentricity systems; e.g., in tidal capture binaries (Kochanek 1992;Mardling 1995) and eccentric neutron star binaries (Vick & Lai 2018). In Wu's calculations, the f -mode's phase is randomly perturbed by the backreaction of the tide on the orbit, causing the mode amplitude to grow diffusively over many pericenter passages. She argued that the f -mode will damp nonlinearly when its amplitude reaches unity and the mode breaks near the planet's surface. Within ∼ 10 4 yr, the planet is transported from a few AU to ∼ 0.2 AU and its eccentricity is decreased from near unity to < 0.9. Such a rapid circularization is equivalent to a remarkably small tidal quality factor of Q ∼ 1, five orders of magnitude smaller than Jupiter's. Wu (2018) showed that an additional feature of the diffusive growth scenario is that a planet that is secularly perturbed to high eccentricity will dynamically decouple from its perturbers when its pericenter distance reaches ∼ 4 tidal radii. All migrating Jupiters will therefore park safely outside the zone of tidal disruption, where they are observed today. This may explain why hot Jupiters are formed more efficiently than previous investigations of high-eccentricity tidal migration found: whereas observations show that the observed ratio of hot to cold Jupiters is ∼ 10%, previous theoretical calculations, which did not consider diffusive growth, yielded a ratio of only ∼ 1% due to their comparatively high rates of tidal disruption (see Table 2 in Dawson & Johnson 2018).
In this paper, we extend the work of Wu (2018) and Vick & Lai (2018) by considering the effects of weakly nonlinear mode interactions on the f -mode's diffusive growth. We show that the random changes in mode phase induced by three-wave nonlinear interactions act in concert with tidal backreaction on the orbit in order to lower the threshold for diffusive growth. We find that for a given orbital period, the diffusive growth can be triggered at a larger pericenter distance (i.e., smaller eccentricity) and hence smaller kick amplitude, as com-pared to calculations including only linear physics. Here we focus on the triggering and early phases of diffusive growth; an investigation of the long term evolution will be left to subsequent papers.
The plan of the paper is as follows. Section 2 contains a review of the coupled equations for the mode and orbit including only linear processes. Following Vick & Lai (2018), an iterative map for the mode amplitude and orbital period is outlined, including the effects of planetary rotation. Nonlinear coupling of the excited fmodes with other f-and p-modes is discussed in Section 3. The non-resonant phase shift and damping rate are derived, and gas giant planet models are used to evaluate the frequencies, damping rates and coupling coefficients. The nonlinear phase shift and damping are incorporated into the iterative mapping algorithm in Section 4. Results for short timescale simulations using the maps are presented in Section 5 and conclusions and discussion are presented in Section 6.

LINEAR PROBLEM
In this section, we review the linear problem and introduce some of the notation and approximations we will use throughout our study. In Sec. 2.1, we present the set of equations needed to construct an iterative map of the coupled mode-orbit evolution in linear theory (including the Doppler shifts caused by rotation). Then in Sec. 2.2, we justify our approximate treatment of angular momentum transfer and the orbital and spin evolution of the planet.
The eigenmodes form a complete basis for the fluid displacements ξ(x, t), which can be expanded as (Schenk et al. 2002) ξ(x, t) = a q a (t)ξ a (x, t) ξ(x, t) = a (−iω a )q a (t)ξ a (x, t), where ω a is the eigenfrequency of a mode and q a its amplitude. The sums run over both radial and angular quantum numbers as well as modes with positive and negative frequencies. We normalize each mode such that where E 0 = GM 2 /R, M is the mass of the planet, R is its radius, and other quantities have their usual meaning. If we ignore nonlinear effects, the equation of motion of a planetary mode in the frame corotating with the planet is (Schenk et al. 2002) where ω a and γ a are the eigenfrequency and linear damping rate of the mode. The amplitude of the tidal force acting on the mode is where D is the orbital separation, Φ is the orbital phase, Ω s is the spin of the planet, the tidal overlap Q a = (M R l ) −1 d 3 rρξ * · ∇(r l Y lm ), and at leading order (l = 2), the nonvanishing W lm coefficients are W 2±2 = 3π/10 and W 20 = − π/5. Note that we include the Doppler shift of frequency due to rotation but ignore corrections to the rotating-frame frequency and eigenfunction. Also note that our sign convention is different from that used in Wu (2018). Specifically, a prograde (retrograde) mode has m a > 0 (m a < 0) in our definition.
The mode amplitudes couple to the orbital motion through the accelerations a D and a φ in the equations of motionD where to linear order is the reduced mass. Throughout our study, we drop the nonlinear tidal back-reaction terms as their effect on the one-kick amplitude is subdominant.

Iterative map including only linear effects
The direct integration of the coupled mode-orbit evolution equations is computationally expensive. To obtain the approximate secular evolution, an iterative mapping procedure has been developed (see, e.g., Vick & Lai 2018), whose key steps we summarize below (see also similar derivations in Vick et al. 2019).
To do so, we first perform a phase shift to transform Eq. (3) from the corotating frame to the inertial framė where q a = q a exp (−im a Ω s t), ω a = ω a + m a Ω s t, and U a = U a exp (−im a Ω s t) are, respectively, the mode amplitude, mode frequency, and tidal driving in the inertial frame. The general solution of q a is q a (t) = e −(iω a +γa)t t iω a U a (τ )e (iω a +γa)τ dτ.
Suppose we know the mode amplitude right before the kth pericenter passage. We can then write the mode amplitude in the kth orbit as (see also Vick et al. 2019) where the superscript (0) and (1) indicate that the amplitudes are respectively evaluated right after and right before a pericenter passage. The quantity ∆q a,1 is the one-kick amplitude the mode receives at the pericenter. It is given by 1 In the equation above, the integration is preformed over one orbital period. Since we care about orbits that are highly eccentric, we make the approximation that the tidal interaction happens only near pericenter. Therefore, the limits of integration in Eq. (13) can be dropped as long as they bracket the pericenter passage. It is convenient to define an orbital integral K lm as 2 (Press & Teukolsky 1977) where D peri ≡ a orb (1 − e orb ) is the pericenter distance, e orb the eccentricity, and ω 0 ≡ GM/R 3 . If we ignore the perturbations on D and Φ, Lai (1997) provide an analytical expression for K 22 (i.e., l = m = 2) assuming the orbit is parabolic, where z ≡ √ 2ω/Ω peri and Ω 2 peri ≡ G(M + M * )/D 3 peri , and in the second line we have expanded the expression around z = 11 to emphasize the steep decline.
1 Formally the integration should be performed from right before the k'th pericenter passage to right before the next passage. We can nonetheless shift the initial time because U a (t) = U a (t + k P orb,k ) in the inertial frame if the pericenter distance stays approximately fixed throughout the evolution. This also is the reason the one-kick amplitude can be treated as a constant for different orbital cycles. 2 Note that a tidal field with spherical degree (l, m) linearly couples to a mode with la = l and ma = m due to the angular integral in Qa.
In terms of K lm , we can write the one-kick amplitude as The damping term entering ∆q a,1 can be safely dropped because γ a Ω peri . Note that for a parabolic orbit, K lm is a real number and therefore ∆q a,1 is purely imaginary. Also note that when calculating K lm one should use the mode frequency in the inertial frame ω a = (ω a + mΩ s ). Combining with the expansion given in Eq. (15), one sees immediately that the spin reduces the one-kick amplitude for a prograde mode with m > 0 in our convention.
To account for the tidal back reaction on the orbit, we adopt an energy conservation argument instead of explicitly coupling the mode amplitude equation and the tidal accelerations a D and a Φ . Upon receiving a kick at pericenter, the energy stored in a stellar mode changes by (including contributions from the mode and its complex conjugate) Since the energy stored in the tidal coupling (the term ∝ Re [q * a U a ]) is small everywhere except for at the pericenter and the spin rate of the planet should stay approximately fixed (which we will justify shortly), the change in the energy of stellar modes needs to be balanced by the orbital energy, where E orb,k = −GM M * /2a orb,k is the orbital energy at the k'th orbit. A direct consequence is that the change in the orbital energy also alters the orbital period P orb = 2π/Ω orb = 2π a 3 Since the value of ∆E a,k is different from orbit to orbit, the orbital period varies. This, in turn, leads to a stochastic evolution of the mode's phase per orbital cycle where we have used a subscript "br" to stand for the fact that this phase is due to the back reaction of the tide. As shown in previous studies (Vick & Lai 2018;Wu 2018), the randomness of the phase shift ∆φ br,k is key to triggering the diffusive growth of a tidally driven mode.
In order to simplify the notation, we will sometime omit the subscript "k" when we do not need the quantity to be evaluated at a specific orbit cycle.

Orbital and spin angular momentum
An energy transfer is typically associated with an angular momentum transfer as well. Nonetheless, since the change in the orbital angular momentum ∆L orb ∆E orb /Ω peri , we have that at high eccentricity (Vick & Lai 2018) Therefore, the orbital angular momentum stays as a constant throughout the evolution to a very good approximation. Suppose, for example, the initial orbit is a orb = 1 AU and e orb = 0.98 and it evolves to a orb = 0.2 AU (with e orb 0.9) due to tidal dissipation. Whereas the orbital energy changes by a factor of 5, the orbital angular momentum changes by only 3% in the process. Furthermore, since the angular momentum is nearly constant and at high-eccentricity, it follows that the pericenter distance D peri is also nearly constant throughout the orbital evolution. Moreover, under the high-eccentricity limit, the linear one-kick amplitude ∆q a,1 depends on the Keplerian elements only through the pericenter distance D peri [which determines Ω peri for fixed (M, M * )]. As the pericenter distance stays nearly unchanged, the one-kick amplitude ∆q a,1 also remains approximately constant.
So far we have left the spin of the planet Ω s as a free parameter. One plausible scenario is that the planet reaches pseudo-synchronization with the orbit via the equilibrium tide, leading to (Hut 1981) Ω s 1 + 15 2 e 2 orb + 45 8 e 4 orb + 5 16 e 6 orb (1 + 3e 2 orb + 3 8 e 4 orb )(1 − e 2 orb ) 3/2 Ω orb , 1.17Ω peri (e orb → 1), where the second line applies in the high-eccentricity limit. A constant pericenter distance (hence constant Ω peri ) would then imply that the spin frequency as set by the pseudo-synchronization condition also stays approximately constant. We note that including a pseudosynchronous rotation of the planet will increase the prograde f-mode frequency in the inertial frame, which will tend to decrease the one-kick amplitude and slow the orbital evolution as compared to the non-rotating planet case. Nevertheless, the one-kick amplitude of the prograde mode is still two orders of magnitude greater then the m a = 0 mode and even more for the retrograde mode. Therefore, we will only consider the prograde mode in the subsequent discussion.
To summarize our proceedure, we discard the evolution of the angular momenta and, self-consistently, treat Ω s and D peri as approximate constants during the circularization process (at least for the initial phase when 1 − e orb 1 is well satisfied). The evolution trajectories will thus reduce to the ones studied by Vick & Lai (2018) and Wu (2018) as long as one uses the mode frequency in the inertial-frame ω a = ω a + m a Ω s and neglects the nonlinear effects described in the next section.

NONLINEAR PROBLEM
We now consider how weakly nonlinear effects modify the problem. At lowest nonlinear order, the amplitude equation of a mode a is (Weinberg et al. 2012) where U ab is the nonlinear tide, κ abc is the three-mode coupling coefficient, and asterisks denote complex conjugation. There is a significant cancellation between the nonlinear tide and three-mode coupling to the equilibrium tide such that U ab + 2 c κ abc U c 0 (Weinberg et al. 2012). By treating the cancelation as perfect, we haveq We will solve this equation (or approximate its solution), in order to determine how nonlinear mode interactions influence the diffusive growth of the f -mode and thereby a planet's high-eccentricity tidal migration.
Since the daughters' direct, linear coupling to the tide is small, we expect the most significant nonlinear effect to be the modification of the parent mode's free evolution away from pericenter (when U a 0). Specifically, we show in Section 3.1 that the nonlinear mode couplings can be viewed as energy-dependent shifts of the parent's eigenfrequency and damping rate. We then derive the time-dependent evolution of such a nonlinear oscillator in Section 3.2. Lastly, we examine the nonlinear effects in a typical Jupiter model in Section 3.3.

Nonlinear frequency shift and effective damping rate
Previous studies have shown that at linear order, the l a = m a = 2 f-mode with ω a > 0 has the greatest energy and that it dominates the orbital evolution (see, e.g., Wu 2018). We will thus focus on a single parent mode (mode a) with l a = m a = 2 and ω a > 0 (and its complex conjugate to get real, physical quantities). We consider the nonlinear effects due to this parent mode coupling to itself and a daughter mode (which can be another fmode or a p-mode), which results in the inhomogeneous driving of the daughter. 3 Once we know the parent mode's angular pattern (l a , m a ), we can further utilize the three-mode angular selection rules and divide up the nonlinear couplings into two categories, which we will refer to as aab and aa * c, respectively.
In the aab case, the driving is formed by mode a coupling to itself. By the angular section rule, only daughter modes with l b = −m b = 4 can couple to this driving. We will refer to such a daughter as mode b and note that it is forced at a frequency −2ω a .
By contrast, in the aa * c case, mode a couples to its complex conjugate a * and drives a daughter mode (mode c) with m c = 0 and l c = 0, 2, 4. In this case, mode c experiences a forcing at zero frequency. 4 To make the abstract problem more transparent, we write out the explicit three-mode coupling equations for both cases. For simplicity, we start by considering only a single mode b and a single mode c. We will perform a summation over modes in the end to obtain the general solution. We will also drop the couplings involving more than one daughter mode for analytical simplicity; all the allowed couplings are included in the numerical calculations when we validate our analytical approximations.
3 As we consider fully convective Jupiter models, there are no lowfrequency g-modes that can parametrically couple to the parent (Weinberg et al. 2012). Therefore, the non-resonant nonlinearity considered in this work should be distinguished from the parametric instabilities considered in, e.g., Essick & Weinberg (2016) for solar-type stars and Yu et al. (2020) for white dwarfs. 4 Care must be taken when such DC forcing is encountered in the nonlinear problem in order to ensure that the forcing represents the physical transfer of energy and angular momentum between distinct oscillation modes, and not just a constant nonlinear modification of the linear mode frequency and eigenfunction. In the present situation, there is additional time-dependence due to the mode amplitudes, which change from one orbit to the next, and such forcing in turn leads to further time-dependent changes in the modes' amplitude and phase.
We then havė Note that the above set of equations describes the evolution away from pericenter, after the parent has received its most recent kick, since here we are interested in following the parent's free (i.e., unforced) evolution leading up to the next pericenter passage. Consequently, we do not include any tidal forcing terms.
To seek the leading-order nonlinear correction, we solve Eqs. (26)-(28) in a perturbative manner. Away from pericenter and without nonlinear couplings, we have q a ∼ exp(−iω a t). Using this as the driving term, the steady-state solutions of the daughters are 5 Plugging the daughter modes above back in to Eq. (26), we obtaiṅ where 6 andẼ a ≡ q * a q a is the dimensionless energy of mode a. 7 The physical mode energy including the contribution from both a and its complex conjugate a * is E a =Ẽ a E 0 in our normalization, where E 0 = GM 2 /R is the natural energy of the planet.
We thus see that the leading-order nonlinear correction corresponds to a shift in the eigenfrequency (conservative part) of the parent mode and an excess damping term (dissipative part), both of which depend linearly on the energy of the parent mode (see also Landau & Lifshitz 1976;Kumar et al. 1994;Kumar & Goodman 1996). We can therefore define where, after putting back the summation over all the daughters that couple to mode a, we have Note that mode c does not contribute to the nonlinear damping Γ. Mathematically, this can be understood by noticing that for every mode c with ω c , there exists a mode c * with −ω c (i.e., the complex conjugate of c; they both have m c = m c * = 0) that has the exact opposite contribution to Γ. Therefore, after summing over the ±ω c pair, the nonlinear dissipation due to mode c cancels exactly. Physically, we can view mode c as a nonlinear modification of the planet's structure, which changes the frequency at which the parent wave propagates (see footnote 4). Nonetheless, mode c is not a wave itself (as it is non-oscillatory) and therefore it does not contribute to the energy dissipation.
By contrast, mode b (corresponding to a nonlinearly excited wave oscillating at 2ω a ) enhances the dissipation, as one would expect physically. After the summation, δγ a is always positive because a mode b with negative (positive) frequency would contribute a positive (negative) dissipation rate, and (2ω a + |ω b |) > (2ω a − |ω b |) as the parent mode has ω a > 0. Consequently, we obtain a net increase in the damping after summing over each ±ω b pair. Now turn to the nonlinear frequency shift Ω. Its sign is not definite. While most of the modes 8 act to reduce the parent mode's frequency, a mode b with ω b < 0 and (2ω a +ω b ) > 0 will increase the parent mode's frequency. In practice, we find that only the l b = −m b = 4 f-mode satisfies the condition ω b (2ω a +ω b ) < 0; this mode can in fact be resonant with the parent, although we find that it is only a mild resonance since the mode spectrum is sparse for low-order (in both n and l) modes. Therefore, in general we would expect Ω < 0. However, in principle it could be positive if there is a rare strong resonance such that (2ω a + ω b ) 0 + . 9

Evolution of the nonlinear oscillator
In the previous section, we showed that nonlinear mode interactions perturb the frequency of the f -mode by δω a (Ẽ a ) and its damping rate by δγ a (Ẽ a ). These cause a dephasing of the f -mode δφ nl = − δω a dt in excess of the back reaction of the f -mode on the tide considered in previous studies of diffusive growth (e.g., Wu 2018). Furthermore, this change in phase varies from orbit to orbit due to the changes in the parent energy. Nonlinear effects can therefore contribute to, and even trigger (as we will show), diffusive growth. Since δω a depends on the energy of the modeẼ a , in order to determine δφ nl as a function of time, we need to determine howẼ a evolves. Note that here we focus on the evolution when the planet is far from pericenter, i.e., of the free oscillator. The goal of this section is therefore to determineẼ a (t) over an orbit, and from it calculate the nonlinear contributions to the dephasing δφ nl . The construction of an iterative mapping from orbit to orbit similar to that of the linear studies will be discussed later in Section 4.
The energy evolution is given by 10 If we substitute in Eq. (35) for δγ a (Ẽ a ), then the above equation can be solved easily as a is the initial mode energy and in the second line we expand exp [2γ a t] (1 + 2γ a t). This is a good approximation because the linear damping of the parent mode is typically small (1/γ a 10 9 yr for the Jupiter model we consider) and over the course of a ∼ 1 yr orbit, the condition γ a t 1 is very well satisfied. The total accumulated phase can be written as φ a (t) = − [ω a + δω a (Ẽ a )]dt. Of particular interest is the excess dephasing due to nonlinear interactions where we use the subscript "nl" to indicate the excess phase due to nonlinear effects (in contrast to tidal back reactions denoted by a subscript "br"), and in the second equality we change variables from time t to energyẼ a . If we use Eq. (34) for δω a and Eq. (38) for dẼ a /dt, then the nonlinear dephasing is We can also write the (leading-order) dephasing as a function of time by plugging in Eq. (39). In the limit that the parent mode's dissipation is small (γ a t → 0), we can cast the dephasing in an intuitive form as where in the second equality we further assumed 2Γ aẼ 10 −3 and t 1 yr, this condition is satisfied if Γ/ω a < 10 −5 . As we will see shortly in the following section and Table 1, Eq. (43) is well satisfied for the Jupiter model we consider in this work as it has a weak damping. Nonetheless, for Jupiters with greater radii, the damping rate can be significantly higher (Arras & Socrates 2009) and Eq. (42) should be used instead. Obviously, Eq. (43) also applies for conservative systems.

Values of Ω and Γ for a Jupiter model
From the discussion above, we see that the leadingorder nonlinear corrections to the parent mode correspond to shifts in both the eigenfrequency and the damping rate that are linearly proportional to the mode energyẼ a [Eqs. (34) and (35)]. The nonlinear dynamics can thus be characterized by the two coefficients Ω and Γ (both having the dimension of frequency in our definition). These coefficients further depend on the parent mode's eigenfrequency (with Ω, Γ ∝ ω a ) and the properties of the daughters. We now describe values for Ω and Γ for a Jupiter model with M = M J and R = 1.1 R J . Table 1. Properties of the Jupiter model considered in our study. We write E0=GM 2 /R=E0,43 × 10 43 erg, ω0= GM/R 3 =ω0,−4 × 10 −4 rad s −1 as well as Γ = Γ−10 × 10 −10 We use the subscript "prnt" to denote the parent mode (the prograde f-mode with la = 2) that directly couples to the tide.

R
E0,43 ω0,−4 ωprnt Qprnt Ω Γ−10 1.1RJ 3.1 5.1 1.1ω0 0.43 −29 ωprnt 8.1 ωprnt Using MESA (version 10398; Paxton et al. 2011Paxton et al. , 2013Paxton et al. , 2015Paxton et al. , 2018Paxton et al. , 2019, we construct a planetary model with a total mass equal to Jupiter's M = M J and a core mass of 5 M ⊕ (though the results should be insensitive to the core as the eigenfunctions of both f-and p-modes are largest near the surface). We then let the model contract until it reaches a desired radius, which we choose to be R = 1.1R J . Irradiation is turned on in this contraction phase with a fixed flux of 6.8×10 6 erg cm −2 s −1 , which corresponds to the average flux the planet receives assuming an orbit with (a orb , e orb ) = (1 AU, 0.98) and a solar-type host star (the equilibrium temperature is 420 K). After the construction of the background model, we find the (adiabatic) eigenmodes (including both the parent f-mode and the leading-order daughter f-and p-modes) using GYRE (Townsend & Teitler 2013;Townsend et al. 2018). The three-mode coupling coefficient is calculated using the expression in Weinberg et al. (2012). We account for the damping of each mode due to turbulent convection using the approach described in Appendix B3 of Burkart et al. (2013).
The key parameters of the Jupiter model are summarized in Table 1. Of particular interest are the values of Ω and Γ. We find that Ω is typically negative and of order |Ω/ω prnt | ∼ 30, where here we will use subscript "a" to stand for a generic mode and "prnt" to indicate the parent mode specifically. The value of Γ will always be positive, as argued in Sec. 3.1.
To illustrate the values of the parameters that enter the calculation of Ω and Γ, we present the eigenfrequency, three-mode-coupling coefficients, and the linear damping rate of each daughter in Figure 1. In addition, we show in Figure 2 each daughter's contribution to Ω (top panel) and Γ (bottom panel). Note that only the negative-frequency, (l a , m a , n a ) = (4, −4, 0) f-mode contributes a positive value to Ω. Although it has the greatest single-mode contribution as it is the most-resonant daughter with respect to the parent's driving, the reso- Here the subscript a in the labels stands for a generic mode, and we specifically label the parent mode's eigenfrequency as ωprnt. For a daughter with |ma| = 4, the coupling we consider is specifically due to parent-parentdaughter, and for an ma = 0 daughter, it is due to parentparent * -daughter. We only show the damping rate for the |ma| = 4 modes as only they contribute to the energy dissipation.
nance is not particularly strong. 11 Instead, upon summing over all couplings, the value of Ω is dominated by the coupling to the many m a = 0 modes ("mode c" in Sec. 3.1) as well as the off-resonant m a = −4 modes ("mode b" with |ω b | > 2|ω a |) and is therefore negative.
By contrast, only the m a = −4 daughters (oscillating at −2ω prnt ) contribute to the damping and Γ is always positive (see both Sec. 3.1). For the model we consider here, Γ has a particularly small value of Γ/ω prnt ∼ 10 −9 . For Jupiters with R > 1.1R J , Γ may be significantly larger. In part, this is because the damping rate due to turbulent convection increases sharply with increasing R (Arras & Socrates 2009). In addition, we find that irradiation causes a thin radiative zone to form near the surface where the daughter p-modes' shears peak. This further reduces the dissipation due to turbulent convection compared to the case without irridiation. . The x-axis is the radial order of the daughter, with na = 0 for f-modes and na > 0 for p-modes. We use different colors to label daughter modes with different (l, m) and solid/dashed lines to label positive/negative-frequency modes. For clarity, we shifted the radial order by +0.1 (−0.1) for modes with positive (negative) frequencies. If a daughter mode contributes a positive value to Ω (Γ), a "+" marker is used; otherwise, we use "-" to show the negative contribution. In the grey, dash-dotted traces, we also show the cumulative values of Ω and Γ obtained by summing over daughter modes with all the possible (la, ma, ωa) and with radial order n a ≥ na. The left-most grey markers thus correspond to the values of Ω and Γ.
From this point onward, we will use values shown in Table 1 as the default parameters for the planetary model. A primary goal of this paper is to develop the theoretical framework for diffusive tidal evolution including nonlinear mode interactions. A more comprehensive survey on how the tidal evolution trajectories depends on different values of (Ω, Γ), and how (Ω, Γ) further vary with respect to (M, R) is deferred to future work.

ITERATIVE MAP INCLUDING NONLINEAR EFFECTS
We now have the ingredients to perform an iterative map similar to the one used by Vick & Lai (2018) but now including nonlinear effects. Since we consider here a single parent mode a, and the daughters' effects are collectively absorbed into Ω and Γ, we will drop the subscript a in mode amplitude and energy from this point onward.
Suppose the parent mode has an amplitude q (1) k−1 in the inertial frame right before the k'th pericenter passage. Its amplitude right after the k'th passage, q (0) k , is given by Eq. (11), just as in the linear case (but see the discussion below of the potential impact of nonlinear effects on the parent's kick). Given q (0) k , the orbital energy and period of the k'th cycle are given by Eqs. (18) and (19), respectively. The next step of the mapping is to relate q (0) k to q (1) k , the amplitude right before the (k + 1)'th passage. While the model we consider in Sec. 3.3 has a particularly weak dissipation, to obtain the accumulated phase over the course of the k'th orbit one needs to account for the gradual decay of the parent's energy due to linear and nonlinear damping. This can be achieved by first obtaining the energyẼ and t = P orb,k . The evolution phase right before the (k + 1)'th passage is then Here δφ nl,k is the excess phase due to the nonlinear frequency shift (the frequency shift due to linear damping is negligible since γ a ω a ), which can be calculated using Eq. (43) withẼ Before we proceed, it is important to point out a few caveats to this approach. First of all, the expressions we derive in this work are only the leading-order nonlinear corrections. They are accurate only when the parent mode's energy satisfiesẼ k 10 −3 . Therefore, in this work our focus will be on the initial triggering of the diffusive tide by the nonlinear mode coupling, particularly the nonlinear phase shift ∝ Ω. The evolution timescale we consider here is thus typically a few hundreds of years or less, when the parent mode is still building up its energy. We defer the examination of tidal evolution over ∼ 10 kyr to follow-up studies in this series, as such a study would require both modifications to our leadingorder expressions, and energy dissipation mechanisms due to both weakly nonlinear damping ∝ Γ and strongly nonlinear wave-breaking as considered in Wu (2018). 13 Secondly, we assumed that the daughters' amplitudes are given by their instantaneous steady-state values, 12 If a system is strictly conservative, then we haveẼ (1) k =Ẽ (0) k , and one can use Eq. (43) with t = P orb,k to obtain the phase. 13 We estimate that one would need Γ/ωa 10 −6 to prevent the parent mode from evolving into the wave-breaking regimeẼ k 0.1 by the weakly-nonlinear damping as we consider here. While this is much greater than the nonlinear damping rate we find for the R = 1.1 R J model, it can nonetheless be achieved if the Jupiter has a greater radius. E.g., we find a Jupiter model R = 2.0 R J can have Γ/ωa 5 × 10 −4 . (30). We show in Appendix A that this may not be strictly true if a daughter mode b (with |m b | = 4) has |2ω a + ω b | Ω peri . Specifically, there should be an additional contribution to the daughter's amplitude that depends on the past history of the mode network. Nonetheless, we drop such corrections for simplicity in the current study. Our numerical experiments suggest this term becomes potentially important only after a few thousand orbital cycles, and therefore should not affect the initial triggering of the diffusive process we consider in this work.

Eqs. (29) and
Lastly, we have assumed that the "kick" at each pericenter passage is always given by the linear calculation. In reality, the kick ∆q 1 depends on the parent's eigenfrequency [through K lm ; Eq. (15)] which changes nonlinearly. Consequently, ∆q 1 should also be modified by the nonlinear frequency shift. However, forẼ k 10 −3 , the fractional decrease of the parent's eigenfrequency ΩEẼ k /ω a is only a few percent. 14 The change in the one-kick amplitude is thus less than 20% according to Eq. (15). Its effect can be more significant as the parent mode's energy further builds up, however; we plan to examine this in follow-up studies.

Relative importance of linear and nonlinear effects for triggering diffusive growth
The main question we want to investigate in this paper is how do nonlinear mode interactions affect the threshold for triggering the diffusive growth of the f-mode? In order to trigger diffusive growth, the phase evolution of the mode must vary randomly from orbit to orbit by an amount (Vick & Lai 2018;Wu 2018) where φ k is given by Eq. (44). In linear theory, this is achieved through the tidal back-reaction on the orbit. At each passage, a random amount of energy ∆Ẽ k is removed from the orbit, which changes the orbital period by ∆P orb,k [Eq. (19)] and consequently the phase by ∆φ br,k = −ω a ∆P orb,k , where the subscript "br" stands for "back-reaction". However, linear theory neglects the fact that the energy ∆Ẽ k gained by the planet's f-mode also changes its eigenfrequency by ∆ω a (∆Ẽ k ) = δω a (Ẽ k ) − δω a (Ẽ k−1 ) Ω∆Ẽ k . This frequency shift induces an additional random phase variation (relative to the previous orbit) ∆φ nl,k − Ω∆Ẽ k P orb,k . The 14 By comparison, the fractional change in Ω peri is only O(10 −5 ) as the mode energy grows from 0 toẼ k 10 −3 ; see Sec. 2.2. nonlinear frequency shift therefore provides another way of triggering the f-mode's diffusive growth.
Quantitatively, Vick & Lai (2018) found that it is sufficient to consider the phase shift after the first pericenter passage in order to determine the boundary for diffusion to happen. Specifically, let ∆Ẽ 1 ≡ |∆q 1 | 2 be the energy gained by the mode after the first pericenter passage (suppose it starts with an amplitude |q 0 | |∆q 1 |). The phase shift caused by the tidal back-reaction after the first pericenter passage is thus 15 where in the second equality we have used the fact that E orb,0 = E orb,0 /E 0 < 0. The threshold for growth is approximately |∆φ br,1 | 1 rad (Vick & Lai 2018), which corresponds to a threshold one-kick energy ∆Ẽ 1 of The nonlinear frequency shift also leads to an excess phase of ∆φ nl,1 − Ω∆Ẽ 1 P orb,0 , where the subscript "nl" stands for "nonlinear" effects, and we have used Eq. (34) for the nonlinear frequency shift. By setting |∆φ nl,1 | = 1 rad we similarly obtain the one-kick energy threshold to trigger diffusive growth soley from nonlinear effects, where we plugged in values representative of the hot Jupiter model described in Section 3.3. Comparing Equations (47) and (49), we see that the nonlinear frequency shift can have a significantly lower one-kick energy threshold than that of tidal back-reaction. It can therefore play a critical role in triggering the diffusive growth of a mode. Furthermore, as we show in Sec. 3.3 that for realistic Jupiter models, Ω < 0 typically and thus ∆φ br,1 ∆φ nl,1 > 0. Intuitively, this can be understood as follows. Suppose a positive amount of energy 15 To be consistent with the indexing convention used in Eq. (19), we should use P orb,1 andẼ orb,1 in Eqs. (46) and (48). Nonetheless, using the quantities evaluated at cycle "0" will only cause a difference of O(∆Ẽ 2 1 ) 1, which can be safely ignored.
is transferred from the orbit to the mode, and as a result the orbital period decreases. Meanwhile, this energy lowers the frequency at which the mode oscillates as Ω < 0 for typical Jupiter models (Table 1). Consequently, both effects make ω a P orb decrease. We thus see the two effects add together to further lower the threshold.
In order to better see how the relative importance of the two effects scales with the various parameters, we take their ratio Note that the ratio is independent of ∆Ẽ 1 and P orb,0 . Instead, it mainly depends on the ratio of orbital energy to binding energy of the planet, ∼ (M * /M )(R/a orb ). Consequently, as a orb decreases during the orbital circularization process, the nonlinear phase shift becomes increasingly dominant over the back reaction shift. This suggests that the nonlinear frequency shift will play a crucial role in maintaining the diffusive energy transfer from the orbit to the planetary mode in the circularization process.

Significance of nonlinear effects in other types of eccentric binaries
We can use Eq. (50) to also estimate the significance of the nonlinear effects in other binary systems with highly eccentric orbits. For a binary of solar-type stars in a highly eccentric orbit with a orb ∼ AU, the nonlinear phase shift is only ∼ 1% of that due to tidal backreaction. 16 A similar ratio of a few percent is also found for a neutron star binary with a orb 1000 km . Indeed, both solar-type stars and neutron stars are more compact (i.e., with smaller M/R) than a typical Jovian planet. Therefore, an eccentric hot Jupiter offers an especially interesting system for studying the impact of nonlinear effects on diffusive growth.

Early orbital evolution following the onset of diffusive growth
16 This is specifically for the nonresonant nonlinear effect of the la = 2 f-mode. A solar-type star also has low-frequency g-modes that have very different values of Ω from that of the f-mode. Those g-modes may allow for a richer family of nonlinear effects, such as the parametric instability.  The top panel shows the mode energy relative to the one kick energy and the bottom panel shows the difference in the mode's excess phase between the k'th and (k − 1)'th orbit. We fix the pericenter distance at Dperi = 3.95 Dt and assume the planet is non-rotating, resulting in a one-kick energy ∆Ẽ1 = 4.5 × 10 −6 . In the linear case (Ω = 0; grey lines) the mode energy and excess phase oscillate periodically and there is no diffusive growth. However, when nonlinear mode interactions are included (Ω/ωa = −29; olive lines) the excess phase is significantly larger and varies randomly, resulting in diffusive growth of the f-mode's energy.
In Fig. 3 we show a representative example of the first 300 orbits of a Jovian planet orbiting a solar-type star at a semi-major axis a orb = 1 AU using the iterative map described in Sec. 4. The planet's parameters are given in Table 1 and we assume here that the planet is not rotating. In the figure, the pericenter distance is set to D peri = 3.95 D t , where D t ≡ R(M * /M ) 1/3 11R J 5.3 × 10 −3 AU is the tidal radius of the planet. The corresponding eccentricity is thus e orb = 1 − D peri /a orb = 0.979. For these parameters, the one-kick energy is ∆Ẽ 1 = |∆q 1 | 2 = 4.5 × 10 −6 [Eq. (16)]. The top panel shows the energy of the parent mode (with l a = m a = 2) and the bottom panel shows the difference of the evolution phase between two adjacent cycles [Eq. (45)]. We see that in the linear case (grey trace), the difference in the mode's propagation phase between adjacent cycles, |∆φ k |, is small (∼ 0.1 rad) and the mode energy just undergoes periodic oscillations (see also the discussions in Vick & Lai 2018). By contrast, when we include nonlinear mode interactions there is an additional contribution to the random phase due to the nonlinear frequency shift [Eq. (48)]. As a result, we see that the f-mode's energy grows diffusively, unlike in the linear case. After 300 cycles, the mode energy grows to about 300∆Ẽ 1 , as one would ex-pect for a diffusive process (i.e., the amplitude grows as the square root of the number of pericenter kicks).
In Fig. 4 we systematically explore some of the conditions necessary to trigger diffusive growth. We show the maximum mode energy achieved after 500 orbital cycles (about 500 years) as a function of the pericenter distance [or equivalently, the orbital eccentricity since D peri = a orb (1 − e orb ) and we set the initial semimajor axis at a orb,0 = 1 AU]. For the given planetary model (Table 1), the one-kick energy ∆Ẽ 1 is shown in the top x-axis for each choice of D peri . In the left panel, we assume the planet is non-rotating, while in the right panel we assume it is pseudo-synchronized with the orbit with Ω s /ω a = 0.19(D peri /0.02 AU) −3/2 and ω a = ω a + 2Ω s . If a mode experiences diffusive growth, then its energy after k pericenter passages is expected to beẼ k ∼ k∆Ẽ 1 on average. Since we set k = 500, we expect max Ẽ k 500∆Ẽ 1 for a mode that grows diffusively (indicated by the black lines), while max Ẽ k ∼ ∆Ẽ 1 for a mode that does not grow (assuming the mode is off resonance with the orbit; in the right panel ω a changes as we vary D peri , allowing it to scan through a series of resonances with different orbital harmonics, thereby causing the excess features to the right of the vertical lines, which we will discuss in Sec. 5.3). We see that each panel in Fig. 4 can be divided up into two regions according to the maximum mode energy achieved. Let us first focus on the left panel (a non-rotating planet with ω a P orb,0 /2π = 2818.73 being a non-integer ). In the linear case (Ω = 0), we find numerically that the boundary where diffusive growth is first triggered is at D peri 3.75D t , corresponding to a one-kick energy of ∆Ẽ 1 1.4 × 10 −5 . The analytical estimate [Eq. (47); see vertical grey line], agrees well with the numerical results but slightly overestimates the threshold value of D peri (and underestimates ∆Ẽ 1 ) because there we simply assumed the threshold phase shift to be 1 rad; in reality a slightly greater phase shift is required. This can also be seen from the bottom panel where we show the fraction of systems undergoing diffusive growth (i.e., the fraction of points around the black lines; the estimate is preformed over a full bin width of 0.1D t ). The error bars around the vertical lines are obtained by setting ∆φ 1 = 0.5 rad and 1.5 rad and then re-evaluating ∆Ẽ 1 using Eqs. (46) and (48).
When we account for the nonlinear frequency shift, we find that the boundary moves to larger D peri (smaller one-kick energies). For the representative value of Ω/ω a −29 (see Table 1), we find that the threshold one-kick energy is lowered to ∆Ẽ 1 2.4×10 −6 , which is a factor of about 6 smaller than the linear case [see the vertical olive line and Eq. (49); note that the threshold is in fact determined by the sum ∆φ nl,1 + ∆φ br,1 , with the latter being 20% of the former for the parameters in Fig. 4]. Because the one-kick energy depends sensitively on the pericenter distance, a factor of six change in ∆Ẽ 1 corresponds to a 10% increase in D peri .

Including spin effects
We consider the effects of planet spin in the right panel of Fig. 4. We assume the planet is rotating at a rate determined by the pseudo-synchronization condition [Eq. (23)]. 17 As we show in Sec. 2.1, the mapping equations including spin are formally the same as the non-spinning equations except that the mode frequency is replaced by the inertial frame value ω a = ω a + m a Ω s (this frequency then enters the calculations of the onekick amplitude and the linear propagation phase). Since we focus on a prograde mode with m a = 2, |∆q 1 | ∝ K 22 decreases sharply as ω a increases [Eq. (15)]. As a result, the D peri boundary where diffusive growth is first triggered is shifted to smaller values.
At the same time, in the right panel the mode can occasionally become resonant with the orbit when ω a P orb /2π = integer, as ω a now varies with D peri due to pseudo-synchronization condition (this is in contrast to the left panel where ω a P orb /2π is fixed at a noninteger value when we vary D peri ). For D peri /D t 3.8, such resonances can bring the mode energy up toẼ (1) k a few × 10 −5 ∆Ẽ 1 . However, as the mode acquires energy from the orbit, the orbital period starts to change (though not by a significant enough amount to trigger diffusion). It thus destroys the resonance and prevents the mode energy from increasing further (see also Vick & Lai 2018). Similarly, the nonlinear frequency shift also destroys the resonance between ω a and P orb and this is why the upper envelope of the olive dots is at a lower value than that of the grey ones.
At 3.2 D peri /D t 3.8, we see that the chance resonance with the orbit may occasionally help a slightly sub-threshold mode to also evolve into the diffusive regime. Suppose the chance resonance helps the mode to initially build up an energyẼ k,0 10 −5 (10 −4 ) with (without) the nonlinear effect (corresponding approximately to the upper envelopes shown in the right panel). The typical energy exchange between a mode and the orbit is then given by Ẽ k,0 ∆Ẽ 1 ∆Ẽ 1 (Wu 2018).  Top panels: maximum f-mode energy achieved after 500 pericenter passages. The grey circles are calculated assuming linear theory while the olive circles also include nonlinear mode interactions (with Ω/ωa = −29). In both plots we set a orb,0 = 1 AU. In the left panel we assume that the planet is not spinning while in the right panel we assume that it spins at a constant rate given by the pseudo-synchronization condition [Eq. (23)]. Bottom panels: the fraction of systems that undergoes the diffusive growth (i.e., points around the black lines) as a function of Dperi, estimated over a full bin width of 0.1 Dt. The vertical lines are the analytic estimates for the diffusive growth threshold based on Eqs. (47) and (49) by setting ∆φ1 = 1 rad. The error bars are obtained if we instead use ∆φ1 = 0.5 rad or 1.5 rad. The region where we expect diffusive growth to occur are also indicated by arrows in the top panel. In the right panel, a mode to the right of the vertical lines can occasionally grow to an intermediate amplitude ofẼ (1) k a few × 10 −5 if its frequency ω a comes to close resonance with one of the orbital harmonics.
If we replace ∆Ẽ 1 by Ẽ k,0 ∆Ẽ 1 in Eqs. (46) and (48), we see the new threshold one-kick energy becomes ∆Ẽ br,1 10 −6 and ∆Ẽ nl,1 3 × 10 −7 for modes that are initially in close resonance with the orbit. As D peri decreases and ∆Ẽ 1 increases, even systems that are not at the upper envelope may enter the chaotic regime and the fraction of diffusive systems increases as indicated by the lower panel. Eventually, when ∆Ẽ 1 reaches the value derived in Eqs. (47) and (49), almost all of the systems will grow diffusively.

Threshold expressed in terms of semi-major axis rather than one-kick energy
An alternative way to consider the problem is to hold D peri and thus ∆Ẽ 1 fixed and instead vary the initial semi-major axis a orb,0 = D peri /(1 − e orb,0 ). By setting |∆φ br,1 | = 1 rad in Eq. (46) as before (Section 5.1) but now solving for a (br) orb,0 , we find that the threshold to trig-ger diffusive growth due to only tidal backreaction is a (br) orb,0 2.5 AU ∆Ẽ 1 10 −6 Similarly, we can use Eq. (48) to find the threshold due to nonlinear mode interactions a (nl) orb,0 In both cases, the threshold a orb,0 increases with decreasing one-kick energy ∆Ẽ 1 (i.e., increasing D peri ). This is because both ∆φ br and ∆φ nl ∝ P orb ∆Ẽ 1 [Eq. (46) and (48)], and a longer P orb ∝ a 3/2 orb is thus required in order for the f-mode to accumulate an excess phase |∆φ k | to O(1) rad. Additionally, Eqs. (51) and (52) (51) and (52), respectively, assuming ∆φ1 = 1 rad is needed to trigger diffusive growth; the error bars show the threshold if instead ∆φ1 = 0.5 rad or 1.5 rad are needed. The fraction of systems undergoing diffusive growth shown in the bottom panel is estimated over logarithmic bins with full width of log 10 (a orb,0 /AU) = 0.4. A system with an a orb,0 that is slightly below the threshold may still trigger diffusive growth if the mode is close to resonance with the orbit.
|E orb | −1 , which reflects the fact that an orbit with greater a orb is less bound and thus sees a greater change in the fractional orbital period. Also note that Eq. (52) overestimates the minimum a orb,0 required to trigger the diffusion because it assumes only the nonlinear contribution to ∆φ k . In reality, the nonlinear frequency shift and the tidal back-reaction both contribute to ∆φ k and they have the same sign (Section 5.1). Similar to the case shown in the right panel of Fig. 4, Eqs. (51) and (52) should be treated as the upper end of the thresholds; almost all of the systems with a orb,0 greater than the values estimated in Eqs. (51) and (52) will trigger the diffusive evolution. On the other hand, if a system is initially close to being resonant with the orbit, then at smaller a orb,0 it may still enter the diffusive regime. We evaluate the boundary in a orb,0 numerically in Fig. 5. Here we fix the pericenter distance to be D peri = 0.02 AU, corresponding to a one-kick energy of ∆Ẽ 1 1.4×10 −5 (∆Ẽ 1 6.7×10 −7 ) for a non-spinning (pseudo-synchronized) planet. The situation is particularly interesting astrophysically in the case where the planet's spin is pseudo-synchronized. For a relatively weak one-kick energy of ∆Ẽ 1 10 −6 , only planets born 2 AU away from the host star can trigger diffusive tidal evolution and form hot Jupiters if only the lin-ear theory is used. On the other hand, when we include nonlinear mode interactions, it can be triggered for planets born with a orb,0 0.7 − 2 AU. Thus, nonlinear mode interactions significantly expand the parameter space allowed for diffusive growth to happen, which not only allows more potential progenitors to form hot Jupiters within the age of the Universe, but also saves more planets from disruption by the host star during the Kozai cycles (see, e.g., Vick et al. 2019). It could thus help alleviate the discrepancy between the predicted and observed hot Jupiter to regular Jupiter occurrence rate (see Dawson & Johnson 2018).

CONCLUSION AND DISCUSSION
We studied the nonlinear interaction between a selfcoupled parent f-mode and daughter f-and p-modes in a Jovian planet. For a parent mode with azimuthal quantum number m a = 2 and frequency ω a , it drives both m b = −4 daughters that correspond to waves oscillating at 2ω a , and non-oscillatory m c = 0 daughters that correspond to a modification of the planet's structure (Sec. 3). We found that at leading order, the interaction leads to a nonlinear shift in the parent mode's eigenfrequency, δω a , as well as a nonlinear increase in the parent mode's damping rate (imaginary part of the frequency), δγ a . Both the nonlinear frequency shift and damping rate follow the scaling δω a , δγ a ∝ ω a E a at leading order [Eqs. (34) and (35); see also Kumar et al. (1994); Kumar & Goodman (1996)]. The modifications are timedependent because we consider planets on highly eccentric orbits with parent mode energies E a that vary at each pericenter passage. Furthermore, we showed that although the frequency shift can, in principle, be either positive or negative, for typical Jupiter models a negative shift is more likely (that is, the parent mode's eigenfrequency decreases with increasing mode energy). The nonlinear damping, on the other hand, strictly increases as the mode energy increases (see Table 1 and Fig. 2).
We then developed the formalism to construct iterative maps including nonlinear effects and applied them to study how nonlinear interactions affect the higheccentricity migration of proto-hot Jupiters. We found that the energy-dependent nonlinear frequency shift leads to an excess phase of the parent mode [Eq. (48)] which is stochastic from orbit to orbit. It thus provides another channel for triggering the diffusive growth of the parent mode in addition to the tidal back-reaction considered in previous studies.
In fact, we found that for typical Jupiter models, the nonlinear phase shift is ≈ 5 times larger than the phase shift due to back-reaction [Eq. (50)]. The two effects add together and lower the threshold one-kick energy required in order to trigger the growth by about a factor of ≈ 6 compared to the case without nonlinear interactions (Fig. 4). Alternatively, if one fixes the one-kick energy, the threshold on the minimum initial orbital semi-major axis can be lowered by a factor of ≈ 2 (Fig. 5). If the one-kick energy is small (due to either a small eccentricity and hence large pericenter distance, or a high spin rate of the planet), then in the linear case only planets born at a orb,0 2 AU can undergo diffusive tidal evolution and form hot Jupiters; however, when nonlinear interactions are accounted for, it is lowered to the interesting range of a orb,0 = 0.7 − 2 AU.
In this paper, we focused on developing the theoretical framework and considered only the evolution over the first O(100) yr. There are several aspects of the problem we think would be interesting to address in future studies.
First, what is the long-term evolution of the system over ∼ 10 kyr? For the Jupiter model we considered in this work (Table 1), the weakly nonlinear damping ∝ Γ is weak and thus the parent mode energy will grow so large that it likely becomes strongly nonlinear, as assumed by Wu (but see discussion below). Wu found that the diffusive process, and hence the orbital evolution, typically stalls when the semi-major axis decays to a orb 0.2 AU while the eccentricity is still high (e orb 0.9), and it was unclear what drives the subsequent orbital circularization. However, nonlinear mode interactions might prevent the circularization from stalling at high e orb because the magnitude of the random phase it induces decreases slower than that due to the tidal back-reaction as the orbit decays; see Eq. (50). This is because the orbit "hardens" (|E orb | increases) as a orb shrinks, which makes it increasingly hard to be perturbed by the tidal back-reaction [Eqs. (19) and (46)]. By contrast, the natural energy scale that enters the nonlinear phase shift is the a orb -independent binding energy of the planet (ignoring the evolution of the planet). An efficient circularization could help explain both the paucity of super-eccentric Jupiters (Socrates et al. 2012) and the relatively young age of hot-Jupiter host stars (Hamer & Schlaufman 2019).
It would also be interesting to investigate how the values of (Ω, Γ) vary for different Jupiter models and how the tidal evolution trajectories depend on (Ω, Γ). We estimate that if Γ 10 −6 ω a , weakly nonlinear damping could be sufficient to prevent the parent mode from evolving into the strongly nonlinear regime. This would lead to another qualitative difference from the trajectories found in Wu (2018), in addition to the excess nonlinear phase shift discussed above. Such large dissipation rates can be achieved by Jupiter models with greater radii and it could have potentially important observational consequences.
A calculation that combines diffusive tidal evolution with the mechanism that drives the eccentricity to large values in the first place (e.g., Lidov-Kozai cycles) would be valuable and help test these ideas further. By including nonlinear effects, it would extend the work of Vick et al. (2019) and thereby provide a more robust estimate of the formation rate of hot Jupiters due to diffusive tidal evolution. Current theories produce too few hot Jupiters relative to regular Jupiters, and it would be interesting to know whether nonlinear effects could help mitigate the tension.
To carry out the studies described above, a few modifications to the current framework would be needed. For instance, as the parent mode's energy builds up and its eigenfrequency decreases, the orbital integral K lm should be modified accordingly. Since the parent's frequency is typically shifted to a lower value (Ω < 0), we would expect the one-kick amplitude ∆q 1 ∝ K lm to increase in magnitude as the parent's energy increases [Eq. (15)]. This would further enhance the significance of the nonlinear effects. On the other hand, we do not expect a linear-in-energy frequency shift [Eq. (34)] to be accurate when |Ω|Ẽ a ω a . Note that this condition can happen at a smaller energy than the wave-breaking energy, and therefore further corrections would be needed. ACKNOWLEDGMENTS  (30) we assumed the daughter modes' amplitudes are given by their steady-state values. While this is a good approximation for 'mode c' (daughters with m c = 0), as we explain here the problem may be more involved for 'mode b' (daughters with l b = −m b = 4).
The equation governing such a mode b is given by [see Eq. (25)] where, as explained in Section 3.1 (also see discussion below), we can ignore the linear tidal forcing on mode b, i.e., the U b term. We can decompose the parent mode (mode a) as q a = q a,dyn + q a,eq = q a,dyn + U a , where q a,eq ≡ U a is the equilibrium tide solution of mode a [which can be obtained from Eq. (24) when we ignore the nonlinear couplings and treat |q a |, |γ a q a | |ω a q a |]. We thus havė In the main text, we focused on the steady-state solution of q b driven by a free-oscillating parent. That is, we assumed q a contains only the dynamical component q a,dyn which oscillates at a single frequency ω a , and found [see Eq. (29)] where ∆ b = (ω b + 2ω a ). Note, however, that the steady state solution q b,ss neglects the U * a U * a term in Eq. A3 and it neglects the 'transient' part of the solution for q b . We will refer to the latter as the history term since it depends on the past history of q b (t) from previous pericenter passages. Here we will show that the U * a U * a term should always be insignificant but not necessarily the history term.
We will make two simplifications in our analysis. First, we do not explicitly solve for the instantaneous value of q b in the vicinity of a pericenter passage for simplicity. As we will see, this does not preclude us from obtaining a qualitative estimate of the history term due to previous pericenter passages. Second, we treat the parent mode as if it is unperturbed by nonlinear interactions. As a result, we assume that the dynamical component of the parent mode, q a,dyn , when far away from the pericenter, oscillates at ω a and not at the nonlinearly shifted frequency (cf. Sec. 3.1).
A.1. The U * a U * a term First we consider the drive due to the U * a U * a term. Similar to the one-kick amplitude of the parent ∆q a,1 , we can define a one-kick amplitude of mode b at each pericenter passage due to the U * a U * a term as where we used the fact that l b = 2l a and m b = −2m a . If we define as a modified temporal overlap, then we can write the one-kick amplitude of the daughter mode as where we have plugged in l a = m a = 2 for the parent and l b = −m b = 4 for the daughter. For typically values (a orb =1 AU, e orb =0.98, ω a =5.6 × 10 −4 rad/s, and ω b = − 1.9ω a ), we find (∆q b ) U 2 a 0.4|∆q 2 a,1 |. Although intially (∆q b,1 ) U 2 a is comparable to the steady-state solution q b,ss [Eq. (A4)], as the system starts to grow diffusively, its effect will soon become subdominant. This can be seen by noticing that even if (∆q b,1 ) U 2 a can grow diffusively itself (e.g., due to a random P orb ), after k pericenter passages, it only increase the amplitude of q b by √ k| (∆q b,1 ) U 2 a | ∝ k 1/2 on average. On the other hand, q b,ss ∼ k|q a,dyn | 2 ∝ k because each |q a,dyn | grows as √ k.
Consequently, the significance of the q b,ss term increases as √ k. In fact, the dominance of q b,ss is further enhanced by the ω b /∆ b factor, especially for the most resonant daughter mode with the smallest |2ω a + ω b |.
We therefore conclude that the modification to q b due to the (U * a ) 2 term can be ignored. It is also worth noting that the drive from (U * a ) 2 is stronger than the direct tidal force on mode b, U b ∝ (M * /M )(R/D peri ) l b +1 , by a factor of (M * /M )(R/D peri ) κ b W 2 22 Q 2 a /W 44 Q b 6 κ b W 2 22 Q 2 a /W 44 Q b . It thus justifies why we can also ignore the daughter modes' linear coupling to the tide.
A.2. The q * a,dyn q * a,dyn term We now consider the effect of the history term on the daughter. If we define c b = q b exp(−2iω a t) then by Eq. (A3), where V b ≡ κ b q * a q * a exp(−2iω a t). Our definition of V b does not include the equilibrium tide contribution U * a U * a and U b since we showed above that they are insignificant. For the same reason, here and below we drop the "dyn" subscript on the parent. Note that if we ignore the parent's nonlinear frequency corrections, then away from pericenter q a ∼ e −iωat and thus V b is a constant. Near pericenter, however, q a has an additional time-dependence due to the kick ∆q a,1 the parent receives over a timescale 1/Ω peri . As we will see, it is this effect that constitutes the history term we are interested in.
The general solution for c b is given by where the initial time is t 0 and we performed integration by parts to get the second line. For future convenience we set V b (t 0 ) = 0 and thus drop the initial condition. Note that the first term in Eq. (A9) recovers the steady-state solution, Eq. (29), and it depends only on the instantaneous value of q a . The second term, on the other hand, captures the past history. For a free oscillator, q a ∼ e −iωat andV b = 0, and thus the value of q b is independent of the past history.
When the system is coupled to the tide, however, we haveV b ∼ Ω peri V b in the vicinity of the pericenter. First consider a mode b for which |∆ b | Ω peri γ b . These inequalities hold for all the daughters in our mode networks with the exception of the l b = −m b = 4, ω b < 0 f-mode, which we consider separately below. For large detuning, if we keep performing integration by parts, we get where we dropped γ b to simplifiy the notation. Note that after the n'th iteration of integration by parts, we have a correction that depends only on the instantaneous value ∝ V , and an integrand ∝ V is the n'th time derivative of V b . Since |∆ b | Ω peri , the history-dependent term gets progressively smaller with increasing n, and the instantaneous corrections to Eq. (29) form a converging series (in fact, the corrections are