Analytical models for secular descents in hierarchical triple systems

Abstract

Triple body systems are prevalent in nature, from planetary to stellar to supermassive black hole scales. In a hierarchical triple system, oscillations of the inner orbit’s eccentricity and inclination can be induced on secular timescales. Over many cycles, the octupole-level terms in the secular equations of motion can drive the system to extremely high eccentricities via the Eccentric Kozai-Lidov (EKL) mechanism. The overall decrease in the inner orbit’s pericenter distance has potentially dramatic effects for realistic systems, such as tidal disruption. We present an analytical approximation in the test particle limit to describe the step-wise eccentricity evolution of the inner orbit. We also integrate the equations of motion and present approximations in both the test particle limit and in the general case to describe the overall octupole-level time evolution of the eccentricity. The analytical approximations are compared to numerical simulations to show that the models accurately describe the form and timescale of the secular descent from large distances to a close-encounter distance (e.g., the Roche limit). By circumventing the need for numerical simulations to obtain the long-term behavior, these approximations can be used to readily estimate descent timescales for populations of systems. We demonstrate by calculating rates of EKL-driven migration for Hot Jupiters in stellar binaries.

1 Introduction

Hierarchical triple systems have been studied extensively in the literature in recent years, with a wide range of applications, from exoplanets to gravitational wave emitting black holes (for a review, see Naoz, 2016). A hierarchical system is defined by having a relatively tight “inner” binary with semi-major axis $a_{1}$ orbited by a faraway “outer” companion with semi-major axis $a_{2}$. In this case, the secular approximation can be applied, where the three-body Hamiltonian can be averaged over the orbital periods and expanded in powers of the small semi-major axis ratio $\alpha=a_{1}/a_{2}$ (e.g., Kozai, 1962; Harrington, 1968; Ford et al., 2000; Naoz et al., 2013a). At the quadrupole level (proportional to $\alpha^{2}$), the inner orbit can undergo oscillations of eccentricity and inclination on timescales much longer than the orbital periods (Kozai, 1962; Lidov, 1962). The quadrupole level of approximation is often insufficient when the outer companion’s orbit is eccentric, or when members have non-negligible mass (Naoz et al., 2013a). In these cases, the octupole contribution (proportional to $\alpha^{3}$) can cause the inner orbit to reach extremely high eccentricities or even flip from prograde to retrograde with respect to the total angular momentum, a phenomenon known as the Eccentric Kozai-Lidov (EKL) mechanism (Naoz, 2016).

These secular perturbations can lead the system to reach extreme eccentricities over a wide range of the parameter space (e.g., Teyssandier et al., 2013; Li et al., 2014a, b). Specifically, Li et al. (2014a) showed that the parameter space can be roughly divided into a chaotic regime and a nearly regular regime. In the latter, the long-term evolution of the quadrupolar cycles is modulated at the octupole level (e.g., Naoz et al., 2011; Katz et al., 2011; Lithwick & Naoz, 2011). The octupolar envelope of the oscillations traces successive maxima of eccentricity reached during individual cycles. As the eccentricity increases, the pericenter distance of the inner orbiting body decreases. Therefore, the overall evolution of the inner orbiting body can be described as a secular descent driven by the octupole-level terms.

Secular descents may play an important role in shaping planetary systems. In a system consisting of an inner binary with a star and Jupiter-like planet orbited by a distant planetary or stellar companion, secular perturbations from the distant body can excite the inner planet’s orbital eccentricity, causing it to undergo close encounters with its host star (e.g., Wu & Murray, 2003; Fabrycky & Tremaine, 2007; Naoz et al., 2011, 2012; Li et al., 2014b; Petrovich, 2015a, b; Petrovich & Tremaine, 2016). When tidal interactions are considered, the inner planet’s semi-major axis can shrink and result in a Hot Jupiter (e.g., Mayor & Queloz, 1995; Naoz et al., 2011, 2012; Petrovich, 2015a, b; Dawson & Johnson, 2018). Furthermore, the planet can be tidally disrupted if its pericenter distance is sufficiently reduced to cross the Roche limit, an effect that is thought to shape the demographics of the Hot Jupiter population (e.g., Guillochon et al., 2011; Naoz et al., 2012; Valsecchi et al., 2015; Petrovich, 2015a; Muñoz et al., 2016).

The EKL mechanism may also drive secular descents in supermassive black hole (SMBH) binaries (e.g., Li et al., 2015; Chen et al., 2009, 2011). The outer SMBH can drive a star towards the inner SMBH to produce a tidal disruption event (TDE). EKL was shown to produce a TDE rate consistent with E+A (post-starburst) galaxies, as well as naturally produce partial, repeated TDEs (Melchor et al., 2023). Furthermore, this mechanism may be responsible for off-nuclear TDEs (Mockler et al., 2023) and may also act to produce extreme-mass ratio inspiral gravitational wave events (Naoz et al., 2022; Naoz & Haiman, 2023).

In each of these physical systems, it is necessary to understand the secular reduction of the inner orbit’s pericenter from large values to a close-encounter distance (e.g., the Roche limit). While the eccentricity evolution can be numerically solved given the initial orbital configuration of a system, an analytical approximation for the long-term eccentricity evolution would be useful to readily estimate the form and timescale of the secular descent. For a large population of systems, having an analytical approximation would circumvent the need for costly numerical integrations and allow for efficient predictions of migration rates and demographic features.

In this paper, we derive an analytical model that describes the secular descent to some minimum pericenter distance. Specifically, we approach the problem from three different directions that provide results consistent with numerical simulations. In §2, we review the physical picture of the EKL mechanism to motivate the models. We present the models for both the test particle case and general case in §3. Specifically, in §3.1.1 we employ a step-like approach to the changes in the z-axis angular momentum. Then, in §3.1.2 and 3.2, we use the equation of motions of the test particle and non-test particle, respectively, to derive the descent timescale. We describe the regimes of applicability for the models and quantify the accuracy of the models by comparing them to numerical simulations. Finally, §4 demonstrates the application of the models to the Hot Jupiter population. We summarize our conclusions in §5.

2 Physical Picture and Motivation

Here, we review the physical picture of the Eccentric Kozai-Lidov Mechanism (see, e.g., Naoz, 2016). The hierarchical three-body system consists of an inner binary (masses $m_{1}$ and $m_{2}$) and a distant third body (mass $m_{3}$). We define the semi-major axis, eccentricity, argument of periastron, and the longitude of ascending nodes of the inner (outer) orbit as $a_{1},e_{1},\omega_{1}$ and $\Omega_{1}$ ($a_{2},e_{2},\omega_{2}$ and $\Omega_{2}$). We choose the z-axis to be parallel to the total angular momentum; thus, the longitude of ascending nodes have the following relations: $\Omega_{1}-\Omega_{2}=\pi$ (e.g., Naoz et al., 2013a). The specific angular momentum is defined by $J_{1}=\sqrt{1-e_{1}^{2}}$ ($J_{2}=\sqrt{1-e_{2}^{2}}$) for the inner (outer) orbit. In this case, the inclination of the inner and outer orbits is defined by the angular momentum components projected on the z-axis. In other words:

$$J_{z,1} = J_{1}\cos i_{1}$$ (1)

$$J_{z,2} = J_{2}\cos i_{2}\ .$$ (2)

The mutual inclination is then simply $i=i_{1}+i_{2}$ (e.g., Ford et al., 2000).

