The Eccentric Kozai-Lidov mechanism for Outer Test Particle

TL;DR

This paper investigates the complex dynamics of a hierarchical three-body system with a test particle, revealing chaotic behavior and orbital flips at the octupole level, with implications for debris disk evolution and retrograde object formation.

Contribution

It extends the secular approximation to the octupole level for a test particle, analyzing chaotic dynamics and long-term evolution including general relativity effects.

Findings

Chaotic orbital flips occur at the octupole level.

Long-term evolution can produce retrograde objects.

Debris disks can become puffed up over Gyr timescales.

Abstract

The secular approximation of the hierarchical three body systems has been proven to be very useful in addressing many astrophysical systems, from planets, stars to black holes. In such a system two objects are on a tight orbit, and the tertiary is on a much wider orbit. Here we study the dynamics of a system by taking the tertiary mass to zero and solve the hierarchical three body system up to the octupole level of approximation. We find a rich dynamics that the outer orbit undergoes due to gravitational perturbations from the inner binary. The nominal result of the precession of the nodes is mostly limited for the lowest order of approximation, however, when the octupole-level of approximation is introduced the system becomes chaotic, as expected, and the tertiary oscillates below and above 90deg, similarly to the non-test particle flip behavior (e.g., Naoz 2016). We provide the Hamiltonian of the system and investigate the dynamics of the system from the quadrupole to the octupole level of approximations. We also analyze the chaotic and quasi-periodic orbital evolution by studying the surfaces of sections. Furthermore, including general relativity, we show case the long term evolution of individual debris disk particles under the influence of a far away interior eccentric planet. We show that this dynamics can naturally result in retrograde objects and a puffy disk after a long timescale evolution (few Gyr) for initially aligned configuration.

Equations

$$e_2,{\omega }_2,\theta ,{\Omega }_2,$$

$$J_2$$

$$J_{2,z}$$

$$f=f_{\mathrm{quad}}+{ϵ}_M{f}_{\mathrm{oct}},$$

$${ϵ}_M=\frac{m_1-m_2}{m_1+m_2}\frac{a_1}{a_2}\frac{e_2}{1-e_2^2}$$

$$f_{\mathrm{quad}}=\frac{(2+3e_1^2)(3{\theta }^2-1)+15e_1^2(1-{\theta }^2)\cos (2{\Omega }_2)}{(1-e_2^2{)}^{3/2}},$$

$$f_{\mathrm{oct}}$$

$$\delta =\frac{m_1-m_2}{m_1+m_2}\frac{a_1}{a_2}{e}_1$$

$${\tilde{f}}_{\mathrm{oct}}=\frac{e_2}{1-e_2^2}\frac{f_{\mathrm{oct}}}{e_1}$$

$$f=f_{\mathrm{quad}}+\delta {\tilde{f}}_{\mathrm{oct}}.$$

$$\frac{d{J}_2}{d\tau }$$

$$\frac{d{J}_{2,z}}{d\tau }$$

$$\frac{d{\omega }_2}{d\tau }$$

$$\frac{d{\Omega }_2}{d\tau }$$

$$\displaystyle\frac{d\theta}{dt}\bigg{|}_{\rm quad}$$

$$\displaystyle\frac{d\Omega_{2}}{dt}\bigg{|}_{\rm quad}$$

$$\frac{d{\Omega }_2}{d\theta }=\theta \frac{2+3e_1^2-5e_1^2\cos 2{\Omega }_2}{5e_1^2(1-{\theta }^2)\sin 2{\Omega }_2}.$$

$${\int }_{{\Omega }_b}^{{\Omega }_a}\frac{\sin 2{\Omega }_2}{2+3e_1^2-5e_1^2\cos 2{\Omega }_2}d{\Omega }_2={\int }_{{\theta }_a}^{{\theta }_b}\frac{\theta }{5e_1^2(1-{\theta }^2)}d\theta ,$$

$$\sin i_{\min }=\sin i_{\max }\sqrt{\frac{1-e_1^2}{1+4e_1^2}}.$$

$$\sin i_{\min }=\sin i_0\sqrt{\frac{(1-e_1^2)}{(1+4e_1^2)}}\mathrm{for}{\Omega }_{2,0}=0^{\circ }.$$

$$\cos 2{\Omega }_{2,\min }=\frac{2+3e_1^2-(2+8e_1^2){\sin }^2{i}_{\min }}{5e_1^2}.$$

$$t_{\mathrm{quad}}\sim \frac{4}{3}{P}_2(1-e_2^2{)}^2\frac{(m_1+m_2{)}^2}{m_1{m}_2}{\left(\frac{a_2}{a_1}\right)}^{2}circulating.$$

$$t_{\mathrm{quad}}$$

$$ϵ=\frac{a_1}{a_2}\frac{e_2}{1-e_2^2}<0.1$$

$$\delta =\frac{m_1-m_2}{m_1+m_2}\frac{a_1}{a_2}{e}_1.$$

