Tides in the High-eccentricity Migration of Hot Jupiters: Triggering Diffusive Growth by Nonlinear Mode Interactions

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 Equation (3) from the corotating frame to the inertial frame

$$\begin{eqnarray}&&{\dot{q}}_{a}^{{\prime} }+\left(i{\omega }_{a}^{{\prime} }+{\gamma }_{a}\right){q}_{a}^{{\prime} }=i{\omega }_{a}{U}_{a}^{{\prime} }(t),\end{eqnarray} \tag{ 9 }$$

where ${q}_{a}^{{\prime} }={q}_{a}\exp \left(-{{im}}_{a}{{\rm{\Omega }}}_{{\rm{s}}}t\right)$, ${\omega }_{a}^{{\prime} }={\omega }_{a}+{m}_{a}{{\rm{\Omega }}}_{{\rm{s}}}t$, and ${U}_{a}^{{\prime} }\,={U}_{a}\exp \left(-{{im}}_{a}{{\rm{\Omega }}}_{{\rm{s}}}t\right)$ are, respectively, the mode amplitude, mode frequency, and tidal driving in the inertial frame. The general solution of ${q}_{a}^{{\prime} }$ is

$$\begin{eqnarray}&&{q}_{a}^{{\prime} }(t)={e}^{-\left(i{\omega }_{a}^{{\prime} }+{\gamma }_{a}\right)t}{\int }^{t}i{\omega }_{a}{U}_{a}^{{\prime} }(\tau ){e}^{\left(i{\omega }_{a}^{{\prime} }+{\gamma }_{a}\right)\tau }d\tau .\end{eqnarray} \tag{ 10 }$$

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)

$$\begin{eqnarray}&&{q}_{a,k}^{{\prime} (0)}={q}_{a,k-1}^{{\prime} (1)}+{\rm{\Delta }}{q}_{a,1}\end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}&&{q}_{a,k}^{{\prime} (1)}={q}_{a,k}^{{\prime} (0)}{e}^{-\left(i{\omega }_{a}^{{\prime} }+{\gamma }_{a}\right){P}_{\mathrm{orb},k}},\end{eqnarray} \tag{ 12 }$$

where the superscripts (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 ⁴

$$\begin{eqnarray}&&{\rm{\Delta }}{q}_{a,1}=\int i{\omega }_{a}{U}_{a}^{{\prime} }(\tau ){e}^{(i{\omega }_{a}^{{\prime} }+{\gamma }_{a})\tau }d\tau .\end{eqnarray} \tag{ 13 }$$

In the equation above, the integration is performed over one orbital period. For the rest of the paper, we will only consider highly eccentric orbits with $\left(1-{e}_{\mathrm{orb}}\right)\ll 1$, where significant tidal interaction happens only near the pericenter. Therefore, the limits of integration in Equation (13) can be dropped as long as they bracket the pericenter passage.

It is convenient to define an orbital integral K_lm as ⁵ (Press & Teukolsky 1977)

$$\begin{eqnarray}&&{K}_{{lm}}(\omega )=\displaystyle \frac{{\omega }_{0}{W}_{{lm}}}{2\pi }\int {\left[\displaystyle \frac{{D}_{\mathrm{peri}}}{D(\tau )}\right]}^{l+1}{e}^{i\left[\omega \tau -m{\rm{\Phi }}(\tau )\right]}d\tau ,\end{eqnarray} \tag{ 14 }$$

where D_peri ≡ a_orb(1 − e_orb) is the pericenter distance, e_orb the eccentricity, and ${\omega }_{0}\equiv \sqrt{{GM}/{R}^{3}}$. If we ignore the perturbations on D and Φ, Lai (1997) provide an analytical expression for K₂₂ (i.e., l = m = 2) assuming that the orbit is parabolic,

$$\begin{eqnarray}\begin{array}{rcl}{K}_{22}(\omega ,{{\rm{\Omega }}}_{\mathrm{peri}}) & \simeq & \displaystyle \frac{2{z}^{3/2}{e}^{-2z/3}}{\sqrt{15}}\left(\displaystyle \frac{{\omega }_{0}}{{{\rm{\Omega }}}_{\mathrm{peri}}}\right)\left(1-\displaystyle \frac{\sqrt{\pi }}{4\sqrt{z}}\right),\\ & \simeq & 1.1\times {10}^{-2}{\left(\displaystyle \frac{z}{11}\right)}^{-5.7}\left(\displaystyle \frac{{\omega }_{0}}{{{\rm{\Omega }}}_{\mathrm{peri}}}\right),\end{array}\end{eqnarray} \tag{ 15 }$$

where $z\equiv \sqrt{2}\omega /{{\rm{\Omega }}}_{\mathrm{peri}}$ and ${{\rm{\Omega }}}_{\mathrm{peri}}^{2}\equiv G(M+{M}_{* })/{D}_{\mathrm{peri}}^{3}$, and in the second line we have expanded the expression around z = 11 to emphasize the steep decline.

In terms of K_lm , we can write the one-kick amplitude as

$$\begin{eqnarray}\begin{array}{rcl}{\rm{\Delta }}{q}_{a,1} & = & i2\pi {Q}_{a}{K}_{{lm}}({\omega }_{a}^{{\prime} },{{\rm{\Omega }}}_{\mathrm{peri}})\\ & & \times \left(\displaystyle \frac{{\omega }_{a}{{\rm{\Omega }}}_{\mathrm{peri}}^{2}}{{\omega }_{0}^{3}}\right)\left(\displaystyle \frac{{M}_{* }}{M+{M}_{* }}\right){\left(\displaystyle \frac{R}{{D}_{\mathrm{peri}}}\right)}^{l-2}.\end{array}\end{eqnarray} \tag{ 16 }$$

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 ${\omega }_{a}^{{\prime} }=({\omega }_{a}+m{{\rm{\Omega }}}_{{\rm{s}}})$. Combining with the expansion given in Equation (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)

$$\begin{eqnarray}&&{\rm{\Delta }}{E}_{a,k}=\left[{\left|{q}_{a,k}^{{\prime} (0)}\right|}^{2}-{\left|{q}_{a,k-1}^{{\prime} (1)}\right|}^{2}\right]{E}_{0}.\end{eqnarray} \tag{ 17 }$$

Since the energy stored in the tidal coupling (the term $\propto \mathrm{Re}\left[{q}_{a}^{* }{U}_{a}\right]$) 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,

$$\begin{eqnarray}&&{\rm{\Delta }}{E}_{\mathrm{orb},k}={E}_{\mathrm{orb},k}-{E}_{\mathrm{orb},k-1}=-{\rm{\Delta }}{E}_{a,k},\end{eqnarray} \tag{ 18 }$$

where E_orb,k = −GMM_* /2a_orb,k is the orbital energy at the kth orbit.

A direct consequence is that the change in the orbital energy also alters the orbital period ${P}_{\mathrm{orb}}=2\pi /{{\rm{\Omega }}}_{\mathrm{orb}}\,=2\pi \sqrt{{a}_{\mathrm{orb}}^{3}/G(M+{M}_{* })}\propto {\left(-{E}_{\mathrm{orb}}\right)}^{-3/2}$, as

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{\Delta }}{P}_{\mathrm{orb},k}}{{P}_{\mathrm{orb},k}}=\displaystyle \frac{3}{2}\displaystyle \frac{{\rm{\Delta }}{E}_{a,k}}{{E}_{\mathrm{orb},k}}.\end{eqnarray} \tag{ 19 }$$

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