In order to motivate the analytical approximation below, we begin by describing the test particle system, which is a 2-degree-of-freedom system and, thus, provides a logical first step to analyze the physical system (see also Lithwick & Naoz, 2011; Katz et al., 2011; Li et al., 2014a, b). Under the test particle approximation ($m_{2}\to 0$), the outer magnitude and organization of the angular momentum remain constant, where $i_{2}\to 0$ and $i_{1}\to i$, thus we define $\theta=\cos i_{1}$ and $J_{z}=J_{z,1}$. We also define $J=J_{1}$. In this case, the Hamiltonian can be written as (Naoz, 2016):

$$\mathcal{H}^{TP}=\frac{3}{8}k^{2}\frac{m_{1}m_{3}}{a_{2}}\left(\frac{a_{1}}{a_{2}}\right)^{2}\frac{1}{\left(1-e_{2}^{2}\right)^{3/2}}F\ ,$$ (3)

where $k^{2}$ is Newton’s constant, $F\equiv F_{\text{quad }}+\epsilon F_{\text{oct }}$ is the energy function, and

$$\epsilon=\frac{a_{1}}{a_{2}}\frac{e_{2}}{1-e_{2}^{2}}$$ (4)

characterizes the strength of the octupole component. The quadrupole component of $F$ is

$$F_{\text{quad }}=-\frac{e_{1}^{2}}{2}+\theta^{2}+\frac{3}{2}e_{1}^{2}\theta^{2}+\frac{5}{2}e_{1}^{2}\left(1-\theta^{2}\right)\cos\left(2\omega_{1}\right)\ ,$$ (5)

and the octupole component is

$$F_{\text{oct}}= \frac{5}{16}\left(e_{1}+\frac{3}{4}e_{1}^{3}\right)\left[\left(1+11\theta-5\theta^{2}-15\theta^{3}\right)\cos(\omega_{1}+\Omega_{1})\right.$$ (6)

$$\left.+\left(1-11\theta-5\theta^{2}+15\theta^{3}\right)\cos(\omega_{1}-\Omega_{1})\right]$$

$$-\frac{175}{64}e_{1}^{3}\left[\left(1-\theta-\theta^{2}+\theta^{3}\right)\cos(3\omega_{1}-\Omega_{1})\right.$$

$$\left.+\left(1+\theta-\theta^{2}-\theta^{3}\right)\cos(3\omega_{1}+\Omega_{1})\right]\ .$$

At the quadrupole level ($\epsilon=0$, corresponding to a circular outer orbit), $F$ and $J_{z}$ are the two constants of motion. There are two classes of trajectories, librating and circulating. Librating trajectories are associated with bound oscillations of $\omega_{1}$. The eccentricity evolves to large values on circulating trajectories, which have $e_{1}$ smallest and $i$ largest at $\omega_{1}=0$ and vice versa for $\omega_{1}=\pm\pi/2$. The separatrix separates these modes of behavior and has $e_{1,0}=0$ and $\theta=\sqrt{3/5}$. As such, large eccentricity oscillations are expected to occur between inclination angles of $\cos^{-1}\sqrt{3/5}\approx 39^{\circ}$ and $141^{\circ}$.

For circulating trajectories with $e_{1,0}\ll 1$ and inclinations between $39^{\circ}$ and $141^{\circ}$, the maximum eccentricity $e_{1,\text{max}}$ is given by (e.g., Lithwick & Naoz, 2011)

$$e_{1,\text{max}}^{2}\approx 1-\frac{5}{3}J_{z}^{2}\ ,$$ (7)

and the minimum eccentricity $e_{1,\text{min}}$ is given by

$$e_{1,\text{min}}^{2}\approx\frac{1}{2}(F-J_{z}^{2})\ .$$ (8)

Antognini (2015) estimated the timescale for each quadrupolar oscillation of $e_{1}$ to be

$$t_{\mathrm{quad}}\approx\frac{16}{15}\frac{a_{2}^{3}\left(1-e_{2}^{2}\right)^{3/2}\sqrt{m_{1}+m_{2}}}{a_{1}^{3/2}m_{3}k}\ .$$ (9)

For an eccentric outer body ($\epsilon>0$), $F$ is now the only conserved quantity, as the octupole component modulates $J_{z}$. Trajectories can be regular or chaotic (e.g., Lithwick & Naoz, 2011; Li et al., 2014b); in this work, we focus on regular trajectories. For these, $J_{z}$ changes in a step-wise fashion, swinging to a new value during each plunge in $e_{1}$. The overall evolution of $J_{z}$ controls the envelope of the $e_{1}$ maxima. We demonstrate this behavior in Figure 1. It is this envelope that we seek to model, as it captures the overall secular descent of the inner body’s pericenter distance $a_{1}(1-e_{1})$.

Figure 1: Time evolution of $1-e_{1}$ for trajectories with $e_{2}=0$ (orange) and $e_{2}=0.17$ (blue) integrated numerically. Increasing the outer orbit’s eccentricity drives $e_{1}$ to evolve to large values over multiple quadrupole cycles. $J_{z}$ controls the envelope of the eccentricity evolution, being constant in the $e_{2}=0$ case (black), and evolving in a step-wise fashion for $e_{2}>0$ (red).

Relaxing the test particle approximation can lead to differences in orbital evolution. The inner orbit can torque the outer orbit, making large eccentricity oscillations occur in different regions of the parameter space (e.g., Teyssandier et al., 2013; Naoz & Fabrycky, 2014; Martin & Triaud, 2015). Below, we provide an analytical description for systems extending beyond the test particle limit.

3 Analytical Approximations of the Secular Descent

Here, we pedagogically describe our analytical approximations of the secular descent. They are divided into two cases: test particle and general.

3.1 Test particle approach

As mentioned above, the test particle approach provides a 2-degree-of-freedom physical system. In the test particle approximation, $m_{2}\to 0$. The full set of equations can be found in Naoz (2016). Here, in Appendix A, we present the equations of motion for completeness. We find that the analytical description can be further divided into two cases: one to describe single steps in $J_{z}$ and the other to describe the overall secular descent of the pericenter distance as eccentricity increases.

3.1.1 Test particle approach, steps in eccentricity

We present an approach in the test particle limit to describe single steps in $J_{z}$. Here, we assume small initial $J_{z}$ ($J_{z,0}\lesssim 0.1$), corresponding to initially high mutual inclination ($i_{0}\gtrsim 85^{\circ})$ for an initially circular orbit. The individual steps can be pieced together to estimate the overall secular descent, as illustrated in Figure 1.

Since most of the change in $J_{z}$ occurs when $\dot{\omega}_{1}\sim 0$ (e.g., Li et al., 2014b), in the limit of small $J_{z}$, $\dot{\omega}_{1}\sim 0$ implies that (see Eq. (A9))

$$\cos{2\omega_{1}}=\frac{1}{5}\frac{5J_{z}^{2}-J^{4}}{J_{z}^{2}-J^{4}}\approx\frac{1}{5}\ .$$ (10)

We seek to characterize the change in $J_{z}$ during a single quadrupolar swing in $J$ (i.e. half a quadrupole cycle). We find the time evolution of $J$ from Eq. (A7), which for $\dot{\omega}_{1}\sim 0$, from Eq. (10), is

$$\frac{dJ}{dt}\bigg{|}_{\dot{\omega}_{1}\sim 0}\approx-2\sqrt{6}(1-J^{2})\left[1-\left(\frac{J_{z}}{J}\right)^{2}\right]\ .$$ (11)

Since $J_{z}$ at the quadrupole level is constant (for a test particle), we find the octupole level time evolution of $J_{z}$ from Eq. (A8), which for $\dot{\omega}_{1}\sim 0$ and taking the case $\Omega_{1}=0$ is