$$\frac{d\omega_{1}}{dt}\bigg{|}_{\rm GR,inner}=\frac{3k^{3}(m_{1}+m_{2})^{3/2}}{a_{1}^{5/2}c^{2}(1-e_{1}^{2})}\ ,$$

$$\frac{d\Omega_{2}}{dt}\bigg{|}_{\rm GR,\Omega_{2}}=-\frac{3k^{3}(m_{1}+m_{2})^{3/2}}{a_{1}^{5/2}c^{2}(1-e_{1}^{2})}\ ,$$

$$t_{\mathrm{GR},{\mathrm{\Omega }}_2}\sim 2\pi \frac{a_1^{5/2}{c}^2(1-e_1^2)}{3k^3(m_1+m_2{)}^{3/2}}.$$

$$\frac{d\omega_{2}}{dt}\bigg{|}_{\rm GR,outer}=\frac{3k^{3}(m_{1}+m_{2})^{3/2}}{a_{2}^{5/2}c^{2}(1-e_{2}^{2})}\ ,$$

$$\tau =\frac{t}{16}\frac{m_1{m}_2}{(m_1+m_2{)}^2}\sqrt{\frac{G_N(m_1+m_2)}{a_2^3}}{\left(\frac{a_1}{a_2}\right)}^2=\frac{t}{16}\frac{m_1{m}_2}{(m_1+m_2{)}^2}\frac{2\pi }{P_2}{\left(\frac{a_1}{a_2}\right)}^2,$$

Full text

1. INTRODUCTION

The hierarchical three body secular dynamics has been studied extensively in the literature and was shown to be very effective in addressing different astrophysical phenomena (see for review Naoz, 2016, and reference therein). In this hierarchical setting the inner binary is orbited by a third body on a much wider orbit, the outer binary, such that the secular approximation can be applied (i.e., phase averaged, long-term interaction). The gravitational potential is then expanded in semi-major axis ratio ($a_{1}/a_{2}$, which, in this approximation, remains constant), where $a_{1}$ ($a_{2}$) is the semi-major axis of the inner (outer) body (Kozai, 1962; Lidov, 1962). This ratio is a small parameter due to the hierarchical configuration. The lowest order of approximation, which is proportional to $(a_{1}/a_{2})^{2}$ is called the quadrupole-level.

Most of these studies focus on the gravitational perturbations that a far away perturber exerts on the inner binary. In early studies of high-inclination secular perturbations (Kozai, 1962; Lidov, 1962), the outer orbit was assumed to be circular and it was assumed that one of the inner binary members is a massless test particle. In this situation, the component of the inner orbit’s angular momentum along the z-axis (which is set to be parallel to the total angular momentum, i.e., the invariable plane) is conserved, and the lowest order of the approximation, the quadrupole approximation, is valid. However, relaxing either one of these assumptions leads to qualitative different behavior (e.g., Naoz et al., 2011; Lithwick & Naoz, 2011; Katz et al., 2011). Considering systems beyond the test particle approximation, or a circular orbit, requires the next level of approximation, called the octupole–level of approximation, which is proportional to $(a_{1}/a_{2})^{3}$ (e.g. Harrington, 1968, 1969; Ford et al., 2000; Blaes et al., 2002). In the octupole level of approximation, the inner orbit eccentricity may reach extreme values (Ford et al., 2000; Naoz et al., 2013a; Li et al., 2014c; Teyssandier et al., 2013). In addition, the inner orbit can flip its orientation, with respect to the total angular momentum (i.e., z-axis), from prograde to retrograde (Naoz et al., 2011).

Here we study the secular evolution of a far away test particle orbiting an inner massive binary. In this case, the inner orbit is fixed, and effectively carries all of the annular momentum of the system, while the outer orbit undergoes a dynamical evolution. This situation has large range of applications from the gravitational perturbations of binary super massive black holes on the surrounding stellar distribution to the effects of planetary orbits on debris disks, Oort cloud and trans-neptunian objects. From N-body simulations, Zanardi et al. (2017) analyzed the long term evolution of test particles in the presence of an interior eccentric planet. Such an study produces particles on prograde and retrograde orbits, as well as particles whose orbital plane flips from prograde to retrograde and back again along their evolution.

We note that Ziglin (1975) investigated the oscillations of an outer circumbinary planet in the context of the restricted elliptical three body problem. Later Verrier & Evans (2009) and Farago & Laskar (2010) studied the stability of high inclined planet around in this situation using a combination of numerical and perturbation theory up to the quadruple level of approximation approaches (see also, Li et al., 2014a; de la Fuente Marcos et al., 2015). Furthermore, Gallardo (2006) and Gallardo et al. (2012) studied the effects of the Kozai-Lidov for trans-neptunian objects near mean motion resonance with Neptune. However, here, we do not allow for mean motion resonances to allow for the double averaging process (see Naoz et al., 2013a, appendix A2 for the canonical transformation which describes the averaging process).We provide a general treatment for the outer-test particle case, up to the octupole–level of approximation in the secular theory.