$$\begin{eqnarray}&&{\rm{\Delta }}{\phi }_{\mathrm{br},k}=-{\omega }_{a}{\rm{\Delta }}{P}_{\mathrm{orb},k}=-\displaystyle \frac{3}{2}{\omega }_{a}{P}_{\mathrm{orb},k}\displaystyle \frac{{\rm{\Delta }}{E}_{a,k}}{{E}_{\mathrm{orb},k}},\end{eqnarray} \tag{ 20 }$$

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.

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)

$$\begin{eqnarray}&&\left|\,\displaystyle \frac{{\rm{\Delta }}{L}_{\mathrm{orb}}}{{L}_{\mathrm{orb}}}\,\right|\simeq \left|\,\displaystyle \frac{{\rm{\Delta }}{E}_{a}(1-{e}_{\mathrm{orb}})}{2\sqrt{2}{E}_{\mathrm{orb}}}\,\right|\ll \left|\,\displaystyle \frac{{\rm{\Delta }}{E}_{a}}{{E}_{\mathrm{orb}}}\,\right|.\end{eqnarray} \tag{ 21 }$$

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) owing 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

$$\begin{eqnarray}\begin{array}{rcl}{L}_{\mathrm{orb}} & = & \mu \sqrt{G(M+{M}_{* }){a}_{\mathrm{orb}}(1-{e}_{\mathrm{orb}}^{2})}\\ & \simeq & \mu \sqrt{2G(M+{M}_{* }){D}_{\mathrm{peri}}}\end{array}\end{eqnarray} \tag{ 22 }$$

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)

$$\begin{eqnarray}\begin{array}{rcl}{{\rm{\Omega }}}_{{\rm{s}}} & \simeq & \displaystyle \frac{1+\tfrac{15}{2}{e}_{\mathrm{orb}}^{2}+\tfrac{45}{8}{e}_{\mathrm{orb}}^{4}+\tfrac{5}{16}{e}_{\mathrm{orb}}^{6}}{\left(1+3{e}_{\mathrm{orb}}^{2}+\tfrac{3}{8}{e}_{\mathrm{orb}}^{4}\right){\left(1-{e}_{\mathrm{orb}}^{2}\right)}^{3/2}}{{\rm{\Omega }}}_{\mathrm{orb}},\\ & \simeq & 1.17{{\rm{\Omega }}}_{\mathrm{peri}}\quad \left({e}_{\mathrm{orb}}\to 1\right),\end{array}\end{eqnarray} \tag{ 23 }$$

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 pseudo-synchronous 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 nonrotating planet case. Nevertheless, the one-kick amplitude of the prograde mode is still two orders of magnitude greater than 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 procedure, 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 ${\omega }_{a}^{{\prime} }={\omega }_{a}+{m}_{a}{{\rm{\Omega }}}_{{\rm{s}}}$ and neglects the nonlinear effects described in the next section.

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 f-mode or a p-mode), which results in the inhomogeneous driving of the daughter. As we consider fully convective Jupiter models, there are no low-frequency g-modes that can parametrically couple to the parent (Weinberg et al. 2012). Therefore, the nonresonant 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. ⁶

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. ⁷

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. We then have

$$\begin{eqnarray}&&{\dot{q}}_{a}+(i{\omega }_{a}+{\gamma }_{a}){q}_{a}=2i{\omega }_{a}{\kappa }_{b}{q}_{b}^{* }{q}_{a}^{* }+2i{\omega }_{a}{\kappa }_{c}{q}_{c}^{* }{q}_{a},\end{eqnarray} \tag{ 26 }$$

$$\begin{eqnarray}&&{\dot{q}}_{b}+(i{\omega }_{b}+{\gamma }_{b}){q}_{b}=i{\omega }_{b}{\kappa }_{b}{q}_{a}^{* }{q}_{a}^{* },\end{eqnarray} \tag{ 27 }$$

$$\begin{eqnarray}&&{\dot{q}}_{c}+(i{\omega }_{c}+{\gamma }_{c}){q}_{c}=2i{\omega }_{c}{\kappa }_{c}{q}_{a}^{* }{q}_{a}.\end{eqnarray} \tag{ 28 }$$

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 Equations (26)–(28) in a perturbative manner. Away from pericenter and without nonlinear couplings, we have ${q}_{a}\sim \exp (-i{\omega }_{a}t)$. Using this as the driving term, the steady-state solutions of the daughters are ⁸

$$\begin{eqnarray}&&{q}_{b}=\displaystyle \frac{{\omega }_{b}\left(2{\omega }_{a}+{\omega }_{b}\right)+i{\omega }_{b}{\gamma }_{b}}{{\left(2{\omega }_{a}+{\omega }_{b}\right)}^{2}+{\gamma }_{b}^{2}}{\kappa }_{b}{q}_{a}^{* }{q}_{a}^{* },\end{eqnarray} \tag{ 29 }$$

$$\begin{eqnarray}&&{q}_{c}=\displaystyle \frac{2{\omega }_{c}^{2}+2i{\omega }_{c}{\gamma }_{c}}{{\omega }_{c}^{2}+{\gamma }_{c}^{2}}2{\kappa }_{c}{q}_{a}^{* }{q}_{a}.\end{eqnarray} \tag{ 30 }$$

Plugging the daughter modes above back into Equation (26), we obtain

$$\begin{eqnarray}&&{\dot{q}}_{a}+\left[i\left({\omega }_{a}+\delta {\omega }_{a}\right)+\left({\gamma }_{a}+\delta {\gamma }_{a}\right)\right]{q}_{a}=0,\end{eqnarray} \tag{ 31 }$$

where ⁹