$$\frac{dJ_{z}}{dt}\bigg{|}_{\dot{\omega}_{1}\sim 0} \approx -\frac{5}{32}\sqrt{\frac{2}{5}}\epsilon\sqrt{1-J^{2}}\left(7-3J^{2}\right)\frac{J_{z}}{J}\left(11-15\frac{J_{z}^{2}}{J^{2}}\right)$$ (12)

$$+ \frac{245}{32}\sqrt{\frac{2}{5}}\epsilon\left(1-J^{2}\right)^{3/2}\frac{J_{z}}{J}\left(1-\frac{J_{z}^{2}}{J^{2}}\right)\ .$$

To characterize the step in $J_{z}$ during a quadrupole cycle, we divide $d{J_{z}}/dt$ by $dJ/dt$, which gives

$$\frac{dJ_{z}}{J_{z}}\approx\epsilon\frac{dJ}{J}\left[\vphantom{\frac{1}{5}}\right.\frac{5}{64\sqrt{15}}\frac{\left(7-3J^{2}\right)}{\sqrt{1-J^{2}}}\left(11-15\frac{J_{z}^{2}}{J^{2}}\right)\times$$

$$\left(1-\frac{J_{z}^{2}}{J^{2}}\right)^{-1}-\frac{245}{64\sqrt{15}}\sqrt{1-J^{2}}\left.\vphantom{\frac{1}{2}}\right]\ .$$ (13)

At face value, integrating both sides is challenging because $J_{z}$ appears both in the left and right-hand side of the equation. However, recall that over a quadrupole cycle, $J$ changes significantly, but $J_{z}$ remains approximately constant (Lithwick & Naoz, 2011). Thus, we treat $J_{z}$ as a constant and find the approximation

$$\ln J_{z} \approx \frac{35}{16\sqrt{15}}\epsilon\ln\left({\frac{1-\sqrt{1-J^{2}}}{{1+\sqrt{1-J^{2}}}}}\right)$$ (14)

$$+ \frac{5}{32\sqrt{15}}\epsilon\frac{(7-3J_{z}^{2})}{\sqrt{1-J_{z}^{2}}}\ln\left(\frac{-\sqrt{1-J^{2}}+\sqrt{1-J_{z}^{2}}}{\sqrt{1-J^{2}}+\sqrt{1-J_{z}^{2}}}\right)$$

$$- \frac{5}{4\sqrt{15}}\epsilon\sqrt{1-J^{2}}+\text{const. }\ .$$

From here we find the relation between the final ($J_{z,f}$) and initial ($J_{z,0}$) values of $J_{z}$ during a step in $J$ (half a quadrupole cycle) to be

$$\frac{J_{z,f}}{J_{z,0}} \approx \exp\left(-\frac{5}{4\sqrt{15}}\epsilon\sqrt{1-J^{2}}\right)$$ (15)

$$\times \left({\frac{1-\sqrt{1-J^{2}}}{{1+\sqrt{1-J^{2}}}}}\right)^{\zeta_{1}}$$

$$\times \left(\frac{-\sqrt{1-J^{2}}+\sqrt{1-J_{z}^{2}}}{\sqrt{1-J^{2}}+\sqrt{1-J_{z}^{2}}}\right)^{\zeta_{2}}\ ,$$

where

$$\zeta_{1}=\frac{35}{16\sqrt{15}}\epsilon\ ,$$ (16)

and

$$\zeta_{2}=\frac{5}{32\sqrt{15}}\epsilon\frac{(7-3J_{z}^{2})}{\sqrt{1-J_{z}^{2}}}\ .$$ (17)

We use this expression to approximate the step in $J_{z}$ from low to high eccentricity. For the step from high to low eccentricity, the system swings to the other $\dot{\omega_{1}}=0$ root, reversing the sign of $dJ/dt$.

The timescale $t_{\text{step}}$ between a local maximum and a minimum in $J$ can then be found by integrating Eq. (11) between the maximum eccentricity in Eq. (7) and the minimum eccentricity in Eq. (8), giving

$$t_{\text{step}}\approx\\ \int_{\sqrt{1-\frac{1}{2}(F-J_{z}^{2})}}^{\sqrt{\frac{5}{3}}J_{z}}\left(-2\sqrt{6}[1-J^{2}]\left[1-\left(\frac{J_{z}}{J}\right)^{2}\right]\right)^{-1}dJ\\ =-\frac{1}{4\sqrt{6}(J_{z}^{2}-1)}\left[\vphantom{\frac{1}{2}}\right.J_{z}\ln(|J+J_{z}|)-J_{z}\ln(|J-J_{z}|)\\ -\ln(|J+1|)+\ln(|J-1|)\left.\vphantom{\frac{1}{2}}\right]\bigg{|}_{\sqrt{1-\frac{1}{2}(F-J_{z}^{2})}}^{\sqrt{\frac{5}{3}}J_{z}}\ ,$$ (18)

where the conversion to true time is described in Eq. (A11). See the annotation on the right-hand side of Figure 2. We note that this timescale is similar to half the quadrupole timescale given in Eq. (9), but it is more accurate for our needs.

Figure 2: Comparison of the numerical evolution to the analytical approximation for the step-wise eccentricity evolution. Here we consider a test particle at 1 AU from a $1M_{\odot}$ star orbited by a distant perturber of mass $0.16M_{\odot}$ at 20 AU. The system is initially set with $e_{1}=0.01$, $e_{2}=0.17$, $\omega_{1}=0^{\circ}$, $\omega_{2}=0^{\circ}$, and $i=86.5^{\circ}$. Left: We compare one step in $J_{z}$ (orange) modeled by Eq. (15) to the numerical evolution (blue). We find that the evolution of $J_{z}$ closely agrees for much of the step. The lower accuracy at low $J$ is explained by $\omega_{1}$ varying here and the assumption of $\dot{\omega}_{1}=0$ breaking down. Right: We stitch together the individual steps (orange dots represent each step, the orange curve represents the overall secular descent) as described in §3.1.1. The zoom-in panel shows one step with timescale estimated from Eq. (18). Comparing to the numerical evolution (blue), we see that the overall time evolution of $J_{z}$ is closely captured.

The total timescale for a descent then is

$$t_{\rm descent,TP}=\sum_{J_{z,0}}^{J_{z,\rm min}}t_{\rm step}\ ,$$ (19)

where $J_{z}$ is updated on each step from Eq. (15). $J_{z,\rm min}$ can be set at a desired minimum pericenter distance (e.g., the Roche limit). In a realistic system, the octupole-level of approximation can drive the eccentricity to extreme values, which is often associated with an orbital flip (e.g., Naoz et al., 2011; Katz et al., 2011; Lithwick & Naoz, 2011; Li et al., 2014a; Teyssandier et al., 2013). We find that for systems that flip, the octupole level maximum eccentricity reached at the flip is roughly at the range of $1-e_{1}\sim 10^{-4}$ to $10^{-6}$, which agrees with values obtained in previous investigations (e.g., Katz et al., 2011; Li et al., 2014b; Liu et al., 2015). In this work, we specify the maximum eccentricity that we allow the system to evolve to.

In the left panel of Figure 2, we show the changes of $J_{z}$ as a function of $J$ for a representative system calculated numerically (Naoz et al., 2013a), overplotting the analytical model in Eq. (15). In the right panel, we present the time evolution of $1-e_{1}$ of the same system, overplotting the analytical secular descent as a function of time. To plot this line, we stitch together individual $J_{z}$ steps and estimate the timescale for each step from Eq. (18). The analytical model and numerical evolution agree to within a small factor. Therefore, this analytical approximation can be used to closely estimate the timescale to reach some eccentricity maximum, as well as the step-wise increase of the pericenter distance.