The paper is organized as follow: We begin by describing the outer test particle Hamiltonian and equations of motion (§2), and continue to discuss the quadrupole-level of approximation (§3) where we also drive the relevant timescales, and then we study the role of the octupole-level of approximations, and provide surface of section maps (§4). We also discuss the role of general relativity precession in §5. We then consider one study case in the form of the long term evolution debris disk particles in §6. Finally, we offer our discussions in §7.

2. The Equations of Motion

We solve the orbit of an exterior massless test particle to an eccentric planet ($m_{2}$), both orbiting a star ($m_{1}$), including only secular interactions expanded to octupole order. The planet is on a fixed eccentric orbit (i.e., $e_{1}={\rm const}$) and the outer particle’s orbit is specified by four variables:

[TABLE]

where $e_{2}$ is the test particle eccentricity, $\theta=\cos i$ and $i$ is the inclination of the test particle with respect to the inner orbit, and $\omega_{2}$ and $\Omega_{2}$ are the argument of periapse and longitude of ascending node of the outer orbit, relative to the inner planet’s periapse (Murray & Dermott, 2000). Specifically, We set $\varpi_{1}=0$, see Appendix A for the coordinate transformation. We kept the subscript “2” in $\omega_{2}$ and $\Omega_{2}$ for consistency with the comparable masses treatments. From $e_{2}$ and $\theta$ we can define the canonical specific momenta

[TABLE]

The Hamiltonian for which $m_{3}\to 0$ is

[TABLE]

where

[TABLE]

and

[TABLE]

[TABLE]

Note that unlike the inner test particle approximation $\epsilon_{M}$ is not constant during the motion, and the constant parameter during the evolution is:

[TABLE]

and the Hamiltonian up to the octupole level of approximation can be defined as:

[TABLE]

and

[TABLE]

We note that $f_{\rm quad}$ has the same functional form as the inner test particle $F_{\rm quad}$ presented in Lithwick & Naoz (2011) up to the $(1-e_{2}^{2})^{3/2}$ which is not constant in our case.

The equations of motion may be expressed as partial derivatives of an energy function $f(e_{2},\omega_{2},\theta,\Omega_{2})$

[TABLE]

where $\tau$ is proportional to the true time (see Eq. [A1]). Unlike inner orbit test particle $\epsilon_{M}$ is not constant while $e_{1}$ is constant (in other words, the angular momentum of the inner orbit is conserved). The equations of motion were tested successfully compared to the general equations of motions, presented in Naoz et al. (2013a). We also test the evolution compared to N-body in Appendix B.

3. quadrupole-level of approximation

3.1. General Analysis

The quadrupole-level of approximation is integrable and thus provides a good starting point. Unlike the quadrupole-level approximation for the inner orbit test particle, here the z-component of the particle’s angular momentum is not conserved, as the Hamiltonian depends on $\Omega_{2}$. However, at this level $J_{2}$ is conserved, and thus the outer orbit eccentricity $e_{2}$ remains constant as the Hamiltonian does not depends on $\omega_{2}$. The equation of motion for the inclination takes a simple form:

[TABLE]

where we consider the time evolution and not the scaled evolution for completeness, and $P_{2}$ is the period of the outer orbit. We can find the maximum and minimum inclination by setting $\dot{\theta}=0$. Thus, we find that the values of the longitude of ascending nodes that satisfy this condition are $\Omega_{2}=n\pi/2$, where $n=0,1,2..$. In other words, $\Omega_{2}$ has two classes of trajectories, librating and circulating. The trajectories in the librating region are bound between two values of $\Omega_{2}$, while the circulation region represents trajectories where the angles are not constrained between two specific values. On circulating trajectories, at $\Omega_{2}=0$, the inclination ($i<90^{\circ}$) is largest (where the $i>90^{\circ}$ case is a mirror image of the $i<90^{\circ}$ one). The extrema points for the librating mode are located at $\Omega_{2}=90^{\circ}$. The time evolution of $\Omega_{2}$ for the quadrupole-level of approximation is given by:

[TABLE]

In Figure 1 we show the evolution associated for the quadrupole-level of approximation. The two librating and circulating trajectories are considered, where we folded the $\Omega_{2}$ angle to be between $0-180^{\circ}$. The librating mode gives the nominal precession of the nodes, at which the inclination oscillates between the $i_{90}$ inclination (the inclination for which $\Omega_{2}=90^{\circ}$) and $180^{\circ}-i_{90}$. The precession of the nodes was noted before in the literature (e.g., Innanen et al., 1997).

From the latter Equation and Equation (3.1) we have

[TABLE]

Integrating the two sides we have:

[TABLE]

where $\Omega_{a,b}$ is the longitude of ascending nodes that is associated with the inclination value of $\theta_{a,b}=\cos i_{a,b}$.

For the circulating mode, we find that, after integrating over Equation (17) from $i_{\rm max}$ to $i_{\rm min}$ (and $\Omega_{2}=0$ to $\Omega_{2}=90^{\circ}$, respectively) we get

[TABLE]