$$\begin{eqnarray}&&\delta {\omega }_{a}=-{\omega }_{a}\left[\displaystyle \frac{2{\omega }_{b}\left(2{\omega }_{a}+{\omega }_{b}\right)}{{\left(2{\omega }_{a}+{\omega }_{b}\right)}^{2}+{\gamma }_{b}^{2}}{\kappa }_{b}^{2}+\displaystyle \frac{4{\omega }_{c}^{2}}{{\omega }_{c}^{2}+{\gamma }_{c}^{2}}{\kappa }_{c}^{2}\right]{\tilde{E}}_{a}\end{eqnarray} \tag{ 32 }$$

$$\begin{eqnarray}&&\delta {\gamma }_{a}=-{\omega }_{a}\left[\displaystyle \frac{2{\omega }_{b}{\gamma }_{b}}{{\left(2{\omega }_{a}+{\omega }_{b}\right)}^{2}+{\gamma }_{b}^{2}}{\kappa }_{b}^{2}+\displaystyle \frac{4{\omega }_{c}{\gamma }_{c}}{{\omega }_{c}^{2}+{\gamma }_{c}^{2}}{\kappa }_{c}^{2}\right]{\tilde{E}}_{a},\end{eqnarray} \tag{ 33 }$$

and ${\tilde{E}}_{a}\equiv {q}_{a}^{* }{q}_{a}$ is the dimensionless energy of mode a. ¹⁰ The physical mode energy including the contribution from both a and its complex conjugate a* is ${E}_{a}={\tilde{E}}_{a}{E}_{0}$ in our normalization, where E₀ = GM²/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

$$\begin{eqnarray}&&\delta {\omega }_{a}({\tilde{E}}_{a})={\rm{\Omega }}{\tilde{E}}_{a},\end{eqnarray} \tag{ 34 }$$

$$\begin{eqnarray}&&\delta {\gamma }_{a}({\tilde{E}}_{a})={\rm{\Gamma }}{\tilde{E}}_{a},\end{eqnarray} \tag{ 35 }$$

where, after putting back the summation over all the daughters that couple to mode a, we have

$$\begin{eqnarray}&&{\rm{\Omega }}=-{\omega }_{a}\left[\sum _{b}\displaystyle \frac{2{\omega }_{b}\left(2{\omega }_{a}+{\omega }_{b}\right)}{{\left(2{\omega }_{a}+{\omega }_{b}\right)}^{2}+{\gamma }_{b}^{2}}{\kappa }_{b}^{2}+\sum _{c}\displaystyle \frac{4{\omega }_{c}^{2}}{{\omega }_{c}^{2}+{\gamma }_{c}^{2}}{\kappa }_{c}^{2}\right],\end{eqnarray} \tag{ 36 }$$

$$\begin{eqnarray}&&{\rm{\Gamma }}=\sum _{b}\displaystyle \frac{-2{\omega }_{a}{\omega }_{b}}{{\left(2{\omega }_{a}+{\omega }_{b}\right)}^{2}+{\gamma }_{b}^{2}}{\gamma }_{b}{\kappa }_{b}^{2}.\end{eqnarray} \tag{ 37 }$$

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 7). Nonetheless, mode c is not a wave itself (as it is nonoscillatory), 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 ¹¹ act to reduce the parent mode's frequency, a mode b with ω_b < 0 and $\left(2{\omega }_{a}+{\omega }_{b}\right)\gt 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 ${\omega }_{b}\left(2{\omega }_{a}+{\omega }_{b}\right)\lt 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 $\left(2{\omega }_{a}+{\omega }_{b}\right)\simeq {0}^{+}$. ¹²