3.1.2 Test particle approach, full descent

Over many quadrupole cycles, the system oscillates between maximum and minimum eccentricity until an octupole level maximum for the eccentricity is reached. This maximum is often associated with a very small angular momentum, which has potentially dramatic effects on realistic systems, like planet consumption by a star or a tidal disruption event. Here, we analytically describe the time evolution of a test particle’s maximum eccentricity, $e_{1,\text{max}}$ until the octupole-level maximum eccentricity, i.e., a full secular decent. The analytical approach, in this case, is different than in §3.1.1. In this case, we adopt the test particle equation of motion and integrate to obtain a single analytical formula for the overall secular descent for systems with $i_{0}\gtrsim 45^{\circ}$. This approach is motivated by Katz et al. (2011), in which the equations of motion are averaged over the quadrupole cycles to yield approximate equations for the long-term octupolar eccentricity evolution.

Figure 3: Comparison of the numerical evolution (blue) with the analytical approximations to the secular descent for an example Jupiter-like planet. The gold-shaded region corresponds to the test particle case from Eq. (26), and the grey-shaded region corresponds to the general case from Eq. (34) (see §3.2). The spread in the shaded regions represents the $\sim$25% variations in the numerical factor $\Upsilon$ obtained from Monte Carlo simulations (see Appendix B), corresponding to varying evolutions of $\omega_{2}$. Here we consider a $1M_{\text{J}}$ planet at 1 AU from a $1M_{\odot}$ star orbited by a distant perturber of mass $0.1M_{\odot}$ at 20 AU. The system is initially set with $e_{1}=0.01$, $e_{2}=0.4$, $\omega_{1}=0^{\circ}$, $\omega_{2}=0^{\circ}$, and $i=70^{\circ}$. The approximations closely track the envelope of the eccentricity oscillations through the descent.

Here, we assume small initial eccentricity, $e_{1,0}\ll 1$. In this case, the maximum eccentricity for a single quadrupole cycle is approximated by Eq. (7). Differentiating the expression with respect to $t$, we obtain

$$\frac{d{e_{1,\text{max}}}}{dt}\approx-\frac{5}{3}\theta\frac{1}{e_{1,\text{max}}}\sqrt{1-e_{1,\text{max}}^{2}}\frac{dJ_{z}}{dt}\ .$$ (20)

The quadrupole component of ${dJ_{z}}/{dt}$ is zero, and the octupole component is obtained from Equations (6) and (A8). For trajectories in which $\omega_{1}$ is circulating, the eccentricity is largest at $\omega_{1}\sim\pm\pi/2$ and $\theta\sim\sqrt{3/5}$. Thus, plugging in these values, we obtain an expression for the time evolution of $e_{1,\text{max}}$:

$$\frac{d{e_{1,\text{max}}}}{dt}f({e_{1,\text{max}}}) \approx$$

$$\frac{15}{64}k\frac{a_{1}^{5/2}}{a_{2}^{4}}\frac{m_{3}}{\sqrt{m_{1}}}\frac{e_{2}}{(1-e_{2}^{2})^{5/2}}\cos{\omega_{2}} \ ,$$ (21)

where

$$f({e_{1,\text{max}}})=\frac{1}{(2+5e_{1,\text{max}}^{2})\sqrt{1-e_{1,\text{max}}^{2}}}\ .$$ (22)

We are interested in the pericenter distance of the inner body and thus adopt $r_{p}=1-e_{1,\text{max}}$. Therefore, we can write Eq. (3.1.2) as

$$\frac{d{r_{p}}}{dt}f(r_{p})\approx-\frac{15}{64}k\frac{a_{1}^{5/2}}{a_{2}^{4}}\frac{m_{3}}{\sqrt{m_{1}}}\frac{e_{2}}{(1-e_{2}^{2})^{5/2}}\cos{\omega_{2}}\ ,$$ (23)

where $f(r_{p})$ is expanded by taking $r_{p}\ll 1$

$$f(r_{p})=\frac{1}{(5r_{p}^{2}-10r_{p}+7)\sqrt{-r_{p}^{2}+2r_{p}}}\approx$$

$$\frac{1}{7\sqrt{2}}\frac{1}{\sqrt{r_{p}}}+\frac{47}{196\sqrt{2}}\sqrt{r_{p}}+\frac{2787}{10976\sqrt{2}}r_{p}^{3/2}+$$

$$\left(\frac{3165}{43904\sqrt{2}}+\frac{150\sqrt{2}}{2401}\right)r_{p}^{5/2}+\mathcal{O}(r_{p}^{7/2})\ .$$ (24)

Eq. (23) can be integrated to obtain the time for a minimum $r_{p,\rm min}$. For the limits on the integral, the eccentricity at the first quadrupole maximum is estimated from Eq. (7) and the quadrupole timescale to reach the first maximum from Eq. (9). One can stop the integration once a desired pericenter distance is reached, such as the Roche limit (see discussion of octupole level maximum eccentricity in §3.1.1). The integral is

$$\Upsilon_{\rm TP}\int_{1-\sqrt{1-\frac{5}{3}J_{z,0}^{2}}}^{r_{p,\rm min}}{d{r_{p}^{\prime}}}f(r_{p}^{\prime})\approx$$

$$-\frac{15}{64}k\frac{a_{1}^{5/2}}{a_{2}^{4}}\frac{m_{3}}{\sqrt{m_{1}}}\frac{e_{2}}{(1-e_{2}^{2})^{5/2}}\int_{t_{\text{quad}}}^{t_{r_{p,\rm min}}}dt^{\prime}\ ,$$ (25)

where $\Upsilon_{\rm TP}$ is a numerical factor that captures the evolution of $\omega_{2}$, see below. The solution of this integral provides the relevant timescale to descend to a minimum $r_{p}$ over an octupole cycle:

$$t_{\rm descent}=t_{\text{quad}}+\Upsilon_{\rm TP}t_{\rm oct}\eta(r_{p,\rm min})\ ,$$ (26)

where

$$t_{\rm oct}=\frac{64}{15}k^{-1}\frac{a_{2}^{4}}{a_{1}^{5/2}}\frac{\sqrt{m_{1}}}{m_{3}}\frac{(1-e_{2}^{2})^{5/2}}{e_{2}}\ ,$$ (27)

and

$$\eta(r_{p,\rm min})=\left[\vphantom{\frac{1}{2}}\right.\frac{\sqrt{2}}{7}\sqrt{r_{p}}+\frac{47}{294\sqrt{2}}r_{p}^{3/2}+\frac{2787}{27440\sqrt{2}}r_{p}^{5/2}+$$

$$\frac{60555}{1075648\sqrt{2}}r_{p}^{7/2}+\mathcal{O}(r_{p}^{9/2})\left.\vphantom{\frac{1}{2}}\right]\bigg{|}_{r_{p,\rm min}}^{1-\sqrt{1-\frac{5}{3}J_{z,0}^{2}}}\ .$$ (28)

Figure 3 shows a comparison of the numerical evolution with the analytical model for an example system.

To evaluate the numerical factor $\Upsilon_{\rm TP}$, we run 500 Monte Carlo realizations of systems, changing the mass ratio, semi-major axis ratio, inclination, and $e_{2}$ (see Appendix B for procedure). We find that the initial mutual inclination is the main parameter that affects the value of $\Upsilon_{\rm TP}$ (see the left panel in Figure 7, in Appendix B). Thus, given the initial mutual inclination, the descent to a minimum $r_{p}$ is well characterized. The variations in our estimation of $\Upsilon_{\rm TP}$ are $\sim$25%, depicted as the shaded area in Figure 3.