Note that this expression can be also achieved by considering the conservation of energy between the minimum and maximum cases. Setting the initial conditions for the energy, Equation (6), we can find the extrema points as a function of the initial conditions. A special case can be considered when $\Omega_{2}$ is initially set to be zero, and then the maximum inclination is the initial inclination $i_{0}$, in other words:

[TABLE]

In other words we can set $i_{\rm max}\to 90^{\circ}$ and for a given $e_{1}$ find the largest $i_{\rm min}$ allowed, which corresponds to the separatrix. This relationship is also apparent in figure 14 in Zanardi et al. (2017) numerical results. As can be seen, for the circulating mode depicted in Figure 1, setting initially $i=20^{\circ}$ corresponds to $i_{\rm min}=16.25^{\circ}$ and $10.56^{\circ}$ for $e_{1}=0.3$ and $e_{1}=0.9$, respectively, consistent with Equation (20). A comparison between the flip criterion and the numerical results can be seen in Figure 9.

For the librating mode we set $i_{\rm min}$ which associated with the $\Omega_{2}=90^{\circ}$ case, and thus integrating over Equation (17) between $i=90^{\circ}$ to $i_{\rm min}$ (and $\Omega_{2}=\Omega_{2,{\rm min}}$ to $\Omega_{2}=90^{\circ}$, respectively), we get,

[TABLE]

Thus, setting initially $\Omega_{2}=90^{\circ}$ as in the examples depicted in Figure 1 we find that the minimum value that $\Omega_{2}$ can achieve in the $e_{1}=0.9$ case is $59.23^{\circ}$ for the $i_{\rm min}=60^{\circ}$ and $15.95^{\circ}$ for the $i_{\rm min}=20^{\circ}$ example, consistent with the numerical results. Hence, the range of which $\Omega_{2}$ is librating on is $2\times(90^{\circ}-\Omega_{2,{\rm min}})$

3.2. Timescale

The timescale associated with the evolution can be estimated from the equation of motion for $\Omega_{2}$ at the quadrupole-level [Equation (3.1)], for the circulating mode, by setting $d\Omega_{2}\to\pi$ and taking the terms in the parenthesis to be roughly order of unity (which is achieved, by setting $\Omega_{2}\to 0$):

[TABLE]

For the example system depicted in Figure 1, this equation gives a timescale of about $6\times 10^{7}$ yr, for initial inclination of $60^{\circ}$, which agrees with the circulating mode (although we note that different $e_{1}$ give slightly different timescales). We also estimate the timescale in the librating mode by setting $d\Omega_{2}\to 2\times(90^{\circ}-\Omega_{2,{\rm min}})$, which we have found earlier. Thus,

[TABLE]

The $2$ pre-factor of here comes from numerical comparisons to the examples depicted in Figure 1. Note the $e_{1}$ dependency that rises from $\Omega_{2,{\rm min}}$. This timescale is consistent with the examples depicted in Figure 1 by less than a factor of two.

4. The role of the octupole level of approximation

The octupole level of approximation can significantly affect the overall dynamics of the general hierarchical three body system (see Naoz, 2016, and reference therein). Specifically, in the inner test particle case, the inner orbit’s z component of the angular momentum is not conserved anymore and the orbit is allowed to flip (for large range of initial inclinations Lithwick & Naoz, 2011; Li et al., 2014c, b). In our case, the $J_{2,z}$ is not conserved at the quadrupole-level, but $J_{2}$ is. Thus, the octupole level of approximation in this case allows for variations of $e_{2}$ and introduces higher level resonances, which may result in a chaotic behavior (see below).

In Figure 2 we consider two representative example for which we compare the quadrupole (blue lines) and octupole (red lines) levels of approximation, where we consider the time evolution of the eccentricity and inclination. In both of these examples we consider a $1$ M⊙ star orbited by an eccentric Jupiter, with a test particle on a far away orbit. One can consider such a setting to represent a result of a scattering event for example.

On the left set of panels of Figure 2 we consider a Jupiter at $0.4$ AU with $e_{1}=0.65$ and a test particle at $7$ AU with $e_{2}=0.4$, initialized on a retrograde orbit ($i=91^{\circ}$). With the introduction of the octupole level of approximation to the calculation, the test particle eccentricity starts to oscillate, though in this case it never increases pass its initial value (due to choice of initial conditions here). More notably, the test particle inclination, with respect to the total angular momentum, oscillates from retrograde ($>90^{\circ}$ which was the initial condition) to prograde ($<90^{\circ}$). As in the more general case, there is no apparent associated timescale for this flipping modulation and it seems chaotic in nature (see below). While the quadrupole-level is circulatory in nature (see $i-\Omega_{2}$ plot) a libration behavior emerges at the octupole level.