In the previous section, we showed that nonlinear mode interactions perturb the frequency of the f-mode by $\delta {\omega }_{a}({\tilde{E}}_{a})$ and its damping rate by $\delta {\gamma }_{a}({\tilde{E}}_{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 owing 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 ${\tilde{E}}_{a}$, in order to determine δϕ_nl as a function of time, we need to determine how ${\tilde{E}}_{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 ${\tilde{E}}_{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 ¹³

$$\begin{eqnarray}&&{\dot{\tilde{E}}}_{a}+2\left[{\gamma }_{a}+\delta {\gamma }_{a}({\tilde{E}}_{a})\right]{\tilde{E}}_{a}=0.\end{eqnarray} \tag{ 38 }$$

If we substitute in Equation (35) for $\delta {\gamma }_{a}({\tilde{E}}_{a})$, then the above equation can be solved easily as

$$\begin{eqnarray}\begin{array}{rcl} & {\tilde{E}}_{a}(t) & =\,\frac{{\gamma }_{a}{\tilde{E}}_{a}^{(0)}}{-{\rm{\Gamma }}{\tilde{E}}_{a}^{(0)}+\left[{\gamma }_{a}+{\rm{\Gamma }}{\tilde{E}}_{a}^{(0)}\right]{e}^{2{\gamma }_{a}t}}\\ & & \simeq \frac{{\tilde{E}}_{a}^{(0)}}{1+2\left[{\gamma }_{a}+{\rm{\Gamma }}{\tilde{E}}_{a}^{(0)}\right]t},\end{array}\end{eqnarray} \tag{ 39 }$$

where ${\tilde{E}}_{a}^{(0)}$ is the initial mode energy and in the second line we expand $\exp \left[2{\gamma }_{a}t\right]\simeq (1+2{\gamma }_{a}t)$. This is a good approximation because the linear damping of the parent mode is typically small and over the course of a ∼1 yr orbit the condition γ_a t ≪ 1 is very well satisfied (for the Jupiter model considered here, we have 1/γ_a ≃ 10⁹ yr or a quality factor QF ∼ ω_a /γ_a ≃ 10¹³; see also Goldreich & Nicholson 1977); it is well satisfied even if one uses the empirically determined value of QF ≃ 10⁵ based on the Jupiter−Io system (see, e.g., Lainey et al. 2009).

The total accumulated phase can be written as ${\phi }_{a}(t)\,=-\int \left[{\omega }_{a}^{{\prime} }+\delta {\omega }_{a}({\tilde{E}}_{a})\right]{dt}$. Of particular interest is the excess dephasing due to nonlinear interactions

$$\begin{eqnarray}&&\delta {\phi }_{\mathrm{nl}}=-\int \delta {\omega }_{a}\left[{\tilde{E}}_{a}(t)\right]{dt}=-\int \displaystyle \frac{\delta {\omega }_{a}({\tilde{E}}_{a})}{d{\tilde{E}}_{a}/{dt}}d{\tilde{E}}_{a},\end{eqnarray} \tag{ 40 }$$

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 ${\tilde{E}}_{a}$. If we use Equation (34) for δω_a and Equation (38) for $d{\tilde{E}}_{a}/{dt}$, then the nonlinear dephasing is

$$\begin{eqnarray}&&\delta {\phi }_{\mathrm{nl}}({\tilde{E}}_{a})=\displaystyle \frac{{\rm{\Omega }}}{2{\rm{\Gamma }}}\mathrm{ln}\left[\displaystyle \frac{{\gamma }_{a}+{\rm{\Gamma }}{\tilde{E}}_{a}}{{\gamma }_{a}+{\rm{\Gamma }}{\tilde{E}}^{(0)}}\right].\end{eqnarray} \tag{ 41 }$$

We can also write the (leading-order) dephasing as a function of time by plugging in Equation (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

$$\begin{eqnarray}&&\delta {\phi }_{\mathrm{nl}}(t)\simeq -\displaystyle \frac{{\rm{\Omega }}}{2{\rm{\Gamma }}}\mathrm{ln}\left[1+2{\rm{\Gamma }}{\tilde{E}}_{a}^{(0)}t\right],\end{eqnarray} \tag{ 42 }$$

$$\begin{eqnarray}&&\simeq \ -{\rm{\Omega }}{\tilde{E}}_{a}^{(0)}t,\end{eqnarray} \tag{ 43 }$$

where in the second equality we further assumed $2{{\rm{\Gamma }}}_{a}{\tilde{E}}_{a}^{(0)}t\ll 1$. For ${\tilde{E}}_{a}^{(0)}\simeq {10}^{-3}$ and t ≃ 1 yr, this condition is satisfied if Γ/ω_a < 10⁻⁵. As we will see shortly in the following section and Table 1, Equation (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 Equation (42) should be used instead. Obviously, Equation (43) also applies for conservative systems.

From the discussion above, we see that the leading-order 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 ${\tilde{E}}_{a}$ (Equations (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.

Using MESA (version 10398; Paxton et al. 2011, 2013, 2015, 2018, 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⁶ erg cm⁻² s⁻¹, 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 Section 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 resonance is not particularly strong. ¹⁴ Instead, upon summing over all couplings, the value of Ω is dominated by the coupling to the many m_a = 0 modes ("mode c" in Section 3.1), as well as the off-resonant m_a = −4 modes ("mode b" with ∣ω_b ∣ > 2∣ω_a ∣), and is therefore negative.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

By contrast, only the m_a = −4 daughters (oscillating at −2ω_prnt) contribute to the damping, and Γ is always positive (see both Section 3.1). For the model we consider here, Γ has a particularly small value of Γ/ω_prnt ∼ 10⁻⁹. 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 irradiation.

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 depend on different values of (Ω, Γ) and how (Ω, Γ) further vary with respect to (M, R) is deferred to future work.

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)

$$\begin{eqnarray}&&| {\rm{\Delta }}{\phi }_{k}| =| {\phi }_{k}-{\phi }_{k-1}| \gt { \mathcal O }(1)\mathrm{rad},\end{eqnarray} \tag{ 45 }$$

where ϕ_k is given by Equation (44). In linear theory, this is achieved through the tidal back-reaction on the orbit. At each passage, a random amount of energy ${\rm{\Delta }}{\tilde{E}}_{k}$ is removed from the orbit, which changes the orbital period by ΔP_orb,k (Equation (19)) and consequently the phase by ${\rm{\Delta }}{\phi }_{\mathrm{br},k}=-{\omega }_{a}^{{\prime} }{\rm{\Delta }}{P}_{\mathrm{orb},k}$, where the subscript "br" stands for "back-reaction." However, linear theory neglects the fact that the energy ${\rm{\Delta }}{\tilde{E}}_{k}$ gained by the planet's f-mode also changes its eigenfrequency by ${\rm{\Delta }}{\omega }_{a}({\rm{\Delta }}{\tilde{E}}_{k})\,=\delta {\omega }_{a}({\tilde{E}}_{k})-\delta {\omega }_{a}({\tilde{E}}_{k-1})\simeq {\rm{\Omega }}{\rm{\Delta }}{\tilde{E}}_{k}$. This frequency shift induces an additional random phase variation (relative to the previous orbit) ${\rm{\Delta }}{\phi }_{\mathrm{nl},k}\simeq -\left({\rm{\Omega }}{\rm{\Delta }}{\tilde{E}}_{k}\right){P}_{\mathrm{orb},k}$. The 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 ${\rm{\Delta }}{\tilde{E}}_{1}\equiv | {\rm{\Delta }}{q}_{1}{| }^{2}$ be the energy gained by the mode after the first pericenter passage (suppose it starts with an amplitude ∣q₀∣ ≪ ∣Δq₁∣). The phase shift caused by the tidal back-reaction after the first pericenter passage is thus ¹⁸

$$\begin{eqnarray}&&{\rm{\Delta }}{\phi }_{\mathrm{br},1}=-{\omega }_{a}^{{\prime} }{\rm{\Delta }}{P}_{\mathrm{orb},1}=\displaystyle \frac{3}{2}{\omega }_{a}^{{\prime} }{P}_{\mathrm{orb},0}\displaystyle \frac{{\rm{\Delta }}{\tilde{E}}_{1}}{\left|{\tilde{E}}_{\mathrm{orb},0}\right|},\end{eqnarray} \tag{ 46 }$$

where in the second equality we have used the fact that ${\tilde{E}}_{\mathrm{orb},0}={E}_{\mathrm{orb},0}/{E}_{0}\lt 0$. The threshold for growth is approximately ∣Δϕ_br,1∣ ≃ 1 rad (Vick & Lai 2018), which corresponds to a threshold one-kick energy ${\rm{\Delta }}{\tilde{E}}_{1}$ of

$$\begin{eqnarray}\begin{array}{rcl}{\rm{\Delta }}{\tilde{E}}_{\mathrm{br},1} & \simeq & 1.0\times {10}^{-5}{\left(\displaystyle \frac{{\omega }_{a}^{{\prime} }}{{\omega }_{a}}\right)}^{-1}{\left(\displaystyle \frac{{\omega }_{a}}{1.1{\omega }_{0}}\right)}^{-1}{\left(\displaystyle \frac{R}{1.1{R}_{{\rm{J}}}}\right)}^{5/2}\\ & & \times {\left(\displaystyle \frac{M}{{M}_{{\rm{J}}}}\right)}^{-3/2}{\left(\displaystyle \frac{{M}_{* }}{{M}_{\odot }}\right)}^{3/2}{\left(\displaystyle \frac{{a}_{\mathrm{orb},0}}{\mathrm{au}}\right)}^{-5/2}.\end{array}\end{eqnarray} \tag{ 47 }$$

The nonlinear frequency shift also leads to an excess phase of

$$\begin{eqnarray}&&{\rm{\Delta }}{\phi }_{\mathrm{nl},1}\simeq -\left({\rm{\Omega }}{\rm{\Delta }}{\tilde{E}}_{1}\right){P}_{\mathrm{orb},0},\end{eqnarray} \tag{ 48 }$$

where the subscript "nl" stands for "nonlinear" effects, and we have used Equation (34) for the nonlinear frequency shift. By setting ∣Δϕ_nl,1∣ = 1 rad, we similarly obtain the one-kick energy threshold to trigger diffusive growth solely from nonlinear effects,

$$\begin{eqnarray}\begin{array}{rcl}{\rm{\Delta }}{\tilde{E}}_{\mathrm{nl},1} & \simeq & 1.8\times {10}^{-6}{\left(\displaystyle \frac{| {\rm{\Omega }}| }{30{\omega }_{a}}\right)}^{-1}{\left(\displaystyle \frac{{\omega }_{a}}{1.1{\omega }_{0}}\right)}^{-1}{\left(\displaystyle \frac{R}{1.1{R}_{{\rm{J}}}}\right)}^{3/2}\\ & & \times {\left(\displaystyle \frac{M}{{M}_{{\rm{J}}}}\right)}^{-1/2}{\left(\displaystyle \frac{{M}_{* }}{{M}_{\odot }}\right)}^{1/2}{\left(\displaystyle \frac{{a}_{\mathrm{orb},0}}{\mathrm{au}}\right)}^{-3/2},\end{array}\end{eqnarray} \tag{ 49 }$$

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 Section 3.3, 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 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

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{| {\rm{\Delta }}{\phi }_{\mathrm{nl},1}| }{| {\rm{\Delta }}{\phi }_{\mathrm{br},1}| } & = & \displaystyle \frac{2}{3}\displaystyle \frac{| {\rm{\Omega }}| }{{\omega }_{a}}\displaystyle \frac{{\omega }_{a}}{{\omega }_{a}^{{\prime} }}\displaystyle \frac{| {E}_{\mathrm{orb},0}| }{{E}_{0}},\\ & = & 5.6\left(\displaystyle \frac{| {\rm{\Omega }}| }{30{\omega }_{a}}\right){\left(\displaystyle \frac{{\omega }_{a}^{{\prime} }}{{\omega }_{a}}\right)}^{-1}\left(\displaystyle \frac{R}{1.1{R}_{{\rm{J}}}}\right)\\ & & \times {\left(\displaystyle \frac{M}{{M}_{{\rm{J}}}}\right)}^{-1}\left(\displaystyle \frac{{M}_{* }}{{M}_{\odot }}\right){\left(\displaystyle \frac{{a}_{\mathrm{orb},0}}{\mathrm{au}}\right)}^{-1}.\end{array}\end{eqnarray} \tag{ 50 }$$

Note that the ratio is independent of ${\rm{\Delta }}{\tilde{E}}_{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.

5.1.1. Significance of Nonlinear Effects in Other Types of Eccentric Binaries

In Figure 3 we show a representative example of the first 300 orbits of a Jovian planet orbiting a solar-type star at a semimajor axis a_orb = 1 au using the iterative map described in Section 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⁻³ 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 ${\rm{\Delta }}{\tilde{E}}_{1}=| {\rm{\Delta }}{q}_{1}{| }^{2}=4.5\times {10}^{-6}$ (Equation (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 (Equation (45)). We see that in the linear case (gray 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 (Equation (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{\rm{\Delta }}{\tilde{E}}_{1}$, as one would expect for a diffusive process (i.e., the amplitude grows as the square root of the number of pericenter kicks).

Zoom In Zoom Out Reset image size

In Figure 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 yr) 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 ${\rm{\Delta }}{\tilde{E}}_{1}$ is shown in the top x-axis for each choice of D_peri. In the left panel, we assume that the planet is nonrotating, while in the right panel we assume that it is pseudo-synchronized with the orbit with Ω_s/ω_a = 0.19 (D_peri/0.02 au)^−3/2 and ${\omega }_{a}^{{\prime} }={\omega }_{a}+2{{\rm{\Omega }}}_{{\rm{s}}}$. If a mode experiences diffusive growth, then its energy after k pericenter passages is expected to be ${\tilde{E}}_{{\rm{k}}}\sim k{\rm{\Delta }}{\tilde{E}}_{1}$ on average. Since we set k = 500, we expect $\max \left[{\tilde{E}}_{k}\right]\,\simeq 500{\rm{\Delta }}{\tilde{E}}_{1}$ for a mode that grows diffusively (indicated by the black lines), while $\max \left[{\tilde{E}}_{k}\right]\sim {\rm{\Delta }}{\tilde{E}}_{1}$ for a mode that does not grow (assuming that the mode is off resonance with the orbit; in the right panel ${\omega }_{a}^{{\prime} }$ 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 Section 5.3).

Zoom In Zoom Out Reset image size

We see that each panel in Figure 4 can be divided up into two regions according to the maximum mode energy achieved. Let us first focus on the left panel (a nonrotating planet with ω_a P_orb,0/2π = 2818.73 being a noninteger). 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 ${\rm{\Delta }}{\tilde{E}}_{1}\simeq 1.4\times {10}^{-5}$. The analytical estimate (Equation (47); see vertical gray line) agrees well with the numerical results but slightly overestimates the threshold value of D_peri (and underestimates ${\rm{\Delta }}{\tilde{E}}_{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 performed over a full bin width of 0.1D_t). The error bars around the vertical lines are obtained by setting Δϕ₁ = 0.5 and 1.5 rad and then reevaluating ${\rm{\Delta }}{\tilde{E}}_{1}$ using Equations (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 ${\rm{\Delta }}{\tilde{E}}_{1}\simeq 2.4\times {10}^{-6}$, which is a factor of about 6 smaller than the linear case (see the vertical olive line and Equation (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 Figure 4). Because the one-kick energy depends sensitively on the pericenter distance, a factor of six change in ${\rm{\Delta }}{\tilde{E}}_{1}$ corresponds to a ≃10% increase in D_peri.

We consider the effects of planet spin in the right panel of Figure 4. We assume that the planet is rotating at a rate determined by the pseudo-synchronization condition (Equation (23)). ²⁰ As we show in Section 2.1, the mapping equations including spin are formally the same as the nonspinning equations except that the mode frequency is replaced by the inertial-frame value ${\omega }_{a}^{{\prime} }={\omega }_{a}+{m}_{a}{{\rm{\Omega }}}_{{\rm{s}}}$ (this frequency then enters the calculations of the one-kick amplitude and the linear propagation phase). Since we focus on a prograde mode with m_a = 2, ∣Δq₁∣ ∝ K₂₂ decreases sharply as ${\omega }_{a}^{{\prime} }$ increases (Equation (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 ${\omega }_{a}^{{\prime} }$ now varies with D_peri owing to the 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 ${\tilde{E}}_{k}^{(1)}\simeq {\rm{a}}\ \mathrm{few}\times {10}^{-5}\gg {\rm{\Delta }}{\tilde{E}}_{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 ${\omega }_{a}^{{\prime} }$ and P_orb, and this is why the upper envelope of the olive circles is at a lower value than that of the gray ones.

At 3.2 ≲ D_peri/D_t ≲ 3.8, we see that the chance resonance with the orbit may occasionally help a slightly subthreshold mode to also evolve into the diffusive regime. Suppose that the chance resonance helps the mode to initially build up an energy ${\tilde{E}}_{k,0}\simeq {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 $\sqrt{{\tilde{E}}_{k,0}{\rm{\Delta }}{\tilde{E}}_{1}}\gg {\rm{\Delta }}{\tilde{E}}_{1}$ (Wu 2018). If we replace ${\rm{\Delta }}{\tilde{E}}_{1}$ by $\sqrt{{\tilde{E}}_{k,0}{\rm{\Delta }}{\tilde{E}}_{1}}$ in Equations (46) and (48), we see that the new threshold one-kick energy becomes ${\rm{\Delta }}{\tilde{E}}_{{br},1}\simeq {10}^{-6}$ and ${\rm{\Delta }}{\tilde{E}}_{{nl},1}\simeq 3\times {10}^{-7}$ for modes that are initially in close resonance with the orbit. As D_peridecreases and ${\rm{\Delta }}{\tilde{E}}_{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 bottom panel. Eventually, when ${\rm{\Delta }}{\tilde{E}}_{1}$ reaches the value derived in Equations (47) and (49), almost all of the systems will grow diffusively.

An alternative way to consider the problem is to hold D_peri and thus ${\rm{\Delta }}{\tilde{E}}_{1}$ fixed and instead vary the initial semimajor axis a_orb,0 = D_peri/(1 − e_orb,0). By setting ∣Δϕ_br,1∣ = 1 rad in Equation (46) as before (Section 5.1) but now solving for ${a}_{\mathrm{orb},0}^{(\mathrm{br})}$, we find that the threshold to trigger diffusive growth due to only tidal back-reaction is

$$\begin{eqnarray}\begin{array}{rcl}{a}_{\mathrm{orb},0}^{(\mathrm{br})} & \simeq & 2.5\,\mathrm{au}{\left(\displaystyle \frac{{\rm{\Delta }}{\tilde{E}}_{1}}{{10}^{-6}}\right)}^{-2/5}{\left(\displaystyle \frac{{\omega }_{a}^{{\prime} }}{{\omega }_{a}}\right)}^{-2/5}{\left(\displaystyle \frac{{\omega }_{a}}{1.1{\omega }_{0}}\right)}^{-2/5}\\ & & \times \left(\displaystyle \frac{R}{1.1{R}_{{\rm{J}}}}\right){\left(\displaystyle \frac{M}{{M}_{{\rm{J}}}}\right)}^{-3/5}{\left(\displaystyle \frac{{M}_{* }}{{M}_{\odot }}\right)}^{3/5}.\end{array}\end{eqnarray} \tag{ 51 }$$

Similarly, we can use Equation (48) to find the threshold due to nonlinear mode interactions

$$\begin{eqnarray}\begin{array}{rcl}{a}_{\mathrm{orb},0}^{(\mathrm{nl})} & \simeq & 1.5\,\mathrm{au}{\left(\displaystyle \frac{{\rm{\Delta }}{\tilde{E}}_{1}}{{10}^{-6}}\right)}^{-2/3}{\left(\displaystyle \frac{| {\rm{\Omega }}| }{30{\omega }_{a}}\right)}^{-2/3}{\left(\displaystyle \frac{{\omega }_{a}}{1.1{\omega }_{0}}\right)}^{-2/3}\\ & & \times \left(\displaystyle \frac{R}{1.1{R}_{{\rm{J}}}}\right){\left(\displaystyle \frac{M}{{M}_{{\rm{J}}}}\right)}^{-1/3}{\left(\displaystyle \frac{{M}_{* }}{{M}_{\odot }}\right)}^{1/3}.\end{array}\end{eqnarray} \tag{ 52 }$$

In both cases, the threshold a_orb,0 increases with decreasing one-kick energy ${\rm{\Delta }}{\tilde{E}}_{1}$ (i.e., increasing D_peri). This is because both Δϕ_br and ${\rm{\Delta }}{\phi }_{\mathrm{nl}}\propto {P}_{\mathrm{orb}}{\rm{\Delta }}{\tilde{E}}_{1}$ (Equations (46) and (48)), and a longer ${P}_{\mathrm{orb}}\propto {a}_{\mathrm{orb}}^{3/2}$ is thus required in order for the f-mode to accumulate an excess phase ∣Δϕ_k ∣ to ${ \mathcal O }(1)\mathrm{rad}$. Additionally, Equations (51) and (52) scale differently with ${\rm{\Delta }}{\tilde{E}}_{1}$ because Δϕ_br ∝ ∣E_orb∣⁻¹, 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 Equation (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 Figure 4, Equations (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 Equations (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 Figure 5. Here we fix the pericenter distance to be D_peri = 0.02 au, corresponding to a one-kick energy of ${\rm{\Delta }}{\tilde{E}}_{1}\simeq 1.4\times {10}^{-5}$ (${\rm{\Delta }}{\tilde{E}}_{1}\simeq 6.7\times {10}^{-7}$) for a nonspinning (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 ${\rm{\Delta }}{\tilde{E}}_{1}\lesssim {10}^{-6}$, only planets born ≳2 au away from the host star can trigger diffusive tidal evolution and form hot Jupiters if only the linear 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/regular Jupiter occurrence rate (see Dawson & Johnson 2018).

Zoom In Zoom Out Reset image size

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

$$\begin{eqnarray}\begin{array}{rcl}{\left({\rm{\Delta }}{q}_{b,1}\right)}_{{U}_{a}^{2}} & = & \displaystyle \int i{\omega }_{b}{\kappa }_{b}{W}_{{l}_{a}{m}_{a}}^{2}\\ & & \times {Q}_{a}^{2}{\left(\displaystyle \frac{{M}_{* }}{M}\right)}^{2}{\left(\displaystyle \frac{R}{D}\right)}^{{l}_{b}+2}{e}^{i\left[\left({\omega }_{b}+{m}_{b}{{\rm{\Omega }}}_{{\rm{s}}}\right)\tau -{m}_{b}{\rm{\Phi }}\right]},\end{array}\end{eqnarray} \tag{ A5 }$$

where we used the fact that l_b = 2l_a and m_b = −2m_a . If we define

$$\begin{eqnarray}&&{K}_{{lm}}^{{\prime} }(\omega )=\displaystyle \frac{{\omega }_{0}}{2\pi }\int {\left[\displaystyle \frac{{D}_{\mathrm{peri}}}{D(\tau )}\right]}^{l+1}{e}^{i\left[\omega \tau -m{\rm{\Phi }}(\tau )\right]}d\tau ,\end{eqnarray} \tag{ A6 }$$

as a modified temporal overlap, then we can write the one-kick amplitude of the daughter mode as

$$\begin{eqnarray}\begin{array}{rcl}{\left({\rm{\Delta }}{q}_{b,1}\right)}_{{U}_{a}^{2}} & = & i2\pi {\kappa }_{b}{W}_{22}^{2}{Q}_{a}^{2}\left(\displaystyle \frac{{\omega }_{b}}{{\omega }_{0}}\right){\left(\displaystyle \frac{{M}_{* }}{M}\right)}^{2}\\ & & \times {\left(\displaystyle \frac{R}{{D}_{\mathrm{peri}}}\right)}^{6}{K}_{5,-4}^{{\prime} }({\omega }_{b}-4{{\rm{\Omega }}}_{{\rm{s}}}),\end{array}\end{eqnarray} \tag{ A7 }$$

where we have plugged in l_a = m_a = 2 for the parent and l_b = −m_b = 4 for the daughter. For typical values (a_orb = 1 au, e_orb = 0.98, ω_a = 5.6 × 10⁻⁴ rad s⁻¹, and ω_b = −1.9ω_a ), we find $\left|{\left({\rm{\Delta }}{q}_{b}\right)}_{{U}_{a}^{2}}\right|\simeq 0.4| {\rm{\Delta }}{q}_{a,1}^{2}|$.

Although initially ${\left({\rm{\Delta }}{q}_{b,1}\right)}_{{U}_{a}^{2}}$ is comparable to the steady-state solution q_b,ss (Equation (A4)), as the system starts to grow diffusively, its effect will soon become subdominant. This can be seen by noticing that even if ${\left({\rm{\Delta }}{q}_{b,1}\right)}_{{U}_{a}^{2}}$ can grow diffusively itself (e.g., due to a random P_orb), after k pericenter passages it only increases the amplitude of q_b by $\sqrt{k}| {\left({\rm{\Delta }}{q}_{b,1}\right)}_{{U}_{a}^{2}}| \propto {k}^{1/2}$ on average. On the other hand, ${q}_{b,\mathrm{ss}}\sim k{\left|{q}_{a,\mathrm{dyn}}\right|}^{2}\propto k$ because each ∣q_a,dyn∣ grows as $\sqrt{k}$. Consequently, the significance of the q_b,ss term increases as $\sqrt{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 ${\left({U}_{a}^{* }\right)}^{2}$ term can be ignored. It is also worth noting that the drive from ${\left({U}_{a}^{* }\right)}^{2}$ is stronger than the direct tidal force on mode b, ${U}_{b}\propto ({M}_{* }/M){\left(R/{D}_{\mathrm{peri}}\right)}^{{l}_{b}+1}$, by a factor of $({M}_{* }/M)(R/{D}_{\mathrm{peri}})\left({\kappa }_{b}{W}_{22}^{2}{Q}_{a}^{2}/{W}_{44}{Q}_{b}\right)\simeq 6\left({\kappa }_{b}{W}_{22}^{2}{Q}_{a}^{2}/{W}_{44}{Q}_{b}\right)$. It thus justifies why we can also ignore the daughter modes' linear coupling to the tide.

We now consider the effect of the history term on the daughter.

If we define ${c}_{b}={q}_{b}\exp (-2i{\omega }_{a}t)$, then by Equation (A3)

$$\begin{eqnarray}&&{\dot{c}}_{b}+\left(i{{\rm{\Delta }}}_{b}+{\gamma }_{b}\right){c}_{b}=i{\omega }_{b}{V}_{b},\end{eqnarray} \tag{ A8 }$$

where ${V}_{b}\equiv {\kappa }_{b}{q}_{a}^{* }{q}_{a}^{* }\exp (-2i{\omega }_{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}\sim {e}^{-i{\omega }_{a}t}$, and thus V_b is a constant. Near pericenter, however, q_a has an additional time dependence owing 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

$$\begin{eqnarray}\begin{array}{rcl}{c}_{b}(t) & = & {e}^{-\left(i{{\rm{\Delta }}}_{b}+{\gamma }_{b}\right)t}{\displaystyle \int }_{{t}_{0}}^{t}i{\omega }_{b}{V}_{b}(\tau ){e}^{\left(i{{\rm{\Delta }}}_{b}+{\gamma }_{b}\right)\tau }d\tau .\\ & = & {\left.\displaystyle \frac{{\omega }_{b}}{{{\rm{\Delta }}}_{b}-i\gamma }{V}_{b}\right|}_{{t}_{0}}^{t}-{e}^{-\left(i{{\rm{\Delta }}}_{b}+{\gamma }_{b}\right)t}{\omega }_{b}\\ & & \times {\displaystyle \int }_{{t}_{0}}^{t}\displaystyle \frac{{\dot{V}}_{b}(\tau )}{{{\rm{\Delta }}}_{b}-i{\gamma }_{b}}{e}^{\left(i{{\rm{\Delta }}}_{b}+{\gamma }_{b}\right)\tau }d\tau ,\end{array}\end{eqnarray} \tag{ A9 }$$

where the initial time is t₀ and we performed integration by parts to get the second line. For future convenience we set V_b (t₀) = 0 and thus drop the initial condition. Note that the first term in Equation (A9) recovers the steady-state solution, Equation (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}\sim {e}^{-i{\omega }_{a}t}$ and ${\dot{V}}_{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 have ${\dot{V}}_{b}\sim {{\rm{\Omega }}}_{\mathrm{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

$$\begin{eqnarray}\begin{array}{rcl}{c}_{b}(t) & \simeq & \displaystyle \frac{{\omega }_{b}}{{{\rm{\Delta }}}_{b}}{V}_{b}(t)+i\displaystyle \frac{{\omega }_{b}}{{{\rm{\Delta }}}_{b}}\displaystyle \frac{{\dot{V}}_{b}(t)}{{{\rm{\Delta }}}_{b}}\\ & & +{{ie}}^{-i{{\rm{\Delta }}}_{b}t}{\omega }_{b}{\displaystyle \int }^{t}\displaystyle \frac{{\ddot{V}}_{b}(\tau )}{{{\rm{\Delta }}}_{b}^{2}}{e}^{i{{\rm{\Delta }}}_{b}\tau }d\tau =...,\end{array}\end{eqnarray} \tag{ A10 }$$

where we dropped γ_b to simplify the notation. Note that after the nth iteration of integration by parts, we have a correction that depends only on the instantaneous value $\propto {V}_{b}^{(n-1)}/{{\rm{\Delta }}}_{b}^{n-1}$ and an integrand $\propto {V}_{b}^{(n)}/{{\rm{\Delta }}}_{b}^{n}\sim {\left({{\rm{\Omega }}}_{\mathrm{peri}}/{{\rm{\Delta }}}_{b}\right)}^{n}{V}_{b}$, where ${V}_{b}^{(n)}$ is the nth time derivative of V_b . Since ∣Δ_b ∣ ≫ Ω_peri, the history-dependent term gets progressively smaller with increasing n, and the instantaneous corrections to Equation (29) form a converging series (in fact, the corrections are nonzero only around a pericenter passage). This is analogous to the fact that the linear tide can be well approximated by its instantaneous equilibrium component when the tidal forcing frequency is much smaller than the mode frequency.

For the most resonant daughter mode b (i.e., the l_b = −m_b = 4, ω_b < 0 f-mode), it is possible to have ∣Δ_b ∣ < Ω_peri. In this case, the series expansion formed by integration by parts does not converge. Instead, we need to directly solve Equation (A9). To do so, we consider the following simple model of V_b near the kth pericenter passage (corresponding to time t = t_k ):

$$\begin{eqnarray}\,\begin{array}{l}{V}_{b}(t)=\left\{\begin{array}{ll}-{\rm{\Delta }}{V}_{{bk}}, & \mathrm{if}\ t\lt {t}_{k}-\displaystyle \frac{\pi }{\eta {{\rm{\Omega }}}_{{\rm{p}}}},\\ {\rm{\Delta }}{V}_{{bk}}\sin \left[\eta {{\rm{\Omega }}}_{{\rm{p}}}(t-{t}_{k})\right], & \mathrm{if}\ {t}_{k}-\displaystyle \frac{\pi }{\eta {{\rm{\Omega }}}_{{\rm{p}}}}\leqslant t\leqslant {t}_{k}+\displaystyle \frac{\pi }{\eta {{\rm{\Omega }}}_{{\rm{p}}}},\\ {\rm{\Delta }}{V}_{{bk}}, & \mathrm{if}\ t\gt {t}_{k}+\displaystyle \frac{\pi }{\eta {{\rm{\Omega }}}_{{\rm{p}}}},\end{array}\right.\end{array}\,\end{eqnarray} \tag{ A11 }$$

where η ∼ 1 is a correction on the characteristic timescale over which q_a changes, and we rewrote Ω_peri as Ω_p in order to reduce notational clutter. With this simple model we can easily evaluate the integration around t_k as

$$\begin{eqnarray}&&\begin{array}{l}{\displaystyle \int }_{{t}_{k}-\pi /\eta {{\rm{\Omega }}}_{{\rm{p}}}}^{{t}_{k}+\pi /\eta {{\rm{\Omega }}}_{{\rm{p}}}}{\dot{V}}_{b}(\tau ){e}^{\left(i{{\rm{\Delta }}}_{b}+{\gamma }_{b}\right)\tau }d\tau \\ \ \simeq \ \left[\displaystyle \frac{\eta {{\rm{\Omega }}}_{{\rm{p}}}}{\left({{\rm{\Delta }}}_{b}+\eta {{\rm{\Omega }}}_{{\rm{p}}}-i{\gamma }_{b}\right)}\sin \left(\displaystyle \frac{{{\rm{\Delta }}}_{b}+\eta {{\rm{\Omega }}}_{{\rm{p}}}}{\eta {{\rm{\Omega }}}_{{\rm{p}}}}\pi \right)\right.\\ \ \ \left.+\ \displaystyle \frac{\eta {{\rm{\Omega }}}_{{\rm{p}}}}{\left({{\rm{\Delta }}}_{b}-\eta {{\rm{\Omega }}}_{{\rm{p}}}-i{\gamma }_{b}\right)}\sin \left(\displaystyle \frac{{{\rm{\Delta }}}_{b}-\eta {{\rm{\Omega }}}_{{\rm{p}}}}{\eta {{\rm{\Omega }}}_{{\rm{p}}}}\pi \right)\right]{\rm{\Delta }}{V}_{{bk}},\end{array}\end{eqnarray} \tag{ A12 }$$

where we ignored dissipation over a time 2π/ηΩ_p. We can now write mode b's amplitude as ²¹

$$\begin{eqnarray}\begin{array}{rcl}{c}_{b}(t) & \simeq & \displaystyle \frac{{\omega }_{b}}{{{\rm{\Delta }}}_{b}-i{\gamma }_{b}}\left\{{V}_{b}(t)+{e}^{-\left(i{{\rm{\Delta }}}_{b}+{\gamma }_{b}\right)t}\right.\\ & & \left.\times \left[\displaystyle \frac{\eta {{\rm{\Omega }}}_{{\rm{p}}}\sin \left(\tfrac{{{\rm{\Delta }}}_{b}\pm \eta {{\rm{\Omega }}}_{{\rm{p}}}}{\eta {{\rm{\Omega }}}_{{\rm{p}}}}\pi \right)}{\left({{\rm{\Delta }}}_{b}\pm \eta {{\rm{\Omega }}}_{{\rm{p}}}-i{\gamma }_{b}\right)}\right]\sum _{k}^{{t}_{k}\lt t}{\rm{\Delta }}{V}_{{bk}}\right\}.\end{array}\end{eqnarray} \tag{ A13 }$$

Therefore, in addition to the instantaneous term V_b (t), there in principle should also be a history-dependent term ∑ΔV_bk . Physically, this case can be understood by the following. As the parent mode's amplitude changes at each pericenter passage over a timescale 1/ηΩ_peri, its frequency content is broadened from a single delta function at ω_a to a band covering ω_a ± ηΩ_peri/2. Since the parent–daughter detuning is small, ∣Δ_b ∣ < Ω_peri, the broadened drive from the parent can now resonantly excite the daughter. This thus gives the daughter a "dynamical" component that depends on the past history (∑ΔV_bk ). By contrast, when the daughter is not resonant, the instantaneous "equilibrium" component dominates.

If ∑ΔV_bk grows diffusively (since ${\rm{\Delta }}{V}_{{bk}}\propto {q}_{a}^{* }{\rm{\Delta }}{q}_{a,1}$, and the phase of ${q}_{a}^{* }$ can be random), then both V_b and ∑ΔV_bk grow with the number of pericenter passages k as ∝k. Therefore, the history term due to the kick on q_a , unlike the one due to U_a , can be potentially important. For simplicity, we drop it in the analysis of this paper and defer its consideration to future work.