We show the percentage error between the numerical evolution and the model in Figure 4 for systematically varied initial $i_{0}$ and $\epsilon$ between $\sim$$45-80^{\circ}$ and $0.006-0.06$, respectively. Here, we fixed the masses and initial eccentricity of both orbits, and set $\omega_{1}=\omega_{2}=0$ initially. We estimated the percentage error between the model timescale depicted in Eq. (26) and the numerical descent, shown as the color code. We also over-plot the analytical flip condition from Katz et al. (2011). We note that the flip of a system is closely associated with the octupolar excitation to extreme eccentricity, and the parameter space in which deep descents occur is very similar to the parameter space in which flips occur (e.g., Lithwick & Naoz, 2011; Katz et al., 2011). In Figure 4, we marked with “X” systems that did not flip within 1000 quadrupole timescales, and we marked with circles systems that did flip, thus displaying extreme eccentricity evolution. Notably, the novel model is able to accurately describe the evolution of systems that undergo extreme eccentricity evolution to within a small factor.

Deep descents in $r_{p}$ tend to occur when $\omega_{1}=\omega_{2}=0$ or $\omega_{1}=\omega_{2}=\pi$, which is not surprising because these values are associated with the resonances (e.g., Hansen & Naoz, 2020). Therefore, in the pedagogical examples, we set our systems nominally with the initial condition $\omega_{1}=\omega_{2}=0$ to produce immediate descents. We later relax this condition and perform numerical studies that show that for most systems with sufficiently high $i$ and $\epsilon$, $\omega_{1}$ and $\omega_{2}$ will quickly synchronize to produce a descent if they are initially unsynchronized. We define the beginning of the descent in the numerical results by working backward from the maximum eccentricity until the first local eccentricity maximum having $e_{1}$ within $5\times 10^{-3}$ of the quadrupole maximum eccentricity in Eq. (7). We find that this criterion effectively captures the numerical deep descents from the local quadrupole maximum to a very large eccentricity for nearly every system.

Figure 4: Comparison of the numerical descent timescale with the analytical test particle approximation of Eq. (26) for varied $\epsilon$ and $i_{0}$. We start each system with $\omega_{1}=\omega_{2}=0$. Darker circles correspond to a lower percentage error in the timescale estimation. Circles show systems that flipped within 1000$t_{\rm quad}$ (thus displaying extreme eccentricity evolution), and we marked with “X” systems that did not flip. We also plot the flip condition from Katz et al. (2011) (red curve), which applies for systems with $i_{0}\gtrsim 60^{\circ}$. We see that the flip condition accurately separates the systems that flip from those that did not. For the systems that flipped and underwent secular descents, the model matches the numerical result to within $25\%$ for most systems, as expected given the variations observed in the Monte Carlo simulations in Appendix B. The model matches the numerical result to a factor of a few for all systems here that undergo a flip.

To highlight this behavior, we run two separate numerical tests. In one, we fix $i_{0}=75^{\circ}$ and $\epsilon=0.024$ and vary the initial $\omega_{1,0}$ and $\omega_{2,0}$ between $0-360^{\circ}$ along a grid. We first examine how long it takes for the descent to begin (as a function of the quadrupole timescale). This is shown in the panels of Figure 5, wherein the top left color codes the timescale. We demonstrate that for a wide range of the parameter space, the descent begins within 100 quadrupole cycles. We also run a Monte Carlo as a proof of concept for systems with $i_{0}=60-80^{\circ}$ and $\epsilon=0.02-0.09$, while varying the initial $\omega_{1,0}$ and $\omega_{2,0}$ between $0-360^{\circ}$ along a grid. These ranges were chosen to encompass systems within the flip parameter space of Katz et al. (2011). The timescales for the Monte Carlo run are shown in the bottom left panel of Figure 5. We find that the majority of systems start to descend within 100 quadrupole cycles.

Given a long period of time, chaos may lead nearly all such systems to eventually produce a deep descent. Once the descent begins, we find that our model is generally able to closely capture the evolution to within a small factor. This is shown in the right-hand panels of Figure 5 where we use the same two numerical runs described above and show the percentage error between the numerical descent and the analytical model, Eq. (26). As depicted, for the majority of systems, the error is less than $20\%$, for a system with fixed $\epsilon$ and fixed $i_{0}$. We further estimate the difference between the numerical and analytical models for the Monte Carlo run with varied $\epsilon$ and $i_{0}$. Only $\sim$$5\%$ of the systems have an error larger than a factor of 2, and the majority of systems have an error smaller than 20%.

Figure 5: Top left: Time until a deep descent in $1-e_{1}$ begins for varied $\omega_{1,0}$ and $\omega_{2,0}$. We show the logarithm of the time in quadrupole cycles. Here, we define a deep descent as one reaching $1-e_{1}=10^{-3}$. Every system here has initially $m_{1}=1M_{\odot},m_{2}=1M_{J},m_{3}=0.1M_{\odot},a_{1}=1$ AU, $a_{2}=20$ AU, $e_{1}=0.01$, $e_{2}=0.4$, and $i=75^{\circ}$. Darker circles correspond to a more rapid time for the descent to begin. We find that given sufficient time ($\sim$1000 quadrupole cycles), every system here eventually undergoes such a descent, with $\sim$60% happening quickly (in $<$100 quadrupole cycles). Top right: For the same set of systems as in the top left, we compare the numerical timescale of the descent once it occurs to the analytical test particle approximation in Eq. (26). Darker circles correspond to lower percentage error in the descent timescale. We find that the model agrees with the numerical evolution to within $20\%$ for $>$80% of systems and agrees within a factor of $\sim$2 for all systems. Bottom left: We vary $\omega_{1,0}$ and $\omega_{2,0}$ as in the top panels, while also varying $\epsilon$ from $0.02-0.09$ and $i_{0}$ from $60-80^{\circ}$. We plot the distribution of times until the descent begins. We again find every system eventually undergoes a descent, with $\sim$60% happening in $<$100 quadrupole cycles. Bottom right: For the same set of systems as the bottom left panel, we plot the distribution of percentage errors comparing the numerical descent timescale to the analytical test particle approximation. We find that the model agrees with the numerical evolution to within $20\%$ for $\sim$50% of systems and to within a factor of a few for all systems.

3.2 Beyond the test particle

Now we move beyond the test particle limit and consider the general case with nonzero $m_{2}$. We begin by taking the full expression for ${de_{1}}/{dt}$ given in Eq. (A6). Following the procedure in §3.1.2, we set $\omega_{1}=\pi/2$ and $\theta=\sqrt{3/5}$ to capture the eccentricity maxima. We obtain

$$\frac{d{e_{1,\text{max}}}}{dt}f({e_{1,\text{max}}}) \approx$$

$$\frac{15}{64}k\frac{a_{1}^{5/2}}{a_{2}^{4}}\frac{(m_{1}-m_{2})m_{3}}{(m_{1}+m_{2})^{3/2}}\frac{e_{2}}{(1-e_{2}^{2})^{5/2}}\cos{\omega_{2}} \ ,$$ (29)

where

$$f({e_{1,\text{max}}})=\frac{1}{(-2+9e_{1,\text{max}}^{2})\sqrt{1-e_{1,\text{max}}^{2}}}\ .$$ (30)

As in §3.1.2, we adopt $r_{p}=1-e_{1,\text{max}}$ to write

$$\frac{d{r_{p}}}{dt}f(r_{p})\approx-\frac{15}{64}k\frac{a_{1}^{5/2}}{a_{2}^{4}}\frac{(m_{1}-m_{2})m_{3}}{(m_{1}+m_{2})^{3/2}}\frac{e_{2}}{(1-e_{2}^{2})^{5/2}}\cos{\omega_{2}}\ ,$$ (31)