On the right set of panels of Figure 2 we consider a Jupiter at $3$ AU with $e_{1}=0.85$ and a test particle at $50$ AU with $e_{2}=0.7$, and the system is initialized on a prograde orbit ($i=20^{\circ}$). Here, like the previous example the outer test particle eccentricity begins to oscillate, and even grows above the initial value. However, unlike the previous example this system does not flip. The inclination does oscillate with a long scale modulation and we show the long scale evolution that captures about four octupole cycles. The system does not exhibit a chaotic behavior in this case, and it remains in a circulatory trajectory even after the inclusion of the octupole level of approximation to the calculation.

In Figure 3 we zoom-in on the evolution of a different example, and also provide the time evolution of $\omega_{2}$ for the octupole level itself. As in the general hierarchical secular three body problem we find the short scale (associated with the quadrupole) oscillations, that are modulated by the higher level octupole approximation. The octupole modulations take place on timescales which is between few$\times t_{\rm quad}$ to few tens$\times t_{\rm quad}$.

Its interesting to note that the inclination flips shown in Figures 2 and 3 are qualitatively different from the ordered, back and forth oscillation of the quadrupole-level of approximate evolution. The latter produces a simple, ordered oscillation of the inclination angle between $i_{0}$ and $180^{\circ}-i_{0}$, for the librating regime. However, in the presence of the octupole-level of approximation the system behaves similarly to the general flips discussed in Naoz et al. (2011), where the inclination oscillates for sometime at the prograde ($i<90^{\circ}$) regime, and then flips to retrograde configuration ($i>90^{\circ}$).

The eccentricity, $e_{2}$, gives rise to an additional complication as in essence the outer test particle eccentricity can grow so much until the orbits will cross. We adopt the nominal stability criterion

[TABLE]

to guid us when the system leaves stability. We discuss this stability criterion in the context of N-body comparisons in appendix B.

To explore the chaotic nature of the system and the different dynamical regimes we use surface of sections. The outer test particle approximation reduces the general hierarchical three body system from six degrees of freedom to four degrees of freedom. In addition, in the test particle limit, the inner orbit is stationary and reduces the system to two degrees of freedom. In this system $f$ and $\delta$ are the only conserved parameters and $\omega_{2}$ and $\Omega_{2}$ are the only coordinates that can change with time. For a two-degrees-of-freedom system, the surface of section projects a four-dimensional trajectory on a two-dimensional surface, where we select intersections of the trajectories on the surfaces when $\omega_{2}$ and $\Omega_{2}$ moves in the positive directions. For simplicity, we separate the two initial conditions into three characteristic parameters: $e_{1}$, the inner orbit eccentricity, which remains constant during the evolution, the energy, or initial value of the reduced Hamiltonian $f$, and

[TABLE]

Note that the energy depends on $\delta$ and $e_{1}$, so in fact, although we choose to characterize the surface of sections by three parameters, there are only two independent ones.

In Figures 4-6 we consider the surface of sections for various values of the $f$, $e_{1}$ and $\delta$, in the $J_{2}-\omega_{2}$ plane (top row in each figure) and in $J_{2,z}-\Omega_{2}$ plane (bottom row in each figure). In both planes we identify the resonances at which the momenta and angles undergo bound oscillations. The trajectories in this region are quasi-periodic, and the system is in the libration mode. The circulation region represents trajectories where the angles are not constrained between two specific values. Both librating and circulatory trajectories are mapped onto a one-dimensional manifold on the surface of section and they form lines on the section. However, chaotic trajectories are mapped onto a higher dimensional manifold and they are filling an area on the surface. We note that in some of the trajectories in the Figures, due to sampling limitation, seem as dashed lines, but they actually represent a one-dimensional manifold. In all of the maps we indicate the instability regime (light orange stripe) for which $\epsilon>0.1$.

We intersect the trajectories at $\Omega_{2}=\pi/2$ to produce the surfaces in the $J_{2}-\omega_{2}$ plane, in order to capture the librating cases. The empty regions at large $J_{2}$ at the parameter spaces in the far left and right panels in Figure 5 (i.e,. for the parameters: $e=0.4$, $\delta=0.02$ and $f=\pm 10$) and the far left panel in Figure 6 ( for the parameters: $e=0.9$, $\delta=0.1$ and $f=-10$) correspond to regions with no physical solutions. The variabilities in $J_{2}$ are mostly small in the stable regime. We see that there are regular behavior, i.e., trajectories which fill one-dimensional lines on the surface of section in most of the stable region. We find the emergence of chaos in parts of the unstable zones, in particular when $J_{2}$ is low ($e_{2}$ is high).

Considering the $J_{2,z}-\Omega_{2}$ plane, we intersect the trajectories at $\omega_{2}=0$. The system exhibits a chaotic behavior across all parameter regime of $e_{1},\delta$ and $f$. Most of the circulation region, associated with curve, non chaotic one-dimensional manifold, are typically associated with |J_{2,z}|\mathrel{\hbox to0.0pt{\lower 4.0pt\hbox{\hskip 1.0pt\sim}\hss}\raise 1.0pt\hbox{>}}0.3. The outer orbits can flip ($J_{2,z}$ shifts signs) in most of the parameter space.

Resonances can be easily identified in a few of the maps. Specifically, in the $J_{2}-\omega_{2}$ plots, the resonances can be found centered near $\bf\omega_{2}=\pi$ (e.g., $e_{1}=0.4$, $\delta=0.02$, $f=-10$, and $e_{1}=0.9$, $\delta=0.1$, $f=-10$, etc.), and $\omega_{2}=\pi/2$ (e.g., $e_{1}=0.4$, $\delta=0.02$, $f=10$). The dynamics is quite complicated when $e_{1}$ is higher and when $\delta$ is larger, and higher order resonances (appearing as small liberating islands) emerge in the surface of section in the $J_{2}-\omega_{2}$ plane when $e_{1}=0.9$, $\delta=0.1$ and $f=0$. Resonances can also be identified in the $J_{2,z}-\Omega_{2}$ plane, such as $e_{1}=0.9$, $\delta=0.1$, $f=10$.

5. The Role of the General Relativity

As noted previously in many studies, General Relativity (GR) precession tends to suppress the inner orbit eccentricity excitations associated with the Eccentric Kozai-Lidov (EKL) mechanism, and thus suppress the flips (e.g., Naoz et al., 2013b). In our secular case, the inner orbit is massive and the outer orbit is a test particle, so practically the inner orbit does not feel the outer orbit gravitational interactions. However, the inner orbit can still precess due to GR with the nominal precession rate (e.g., Naoz et al., 2013b):

[TABLE]

where $k^{2}$ is the gravitational constant and $c$ is the speed of light. However, in our frame of reference, where the inner orbit carries all the angular momentum, we are basically working in the rotating frame of the inner orbit. Therefore, since we set $\omega_{1}=-\pi-\Omega_{2}$, (see Appendix A), GR precession of $\omega_{1}$ translates to a precession of $\Omega_{2}$. Thus, using our coordinate transformation we find:

[TABLE]

which can suppress the inclination oscillations. The timescale associated with that precession is the nominal GR one:

[TABLE]

In Figure 7 we consider two examples, one for which $t_{\rm GR,\Omega_{2}}\sim 5\times 10^{6}$ yr is smaller than $t_{\rm quad}\sim 6\times 10^{6}$ yr (top panel), and the other of which $t_{\rm GR,\Omega_{2}}\sim 4.4\times 10^{8}$ yr is a bit longer than the corresponding $t_{\rm quad}\sim 2\times 10^{7}$ yr (bottom panel). Both of these examples had a Sun size star and a Jupiter planet, orbited by a far away test particle. The Jupiter has a non-negligible eccentricity, that perhaps can be a result of either a scattering event or a high eccentricity migration. In the top panel the Jupiter was set at $a_{1}=0.5$ AU with $e_{1}=0.4$ and the test particle was set at $a_{2}=10$ AU and $e_{2}=0.6$. The system in the absence of GR was in libration mode and as noted before exhibited a chaotic nature. However, the GR precession bound the system into a circulatory regime, and suppressed the flips. In the bottom panel the Jupiter was set at $a_{1}=3$ AU and $e_{1}=0.9$, while the test particle was set at $a_{2}=40$ AU with $e_{2}=0.65$. Note that in this latter case, although the Jupiter is rather far from the host star and the GR precession timescale is longer than the secular precession timescale, the GR changes the dynamics. Specifically, before the inclusion of GR precession, the system was in a libration mode and seemed quasi-periodic, however, after the inclusion of GR the system exhibits both libration and circulation, and the emergence of chaotic behavior seems to take place. The dramatic change in dynamical behavior with the inclusion of GR precession, even if takes place on longer timescales than the secular timescales, was noted previously in the general and inner test particle case in Naoz et al. (2013b).

We note that in all of our calculation below we also take into account the outer orbit GR precession (e.g., Naoz et al., 2013b):

[TABLE]

which typically takes place on much larger timescales.

6. A study case application: individual Debris disk particles

Debris disks mark the late end stages of planet formation, and are made of the leftover material of rocks and ices. The gravitational interactions between these particles and interior or exterior companion can leave a distinct imprint on the morphology of the disk and can cause dust production (e.g., Rodigas et al., 2014; Matthews et al., 2014; Nesvold & Kuchner, 2015; Lee & Chiang, 2016; Nesvold et al., 2016). Many of these studies typically focus on few million years of integration to allow for comparison of observations which usually can detect young systems. Here we allow for longer integration timescales and investigate the evolution of a test particle under the influence of an eccentric Jupiter.

In Figure 8 we show an example system, where we consider an eccentric Jupiter at $5$ AU, with $0.85$ eccentricity orbiting a one solar mass star. The test particle is located at $55$ AU and initialized with an eccentricity of $0.5$. We integrate the octupole level equations of motion in the presence of GR precession for both the inner and outer orbits. As in the example depicted in Figure 7, which also considered an icy body or comet reservoirs analogs, the orbit switches between libration and circulation as can be seen in the bottom panel.