where $f(r_{p})$ is expanded by taking $r_{p}\ll 1$

$$f(r_{p})=\frac{1}{(9r_{p}^{2}-18r_{p}+7)\sqrt{-r_{p}^{2}+2r_{p}}}\approx$$

$$\frac{1}{7\sqrt{2}}\frac{1}{\sqrt{r_{p}}}+\frac{79}{196\sqrt{2}}\sqrt{r_{p}}+\frac{9507}{10976\sqrt{2}}r_{p}^{3/2}+$$

$$\left(\frac{10109}{43904\sqrt{2}}+\frac{1782\sqrt{2}}{2401}\right)r_{p}^{5/2}+\mathcal{O}(r_{p}^{7/2})\ .$$ (32)

Eq. (31) can be integrated

$$\Upsilon\int_{1-\sqrt{1-\frac{5}{3}J_{z,0}^{2}}}^{r_{p,\rm min}}{d{r_{p}^{\prime}}}f(r_{p}^{\prime})\approx$$

$$-\frac{15}{64}k\frac{a_{1}^{5/2}}{a_{2}^{4}}\frac{(m_{1}-m_{2})m_{3}}{(m_{1}+m_{2})^{3/2}}\frac{e_{2}}{(1-e_{2}^{2})^{5/2}}\int_{t_{\text{quad}}}^{t_{r_{p,\rm min}}}dt^{\prime}\ ,$$ (33)

where $\Upsilon$ is a numerical factor that captures the evolution of $\omega_{2}$ as discussed in §3.1.2. The solution of the integral is the relevant timescale to descend to a minimum $r_{p}$ over an octupole cycle:

$$t_{\rm descent}=t_{\text{quad}}+\Upsilon t_{\rm oct}\eta(r_{p,\rm min})\ ,$$ (34)

where

$$t_{\rm oct}=\frac{64}{15}k^{-1}\frac{a_{2}^{4}}{a_{1}^{5/2}}\frac{(m_{1}+m_{2})^{3/2}}{(m_{1}-m_{2})m_{3}}\frac{(1-e_{2}^{2})^{5/2}}{e_{2}}\ ,$$ (35)

and

$$\eta(r_{p,\rm min}) = \left[\vphantom{\frac{1}{2}}\right.\frac{\sqrt{2}}{7}\sqrt{r_{p}}+\frac{79}{294\sqrt{2}}r_{p}^{3/2}+\frac{9507}{27440\sqrt{2}}r_{p}^{5/2}+$$ (36)

$$+ \frac{526955}{1075648\sqrt{2}}r_{p}^{7/2}+\mathcal{O}(r_{p}^{9/2})\left.\vphantom{\frac{1}{2}}\right]\bigg{|}_{r_{p,\rm min}}^{1-\sqrt{1-\frac{5}{3}J_{z,0}^{2}}}\ .$$

In the test particle case (§3.1.2), the dependence of $t_{\rm oct}$ is similar; see Eq. (27).

As in the test particle case, we fit to Monte Carlo simulations to estimate values of $\Upsilon$ (see Appendix B for procedure). We again find that the initial mutual inclination is the main parameter that affects the value of $\Upsilon$ (see the right panel in Figure 7, in Appendix B). We find that the average variations in $\Upsilon$ are $\sim$25%, similar to the test particle case. Figure 3 shows a comparison of the numerical evolution with the analytical model for an example system, also compared with the test particle approximation. We find that the two approaches provide consistent results.

4 Application to Hot Jupiters

Here, we demonstrate the application of our model to a population of astrophysical systems. We estimate the EKL-driven high-eccentricity migration timescale for Hot Jupiters in stellar binaries. We use the test particle approximation in Eq. (26) to estimate the time it takes for a Jupiter-like planet in a stellar binary to descend from large distances to the Roche limit, which we define as (e.g., Guillochon et al., 2011)

$$R_{\rm Roche}=2.7R_{2}\left(\frac{m_{1}+m_{2}}{m_{2}}\right)^{1/3}\ ,$$ (37)

where $R_{2}$ is the radius of the secondary body (in this case, the planet). The Roche limit is $\sim$0.01 AU for a Jupiter-like planet and Sun-like star, requiring an orbital eccentricity of $\sim$0.99 for a planet initially at 1 AU. In Appendix C, we compare the descent timescale calculated with the model with other relevant physical timescales.

We demonstrate that the model can be used to estimate the rate of Hot Jupiter migration for a population of systems, see Figure 6. We generate two ensembles of 500 systems each with initial conditions of $m_{1}=1M_{\odot}$, $m_{2}=1M_{J}$, $m_{3}=1M_{\odot}$, $a_{1}=2$ AU, $a_{2}$ varied uniformly from 100 to 1000 AU, $e_{1}=0.01$, $i=80^{\circ}$, $\omega_{1}$ drawn uniformly from $0^{\circ}$ to $360^{\circ}$, and $\omega_{2}$ drawn uniformly from $0^{\circ}$ to $360^{\circ}$. For one set of systems, we fix $\epsilon=0.02$ and for the other we fix $\epsilon=0.05$. We choose these values because they lie in the regime sensitive to high-eccentricity excitation and produce different descent timescales. Each system is numerically evolved for 1000$t_{\rm quad}$.

We then divide each ensemble of systems into 20 bins with 25 systems each and numerically determine the rate of migration for each bin by taking $\Gamma_{\rm migration}=t_{\rm descent}^{-1}$ for the systems that undergo a descent, then normalize by the number of systems within the bin. The numerical rates are shown as the dots in Figure 6. To determine the rate of migration analytically, we take the reciprocal of Eq. (26), shown as the solid curves in Figure 6. As expected, the rates become lower as the stellar perturber becomes more distant and thus provides weaker perturbations. Importantly, the analytical model closely tracks the numerical results, with slight overestimation due to the fact that some of the numerically integrated systems begin in the chaotic regime before quickly becoming regular (as discussed in §3.1.2). Thus, these models provide an accurate method for calculating rates of migration and tidal disruption for populations of systems with varied initial conditions.

Figure 6: Comparison of the numerical (dots) and analytical (solid curves) rates of migration for a consistent number of planetary systems (see §4 for initial conditions). We see that the analytical prediction closely matches the numerical results. As $a_{2}$ increases, the rates decrease as expected due to the weaker influence of the stellar perturber.

5 Conclusion

Triple body systems are prevalent in nature, from planetary to stellar to supermassive black hole scales. The EKL mechanism seems to largely describe the dynamics of these sets of systems (e.g., Naoz, 2016). In this paper, we have presented analytical approximations for secular descents driven by the EKL mechanism in hierarchical three-body systems. These models can be used to estimate the secular descent of an object’s pericenter and are applicable to a wide variety of astrophysical systems.

EKL causes oscillations of eccentricity and inclination at both the quadrupole and octupole levels in the secular equations of motion. The octupolar envelope of the quadrupole oscillations describes the overall increase in the inner orbit’s eccentricity, i.e., the overall decrease in the inner orbit’s pericenter distance. For sufficiently high initial $i$ and $\epsilon$, the eccentricity can be driven to near unity during an octupole cycle.