We consider the effect of the planet’s eccentricity, $e_{1}$, and the test particle inclination by surveying the parameter space of $e_{1}$ and initial inclination for a given system, where we set $a_{2}=55$ AU, with $e_{2}=0.5$, and set the system initially with $\Omega_{2}=90^{\circ}$ and $\omega_{2}=0^{\circ}$. Many giant exoplanets have high eccentricities, specifically for giant planets ($m\sin i>0.1$ MJ, with separation $>0.05$ AU) the average eccentricity is $\sim 0.2$, and a maximum value of $0.97$111Taken from The Exoplanet Orbit Database (Wright et al., 2011).. Thus, an eccentric Jupiter doesn’t seem like an unlikely configuration for a planetary system. In the example depicted in Figure 8 the time evolution of the test particle’s inclination that starts with a moderate eccentricity oscillates between extreme values $20^{\circ}\sim 160^{\circ}$ without increasing its eccentricity.

The system depicted in Figure 8 is, of course, just one example for a particular choice of the orbital parameters. To study the effects of planet’s eccentricity and initial inclination we have systematically explored the $e_{1}$-$i$ parameter space in Figure 9. We choose a Jupiter like system ($a_{1}=5$ AU) and for a range of eccentricities, with an outer orbit at $55$ AU. The test particle orbit was initialized with $e_{2}=0.5$, $\Omega_{2}=90^{\circ}$, $\omega_{2}=0^{\circ}$ and a range of inclinations. We show the maximum inclination reached during the evolution as a function of the initial inclination and eccentricity in Figure 9. In the left panel we depict the initial inclination vs the initial Jupiter’s eccentricity, where the color code marks the maximum inclination reached. The solid line in the left panel follows Equation (19), which is consistent that the resonance associated with the quadrupole-level of approximation is indeed the main driver for the dynamical evolution of the system. In the figure we depicted the initial inclination regime to the prograde case ($i_{\rm initial}\leq 90^{\circ}$) to avoid clutter. However, we show in the right panel the maximum inclination reached during the evolution as a function of the initial inclination going all the way to $180^{\circ}$ this time.

As a proof of concept we depict in Figure 10 the behavior of a narrow debris disk after $4$ Gyr of evolution. This inclination represents the instantaneous inclination at $4$ Gyr of evolution. The system, of course, continue to oscillate, and the disk of particles will remain puffed. The inclination and eccentricity of the system at this snapshot are qualitatively different from the initial conditions assumed. This hypothetical system shows the orbital configuration of a disk located between $55-65$ AU, with an interior eccentric Jupiter ($e_{1}=0.85$) at $5$ AU around a solar mass star. The system was set initially with a mutual inclination of $20^{\circ}$, $e_{2}=0.3$ and $\Omega_{2}$ and $\omega_{2}$ are chosen from a random uniform distribution between $0-360^{\circ}$. While the particles in the disk have eccentric values, the disk does not appear as a coherent eccentric ring as the values of $\Omega_{2}$ and $\omega_{2}$ are random. At the end of the integration, the particles in the disk became slightly more eccentric (with an average eccentricity of $\sim 0.34$) and there is a clear trend of $\Omega_{2}$ as a function of inclination. The particles with inclination above the initial $20^{\circ}$ have a value of $\Omega_{2}$ close to zero (or with the symmetric value $180^{\circ}$), while the particles with inclination below $20^{\circ}$ have $\Omega_{2}$ values closer to $90^{\circ}$. The behavior is singular to $\Omega_{2}$ and is not manifest itself in $\omega_{2}$ (does not depicted here to avoid clutter). We note that some retrograde particles were formed as well, in line with Zanardi et al. (2017) numerical results of an orbital flip. These behaviors are easily understandable from the surface of section maps depicted above.

Its important to note that during the evolution of the system depicted in Figure 10, the maximum $\epsilon$ achieved was $0.0625$. The average value of $\epsilon$ which corresponds to the maximum $e_{2}$ achieved during the evolution was $\sim 0.4$. Thus, the system is kept stable during this evolution, the secular approximation holds, and we do not expect any scattering event. We also note that we have compared a debris disk particles secular and N-body evolution and found a qualitative agreement, which is similar in behavior to Figure 11 left panel, and this is not shown here to avoid clutter.

7. Discussion

We have studied the secular evolution of an outer test particle hierarchical system. We presented the three body, outer test body Hamiltonian up to the octupole level of approximation in the power series of the semi-major axis ratio. We showed that in the quadrupole-level of approximation, $(a_{1}/a_{2})^{2}$, the system has two distinct behaviors, librating and circulating (see Figure 1), where the librating mode gives the nominal precession of the nodes, results for which the inclination oscillates between the $i_{90}$ inclination (the inclination for which $\Omega_{2}=90^{\circ}$ and $180^{\circ}-i_{90}$). Furthermore, the bound values of the liberating mode have a simple analytical expression, Eq. (21). We also found the minimum and maximum inclination that the system can reach in the circulating mode (see Equations (19) and (20)). These conditions are sensitive to the initial inner orbit’s eccentricity, and are nicely reproduced in numerical testing here (see Figure 9) and in Zanardi et al. (2017) numerical experiments (see their figure 14). We also estimated the timescale for oscillations for the two modes (see Section 3.2).

We then showed that introducing the octupole level of approximation allows for transition between the two libration and circulation modes (see Figures 2 and 3). This yields that the overall dynamics of the system is similar to the behavior of the general flip behavior in the eccentric Kozai-Lidov mechanism. In particular, the dynamics is quite chaotic for parameter regions with high $e_{2}$ and perpendicular mutual inclinations (when $J_{2}$ is low and when $J_{2,z}$ is near zero), as shown in the surface of sections Figures 4-6.

General relativity can play an important role in suppressing or exciting the eccentricities in the hierarchical three body problem (e.g., Naoz et al., 2013b). We find here similar behavior. Specifically, the inclination excitation will be suppressed for systems with GR precession faster than the quadrupole precession. However, when GR precession takes place on similar (or even somewhat larger) timescales to that of the quadrupole precession, the additional precession can produce inclination excitations, in a non-regular manner (see Figure 7).

The dynamics of these type of systems can have a wide range of applications, from stars around supermassive black hole binaries to the evolution of individual debris disk particles. We have chosen the latter as an example and presented a typical example of the evolution of a test particle due to the gravitational perturbations from an eccentric Jupiter (see Figure 8). We systematically varied the Jupiter’s eccentricity and the outer orbit’s inclination, where we found an agreement between the analytical relation for crossing the $90^{\circ}$ threshold, and the numerical tests. This also suggests that eccentric planet can pump up the inclination of icy body or comet reservoir analogs (Figure 9). This was further supported by considering the evolution of an initially narrow, thin disk of test particles exterior to an eccentric planet, with initial mutual inclination of $20^{\circ}$. The disk became puffed with some particles on a retrograde orbits (see Figure 10). A detailed study of the effects of eccentric planets on exterior test particle is presented in Zanardi et al. (2017). This mechanism will also be important for example to circumbinary planetary systems (considered first by Ziglin, 1975) and the stars near mergers of black hole binaries, but detailed studies is beyond the scope of this paper.

We thank the referee for a quick and detailed report, and especially for his/her inquiry about the surface of sections. We also thank Vladislav Sidorenko for pointing out some missing references. SN acknowledges partial support from a Sloan Foundation Fellowship. GL is supported in part by the Harvard William F. Milton Award. MZ, GdE, and RPD acknowledge the financial support given by IALP, CONICET, and Agencia de Promoción Científica, through the PIP 0436/13 and PICT 2014-1292.

Appendix A A. Orbital Parameters and the Scaled Time

One might have expected that $\omega_{1}={\rm const}$ when the outer particle is massless. But in truth $\Omega_{1}$ is undefined because the reference plane is aligned with the inner orbit. Therefore, the inner planet must only have $\omega_{1}+\Omega_{1}={\rm const}$, and we may choose without loss of generality the constant to equal zero. Hence, elimination of the nodes (i.e., $\Omega_{1}-\Omega_{2}=\pi$) implies $\omega_{1}=-\pi-\Omega_{2}$.

Similarly to the treatment done in Lithwick & Naoz (2011) we have rescaled the momenta by an arbitrary constant to achieve the specific angular momentum. The Hamiltonian is rescaled by the same constant, and we find that the rescaled time is:

[TABLE]

where $t$ is the true time. The numerical factor $16$ comes by taking $m_{3}\to 0$ in the general Hamiltonian (see Naoz, 2016, for the general form of the hierarchical three body double averaged hamiltonian). There is a choice to be made, to either scale the Hamiltonian by this numerical factor or $\tau$. Here, we choose to absorb this number in $\tau$ to be consistent with the inner test particle Hamiltonian.

Appendix B B. Comparison with N-Body

In this section, we compare the secular approximation at the octupole level in the test particle limit with the N-body simulation, using Mercury code (Chambers & Migliorini, 1997). In this comparison we did not include GR. Good agreements can be reached when the apocenter distance of the inner binary is much smaller than the pericenter distance of the outer binary. In particular, we include an illustrative example here in Figure 11, where we consider an eccentric ($e=0.9$ Jupiter at $3$ AU around a solar like object. The test particle is set at $40$ AU with $e_{2}=0.65$. The system is set initially with $i=20^{\circ}$, $\omega_{2}=90^{\circ}$ and $\Omega_{2}=100^{\circ}$.

As shown in Figure 11, the secular approximation (shown as red lines) agrees qualitatively well with the N-body results for the eccentricity and inclination oscillations. This is likely due to the double averaging process, but nonetheless, the maximum and minimum of the orbital parameters are conserved in both N-body and secular calculations. It is interesting to note that similarly to (Lithwick & Naoz, 2011) the approximation holds as long as $\epsilon<0.1$. However, unlike the inner-test particle case, $e_{2}$ can change and increase during the evolution which may break the validity of the approximation during the evolution, this is shown in Figures 12.