We divide our discussions into two regimes, the test particle approximation, §3.1, of which $m_{2}\to 0$, and beyond the test particle, §3.2. In the former approximation, we employed two approaches. The first approach, discussed in §3.1.1, assumes small $J_{z}$ (near-perpendicular initial mutual inclination for an initially circular orbit) and provides an analytical description of the step-wise evolution of $J_{z}$. Using the evolution of $J_{z}$ in Eq. (15) and the corresponding step timescale in Eq. (18), individual steps in $J_{z}$ can be stitched together to yield the overall secular descent, see Figure 2. The second test particle approach, discussed in §3.1.2, is valid for inclinations $\gtrsim 45^{\circ}$ and is based on integrating the equation of motions (see Eq. (26)). We present a similar approach extended beyond the test particle case in §3.2, and the analogous result is Eq. (34). As expected, the timescale for descent has a similar functional form between the test and non-test particle cases, with a characteristic octupole timescale

$$t_{\rm oct}=\frac{64}{15}k^{-1}\frac{a_{2}^{4}}{a_{1}^{5/2}}\frac{(m_{1}+m_{2})^{3/2}}{(m_{1}-m_{2})m_{3}}\frac{(1-e_{2}^{2})^{5/2}}{e_{2}}\ .$$ (38)

As depicted in Figure 3, the two analytical approaches (for the test and non-test particle cases) agree well with numerical evolution for a representative system, and the analytical approach agrees with the numerical evolution to within a small factor over the entire parameter space sensitive to deep descents (see Figure 4). While these approximations are valid for nearly regular (non-chaotic) trajectories, we perform numerical tests that show most initially chaotic systems tend to quickly become regular, and the approximations in this work accurately describe the regular portion of the trajectory (see Figure 5).

In §4, we demonstrate the application of the model to Jupiter-like planets being secularly driven from large distances to the Roche limit by a stellar companion. In Figure 6, we show that the model can be used to accurately estimate rates of EKL-driven migration for a population of Hot Jupiters in stellar binaries. The ability to circumvent numerical simulations makes the models especially useful for understanding features of populations of systems. We note that the approximations can be applied to other astrophysical settings in which the EKL mechanism acts to reduce an object’s pericenter distance, such as black hole systems that produce tidal disruption events or gravitational waves.

The authors thank Isabel Angelo, Thea Faridani, and Denyz Melchor for useful discussions. The authors thank the support of NASA XRP grant 80NSSC23K0262. S.N thanks Howard and Astrid Preston for their generous support. The authors also acknowledge the use of the UCLA cluster Hoffman2 for computational resources.

Appendix A Equations of Motion

Here we present the hierarchical three-body Hamiltonian beyond the test particle. The doubly-averaged Hamiltonian (averaged over both orbits) with the nodes eliminated ($\Omega_{1}-\Omega_{2}=\pi$) up to octupole order is (see e.g., Naoz, 2016)

$$\mathcal{H}=\mathcal{H_{\text{quad}}}+\epsilon_{M}\mathcal{H_{\text{oct}}}\ ,$$ (A1)

where the quadrupole component is

$$\mathcal{H}_{\text{quad }}=C_{2}\left\{\left(2+3e_{1}^{2}\right)\left(3\cos^{2}i-1\right)+15e_{1}^{2}\sin^{2}i\cos\left(2\omega_{1}\right)\right\}\ ,$$ (A2)

the octupole component is

$$\mathcal{H}_{\text{oct}} =\frac{15}{4}C_{2}e_{1}\left\{A\cos\phi+10\cos i\sin^{2}i\left(1-e_{1}^{2}\right)\sin\omega_{1}\sin\omega_{2}\right\}\ ,$$ (A3)

and

$$\epsilon_{M}=\frac{m_{1}-m_{2}}{m_{1}+m_{2}}\frac{a_{1}}{a_{2}}\frac{e_{2}}{1-e_{2}^{2}}\ ,$$ (A4)

and

$$L_{1} =\frac{m_{1}m_{2}}{m_{1}+m_{2}}\sqrt{k^{2}\left(m_{1}+m_{2}\right)a_{1}}$$ (A5)

$$L_{2} =\frac{m_{3}\left(m_{1}+m_{2}\right)}{m_{1}+m_{2}+m_{3}}\sqrt{k^{2}\left(m_{1}+m_{2}+m_{3}\right)a_{2}}$$

$$G_{1} =L_{1}\sqrt{1-e_{1}^{2}}$$

$$G_{2} =L_{2}\sqrt{1-e_{2}^{2}}$$

$$C_{2} =\frac{k^{4}}{16}\frac{\left(m_{1}+m_{2}\right)^{7}}{\left(m_{1}+m_{2}+m_{3}\right)^{3}}\frac{m_{3}^{7}}{\left(m_{1}m_{2}\right)^{3}}\frac{L_{1}^{4}}{L_{2}^{3}G_{2}^{3}}$$

$$C_{3} =-\frac{15}{16}\frac{k^{4}}{4}\frac{\left(m_{1}+m_{2}\right)^{9}}{\left(m_{1}+m_{2}+m_{3}\right)^{4}}\frac{m_{3}^{9}\left(m_{1}-m_{2}\right)}{\left(m_{1}m_{2}\right)^{5}}\frac{L_{1}^{6}}{L_{2}^{3}G_{2}^{5}}$$

$$A =4+3e_{1}^{2}-\frac{5}{2}B\sin i^{2}$$

$$B =2+5e_{1}^{2}-7e_{1}^{2}\cos 2\omega_{1}$$

$$\cos\phi =-\cos\omega_{1}\cos\omega_{2}-\cos i\sin\omega_{1}\sin\omega_{2}\ .$$

The full equations of motion can be found in Naoz (2016). In §3.2, we use the eccentricity evolution of the inner orbit given by

$$\frac{de_{1}}{dt} =C_{2}\frac{1-e_{1}^{2}}{G_{1}}\left[30e_{1}\sin^{2}i\sin\left(2\omega_{1}\right)\right]$$ (A6)

$$+C_{3}e_{2}\frac{1-e_{1}^{2}}{G_{1}}\left[35\cos\phi\sin^{2}ie_{1}^{2}\sin\left(2\omega_{1}\right)\right.$$

$$-10\cos i\sin^{2}i\cos\omega_{1}\sin\omega_{2}\left(1-e_{1}^{2}\right)$$

$$\left.-A\left(\sin\omega_{1}\cos\omega_{2}-\cos i\cos\omega_{1}\sin\omega_{2}\right)\right]\ .$$

Moving to the test particle limit ($m_{2}\to 0$), the test particle Hamiltonian is defined in §2. Using the formalism presented there, where $J=\sqrt{1-e_{1}^{2}}$ and $J_{z}=\theta\sqrt{1-e_{1}^{2}}$, the equations of motion for a test particle are

$$\frac{dJ}{dt} = \frac{\partial F}{\partial\omega_{1}}$$ (A7)

$$\frac{dJ_{z}}{dt} = \frac{\partial F}{\partial\Omega_{1}}$$ (A8)

$$\frac{d\omega_{1}}{dt} = \frac{\partial F}{\partial e_{1}}\frac{J}{e_{1}}+\frac{\partial F}{\partial\theta}\frac{\theta}{J}$$ (A9)

$$\frac{d\Omega_{1}}{dt} = -\frac{\partial F}{\partial\theta}\frac{1}{J}\ .$$ (A10)

The conversion from $t$ to true time $T$ is given by

$$t\equiv T\frac{3}{8}\frac{m_{3}}{m_{1}}\Omega_{*}\left(\frac{a_{1}}{a_{2}}\right)^{3}\frac{1}{(1-e_{2}^{2})^{3/2}}\ ,$$ (A11)

where $\Omega_{*}$ is the orbital angular speed of the test particle.

Appendix B Fitting the model

Figure 7: Fits to $\Upsilon$ obtained from Monte Carlo simulations in the test particle case (left) and general case (right). The procedure is described in Appendix B. Fitted values of $\Upsilon$ can be used in e.g., Eq. 3.1.2 to estimate the secular descent. The average error in each case is $\sim$25%. For the general case fit displayed here, we use a fiducial value of $G_{1}/G_{2}=0.01$, but the exact value was calculated for each system to obtain the fit. In these simulations, systems below $\sim$50^∘ do not flip as they require a value of $\epsilon$ that exceeds the stability criterion to induce a flip.

We fit the models in §3.1.2 and §3.2 by comparing to Monte Carlo simulations. For the test particle case, we generate a population of 500 systems uniformly drawn with initial conditions of $m_{1}=1M_{\odot}$, $m_{3}=0.1-1.0M_{\odot}$, $a_{1}=1$ AU, $a_{2}=8-25$ AU, $e_{1}=0.01-0.1$, $e_{2}=0.1-1$, and $i=45-90^{\circ}$. These parameters are chosen to represent a range of values in which deep descents are able to occur. We repeat for the general case and draw $m_{2}$ log-uniformly from $10^{-3}-10^{-0.5}M_{\odot}$. When sampling, we enforce the stability criterion of Mardling & Aarseth (2001):

$$\frac{a_{2}}{a_{1}}>2.8\left(1+\frac{m_{3}}{m_{1}+m_{2}}\right)^{2/5}\frac{\left(1+e_{2}\right)^{2/5}}{\left(1-e_{2}\right)^{6/5}}\left(1-\frac{0.3i}{180^{\circ}}\right)\ .$$ (B1)

For the systems that undergo a flip on the initial descent (and thus display large eccentricity evolution that can be robustly fitted), we find that $\Upsilon$ depends most sensitively on the initial mutual inclination $i_{0}$. We find $\Upsilon$ to be well-characterized by the inclination modification term for the octupole timescale from Teyssandier et al. (2013):

$$\Upsilon=\frac{A}{\frac{G_{1}}{G_{2}}+\cos{i_{0}}},$$ (B2)

where $A$ is a numerical coefficient, and $G_{1}$ and $G_{2}$ are defined in Eq. (A5). For the test particle case, we obtain $A=1.83$, In the test particle case, $G_{1}=0$, so $\Upsilon$ only depends on $i_{0}$. For the general case, we calculate $G_{1}/G_{2}$ for each system and obtain a numerical factor of $A=1.87$. The fits to the simulations are shown in Figure 7. We find that the average error in both cases is $\sim$25%.

Appendix C Other physical timescales

Figure 8: Example calculations of of the time to descend to the Roche limit, $t_{\rm descent,TP}$, from Eq. (26) (solid black curve). We compare with other relevant physical timescales for a Jupiter-like planet ($m_{2}=1M_{J}$) in a stellar binary ($m_{1}=m_{3}=1M_{\odot}$). In both panels, we set $e_{1}=0.01$, $e_{2}=0.5$, and $i=70^{\circ}$. In the left panel, we fix $a_{2}=500$ AU and vary $a_{1}$. In the right panel, we fix $a_{1}=2$ AU and vary $a_{2}$. We compare with the quadrupole timescale from Eq. (9) (dashed black curve), the GR precession timescale from Eq. (C1) (orange curve), the tidal precession timescale from Eq. (C2) (blue curve), the $J_{2}$ precession timescale from Eq. (C4) (yellow curve), a typical disk dissipation timescale of 3 Myr (violet curve), and the main-sequence lifetime of a Sun-like star (green curve). GR dominates when $t_{\rm GR}<t_{\rm quad}$, and tides and $J_{2}$ become important at low $a_{1}$. The descent timescale is comparable to the disk lifetime when $a_{1}/a_{2}$ is small. For most systems, the Jupiter arrives at a close distance after the disk dissipates but during the star’s lifetime. Here, we see that the model provides an efficient method for estimating features of the Hot Jupiter population, such as the age of Hot Jupiters post-migration and the system configurations in which high-eccentricity migration is expected to occur.

In Figure 8, we compare the secular descent timescale (solid black line) with other physical timescales. Specifically, general relativity (GR) causes the inner orbit to precess in the opposite direction of the EKL precession, potentially suppressing the EKL eccentricity excitations. The GR precession timescale of the inner orbit at the first Post-Newtonian level can be estimated as (e.g., Naoz et al., 2013b)

$$t_{\rm GR}\approx 2\pi\frac{a_{1}^{5/2}c^{2}(1-e_{1}^{2})}{3k^{3}(m_{1}+m_{2})^{3/2}}\ ,$$ (C1)

where $c$ is the speed of light. It was noted before, e.g., Naoz et al. (2013b) and Hansen & Naoz (2020), that EKL eccentricity excitations are suppressed roughly when $t_{\rm GR}<t_{\rm quad}$, where $t_{\rm quad}$ is the quadrupole timescale from Eq. (9). This timescale is shown in Figure 8 as a light orange line.

Another physical process that affects the evolution is tidal precession, circulation, and shrinking of the semi-major axis. The tidal precession can suppress EKL eccentricity excitations. The tidal precession timescale (shown as a blue line in Figure 8) is estimated as (e.g., Eggleton et al., 1998; Naoz & Fabrycky, 2014)

$$t_{\rm tide}\approx\frac{a_{1}^{13/2}(1-e_{1}^{2})^{5}m_{1}m_{2}}{15k\sqrt{m_{1}+m_{2}}\left(1+\frac{3}{2}e_{1}^{2}+\frac{1}{8}e_{1}^{4}\right)\Lambda}\ ,$$ (C2)

where

$$\Lambda=m_{2}^{2}k_{L,1}R_{1}^{5}+m_{1}^{2}k_{L,2}R_{2}^{5}\ ,$$ (C3)

where $R_{1}$ is the radius of the primary body and $k_{L,1}$ ($k_{L,2}$) is the apsidal motion constant of the primary (secondary) body. Tides become important when $t_{\rm tide}$ is comparable to the other dynamical timescales, which occurs for very low $a_{1}$. In the high-eccentricity migration formation channel for Hot Jupiters, tidal forces can shrink and circularize the inner orbit.

Additionally, the oblate potential of the star can cause the planetary orbit to precess. This timescale (shown as a yellow line in Figure 8) is estimated as (e.g., Murray & Dermott, 1999)

$$t_{J_{2}}\approx\frac{4\pi}{3}\frac{1}{nJ_{2}}\left(\frac{a_{1}}{R_{1}}\right)^{2},$$ (C4)

where $n$ is the mean motion of the planet and $J_{2}$ is the gravitational quadrupole moment of the primary body. For the solar value of $J_{2}\approx 2\times 10^{-7}$ (e.g., Godier & Rozelot, 1999), the $J_{2}$ precession becomes relevant only for very low $a_{1}$ (for the consequences of evolving $J_{2}$, see e.g., Faridani et al., 2023).

We also compare to relevant lifetimes within a Hot Jupiter system. A major question concerns the role of a gas disk in forming Hot Jupiters, either as a distinct formation channel (e.g., Lin & Papaloizou, 1986) or as an aid to the high-eccentricity migration process (e.g., Vick et al., 2023). The typical timescale for a gas disk to dissipate is $\sim$3 Myr (e.g., Fedele et al., 2010), shown as a magenta line in Figure 8. We find that most Hot Jupiters formed via EKL are delivered to the Roche limit after the disk dissipates. We also compare to the $\sim$10 Gyr lifetime of a Sun-like main sequence star (green line in Figure 8) and see that Hot Jupiters migrate over a large range of timescales depending on the system configuration.