{'fraction_non_alphanumeric': 0.037477558348294436, 'fraction_numerical': 0.019765916309901946, 'mean_word_length': 4.3420324603467355, 'pattern_counts': {'":': 0, '<': 21, '<?xml version=': 0, '>': 16, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Mildly-Hierarchical triple dynamics and applications to the outer solar system

Three-body interactions are ubiquitous in astrophysics. For instance, Kozai-Lidov oscillations in hierarchical triple systems have been studied extensively and applied to a wide range of astrophysical systems. However, mildly-hierarchical triples also play an important role, but they are less explored. In this work we consider the secular dynamics of a test particle in a mildly-hierarchical configuration. We find the limit within which the secular approximation is reliable, present resonances and chaotic regions using surface of sections, and characterize regions of phase space that allow large eccentricity and inclination variations. Finally, we apply the secular results to the outer solar system. We focus on the distribution of extreme trans-neptunian objects (eTNOs) under the perturbation of a possible outer planet (Planet-9), and find that in addition to a low inclination Planet-9, a polar or a counter-orbiting one could also produce pericenter clustering of eTNOs, while the polar one leads to a wider spread of eTNO inclinations.


INTRODUCTION
The three-body problem is one of the oldest problems in astrophysics. No general solution is possible due to its chaotic nature; however, important results can be obtained using perturbative approaches. For instance, when one of the objects is much farther from the other two objects (in a hierarchical configuration), one can treat the influence of the further object as a weak perturbation and discover interesting dynamical properties. Kozai (1962) and Lidov (1962) developed the framework to study such three-body systems. It was found that when the inner binary (composed of the two closely separated objects) are inclined (i m > 40 • ) with respect to the orbit of the farther object (outer orbit), Kozai-Lidov resonances could lead to large amplitude eccentricity and inclination variations of the inner binary.
Due to the prevalence of three-body systems in hierarchical configurations, the Kozai-Lidov mechanism has since been applied to explain a wide variety of astrophysical phenomena. Inspired by the objects in our solar system, the Kozai-Lidov mechanism focused on systems with a circular outer orbit. In this case, expanding the disturbing function to second order in the small parameter (semi-major axis ratio of the inner and outer orbit, a 1 /a 2 ) is sufficient. However, the quadrupole limit cannot describe the dynamics when the perturber is eccentric (e.g., Naoz et al. 2013Naoz et al. , 2011Katz et al. 2011), and the disturbing function needs to be expanded up to the octupole order (a 1 /a 2 ) 3 . With an eccentric outer perturber, the eccentricity of the inner orbit can be excited close to unity. In addition, the inclination of the inner orbit can cross 90 • , even when it starts in a near-coplanar configuration with respect to the outer binary (Li et al. 2014b). This can explain a wide range of astrophysical phenomena (see Naoz 2016 for a review).
While the Kozai-Lidov mechanism and the eccentric Kozai-Lidov mechanism mainly describe the evolution of a test particle perturbed by an outer massive object, the evolution of an outer test particle perturbed by an inner binary has also been studied (e.g., Naoz et al. 2017;Vinson & Chiang 2018;de Ela et al. 2019). In contrast to the inner test particle case, the quadrupole resonance allows the outer test particle's orbit to flip without changing its eccentricity. Higher order resonances can further excite test particle eccentricity and inclination. This has important implications for the dynamical evolution of debris disk surrounding planets in eccentric orbits (e.g., Zanardi et al. 2017).
Beyond the hierarchical limit, non-hierarchical dynamics also has wide applications but has been less explored. Recently, using a large ensemble of N-body simulations, Stone & Leigh 2019 obtained a statistical solution to the chaotic non-hierarchical three-body problem, under the assumption of ergodicity. They found that the nonhierarchical triple interactions almost always lead to a single escaping object and a stable bound binary. In addition, the eccentricity of the surviving binary follows a super thermal distribution.
In mildly hierarchical triples, in which the triple system could still survive, non-secular effects can become important and can enhance the inner binary eccentricity (e.g., Cuk & Burns 2004;Antonini & Perets 2012;Antonini et al. 2014). In particular, Luo et al. 2016 showed that short-time-scale oscillations can accumulate and make the secular results unreliable. They obtained 'corrected double averaging' equations to account for the error in the secular results. Based on this correction, Grishin et al. (2018) obtained analytical results on the maximum eccentricity of the inner binary at the quadruple level.
However, the secular results can still provide a good approximation when the perturber is much less massive than the central object, and the corrected terms at the quadruple level are not sufficient when the three ob-jects are closer to each other (with semi-major axis ratio 0.3). In particular, Gronchi & Milani (1998) studied the dynamics of near earth asteroids, and developed a secular method without expansion of the semi-major axis ratio. More recently, this method has been used to study secular interactions between extreme TNOs and the hypothetical Planet-9 (Beust 2016;Li et al. 2018). These studies mostly focused on the near-coplanar regime, where Planet-9 is located near the ecliptic plane. Nevertheless, a systematic study of dynamics of the mildly-hierarchical systems is missing in the literature, and the effects of an inclined Planet-9 on the extreme TNOs remains poorly explored.
In this paper, we study the secular interactions in mildly-hierarchical systems with either an inner or an outer perturber. We compare the dynamical features to those in the hierarchical limit using the surface of section, which identifies the location of resonant and chaotic regions. We also identify initial conditions which can lead to large eccentricity and inclination variations. Finally, we apply our results to study interactions between objects in the outer solar system and the hypothetical Planet-9. The paper is organized as the following: in §2 we discuss the numerical techniques we adopt in this work, and in §3, we derive the limit where the secular approach is valid. Then, we present the surface of sections in §4, and the maximum eccentricity and inclination variations in §5. In §6, we apply our results to the objects in the outer solar system. Finally, we discuss our results and conclude in §7.

NUMERICAL METHODS
Mildly hierarchical triples can be separated into two binaries (as shown in Figure 1). The closely separated two objects, m 1 and m, form the inner binary, and the outer binary is composed of the farther object, m 2 , orbiting around m. Note that this is different from the generic set up of hierarchical triples, where the outer binary is typically set to originate from the center of mass of the inner binary. This makes little difference here, since the perturbing mass that we consider must be much smaller than the central object for the secular results to be valid.
Here, we use subscript 1 to denote the orbital elements of the inner orbit, and subscript 2 for the outer binary: r 1 and r 2 represent the position vectors of m 1 and m 2 from the central massive object m. In this paper we study the dynamics of a test particle in the triple system, and consider configurations with both inner and outer massive perturber. Thus, depending on the configuration, either m 1 or m 2 could be the test particle. We use parameter α to denote the ratio of the semi-major axis of the test particle with respect to the semi-major axis of the perturber. α < 1 for the inner test particle scenario, and α > 1 for the outer test particle case.
To describe the motion of a test particle under the influence of the central object and a lower mass companion, we follow the expressions below (e.g., Murray & Dermott Figure 1. Configuration of a mildly-hierarchical triple system. It comprises an inner binary, m and m 1 , and an outer binary formed by m 2 orbiting around the central object (m). Orbital elements of the inner binary are denoted by subscript 1, and those of the outer binary are denoted by subscript 2.

2000):¨
Where R i is the disturbing function, r i and r j are the positions of the test-particle and the perturber with respect to the central object respectively, m j is the mass of the perturber, m is the mass of the central object and G is the universal gravitational constant. The disturbing function describes the interaction potential between the test particle and the perturber, as the effect due to the influence of the perturber on the central body. The first term of the disturbing function is called the direct term and the second one is called the indirect term. Indirect terms don't contribute to secular perturbations, as they average to zero (e.g., Murray & Dermott 1999). We can write the secular disturbing function as, which averages out the fast, orbital timescale variations. In the mildly-hierarchical limit, the semi-major axis ratio is no longer a small parameter, so we numerically average the disturbing function instead of expanding it in the ratio of semi-major axes. This allows us to explore regimes which are non-hierarchical. Orbital elements of the test particle are evolved using Hamiltons equations. We use the scale invariant Delaunay conjugate variables to define our phase space. They are composed of conjugate pairs J z − Ω and J − ω, where in the test particle limit J = √ 1 − e 2 and J z = √ 1 − e 2 cos I. The equa-tions of motion are given by: We solve the above set of equations numerically 1 . In non-hierarchical systems, the orbits of the perturber and the test particle can intersect. Although in the unaveraged system the perturber and the test particle may not actually collide, the two bodies may coincide at one or more points during the process of averaging. When they do, the disturbing function and its derivatives have singularities which can be difficult to evaluate. We use the CQUAD integration routine from the GNU Scientific Library(GSL) to perform these integrations. The equations of motion are evolved using GSL's implementation of the Runge-Kutta-Fehlberg (4, 5) method.
Furthermore, Gronchi & Milani (2001) show that when orbits cross (or intersect) the derivatives of the disturbing function are not continuous across the crossing. Extra care needs to be taken when orbits cross. Whenever orbits cross during a single time-step, we abandon the current step and take a specific time-step to land exactly at the intersection. To find the intersection, we use the mutual nodal distance as the independent variable, which allows us to choose a step to reach the intersection (following Henon 1982). This is similar to the approach used by . Instead of using implicit integration schemes, we use forward (and backward) discrete differentiation to calculate derivatives of the disturbing function from the positive (and negative) direction approaching the intersection separately. This provides the equation of motion prior (and post) the intersection.

ACCURACY OF THE SECULAR RESULTS
The double averaging method outlined in Section 2 is limited to low mass perturbers. As the mass of the perturber increases, non-secular oscillatory effects, such as mean motion resonances and evection resonances, become important (Luo et al. 2016;Grishin et al. 2018), and the secular approach is no longer reliable. For instance, in the secular regime the mass of the perturber is just a scaling factor which would only change the dynamical timescale of the system (as shown in Equations 4 and 5). This is not true when the mass of the perturber is high.
To illustrate this, we perform an ensemble of double averaged secular simulations and compare them with Nbody simulations. We focus on the case of an inner test particle under the perturbation from an object on a circular orbit. We run these simulations for three Kozai timescales, t k = (m/m 2 )(a 2 /a 1 ) 3 P , where P is the orbital period of the test particle. We record the maximum eccentricity reached by the test particle. In Figure 2, we plot the maximum eccentricities against the mass of the perturber. The perturber is on a circular orbit, and test particles start on a circular orbit (e 1 = 0) with i 1 = 75 • and Ω 1 = π/4. Three different configurations of the test particle are shown: the panel on the left has a semi-major ratio of α = 0.1, the one in the middle corresponds to α = 0.3, and the one on the right corresponds to α = 0.5. The blue pluses represent the double average results, and the green crosses represent the N-body results.
The maximum eccentricity obtained from the double average results is independent of the perturber's mass, because the masses only determine the dynamical timescale of the system. For the n-body results, as the mass of the perturber increases so do the non-secular effects, causing e max to increase. A purely secular double averaging method fails to capture this effect. Thus, as the perturber become more massive, the maximum eccentricity calculated from the double averaged method (e DA max ) starts to deviate from that calculated using Nbody simulations (e N B max ), as shown in Figure 2. We use black dashed lines to denote the critical perturber mass, the mass of the perturber above which the double averaged results start to deviate from the N-body results. Specifically, the critical mass is defined by the following criteria: As the semi-major axis ratio increases, the non-secular effects of the perturber becomes stronger and the critical mass decreases.
For low values of α, the deviation from the secular results can be explained by single averaged corrections, which are obtained by averaging the disturbing function only over the orbit of the inner test particle (the shortest time scale in the problem). Luo et al. (2016) derived single averaged corrections to the double averaged disturbing function. Using their derivation and the quadrupole order double averaged disturbing function, Grishin et al. (2018) derived single averaged correction (δe SA corr (m 2 , α)) to e DA max . In Figure 2, red dashed lines represent the results by Grishin et al. (2018). While the quadrupole result agree very well with N-body results at α = 0.1, at higher values of α, it fails to agree even for low m 2 . Therefore, we use a semi-analytical method which includes terms up to α 8 to calculateē SA max , while using the same single averaged correction (δe SA corr (m 2 , α)) from Grishin et al. (2018). We focus on an object in the resonant region following a librating trajectory. Thus, ω 1 = π/2 when e 1 = 0 and e 1 = e max . Since j z is constant under the perturbation of objects in circular orbits, we obtain the inclination when the eccentricity assumes it's maximum value: When the orbit of the perturber is circular, the Hamiltonian is independent of Ω, and only terms of even order survive in the Hamiltonian in addition to the the single averaged correction: where the single averaged correction is derived in Eqn. . Maximum eccentricity of the test particle as a function of the mass of the perturber. Left panel corresponds to low α = 0.1, middle panel corresponds to α = 0.3, and right panel corresponds to α = 0.5. Green crosses represent results using N-body simulations, and blue dots represent results using the double averaged secular simulations. The critical perturber masses beyond which the secular results deviate from that of N-body simulations are indicated by the black dashed line. As α increases, the critical mass decreases. Red dashed lines show emax obtained by Grishin et al. (2018) including the single average oscillations, and cyan solid lines represent our results of eemax, calculated using a higher order semi-analytical method (using terms of disturbing function up to α 8 ). This method provides more accurate estimates of emax at higher values of α.
In this case, the maximum eccentricity can be obtained numerically, requiring that energy is conserved: We show the semi-analytical results (cyan lines) in Figure 2. By including the 8 th order, the semi-analytical approach agrees very well with the direct N-body simulations at higher values of α and low values of m 2 . However, the semi-analytical results deviate from the N-body results for larger m 2 when α ∼ 0.5, suggesting higher order single averaged corrections need to be included. In addition, wiggles in the numerical results are mainly due to close encounters, which cannot be captured by the secular approach.
In the end, we derive an analytical expression for the critical perturber mass, below which the double average provides a good approximation, using the results of the maximum eccentricities. When the perturber is on a circular orbit, j z is a constant of motion as a result of the double average over the inner and outer binary orbits. However, single averaged corrections lead to oscillations in j z with maximum fluctuations proportional to SA . Assuming all the deviations between e DA max and e N B max are due to fluctuations in j z and are caused by the single average oscillations, we can obtain the fractional change in the maximum eccentricity due to the single averaged oscillation effects as the following: where f e = 0.01 gives the critical mass of m 2 according to our definition in Eqn. (6), and e SA corr and e SA max in the quadrupole limit are obtained in Grishin et al. (2018) (Eqns. 33 and 38): To the first order in SA and assuming m 2 m 1 , we can get an analytical expression for the critical mass m 2,crit , where f e < C: We compare the expression with N-body simulations in Figure 3, which shows the critical mass as a function of α. The dots represent the N-body results, the dashed line represents the analytical results in Eqn. 14, and the solid red line represents the semi-analytical results based on the 8 th order expansion in semi-major axis ratio ((a in /a out ) 8 ). When the perturber is close to the test particle a 1 /a 2 ∼ 0.5, the perturber needs to be small comparing with the central object (∼ 10 −3.5 that of the central object's mass). When the perturber is farther (a 1 /a 2 ∼ 0.1), the perturber only needs to be lower than ∼ 10 −1.5 times the mass of the central object for the double average to be valid.
Our analytical expression agrees well with the N-body results when α 0.3. When α 0.3, effects of mean motion resonances (P t : P 2 = 1 : 5, 1 : 4, 1 : 3) become important, and this makes the analytical results based on the single averaged method no longer reliable. The semi-analytical results up to the 8 th power in α agree N-body Semi-analytical(n=8) Eqn. 14 Figure 3. Critical mass as a function of the ratio of semi-major axes (α). The blue dots show the results using N-body simulations. We define critical mass as the minimum mass of the perturber above which the maximum eccentricity derived using secular approximation differs from that derived from N-body simulations (Equation 6). As the semi-major axis ratio (α) increases, critical mass decreases. We derived an analytical expression for critical mass (green dashed line, Equation 14), which agrees well with the N-body results when α 0.3. The red solid line represents the semi-analytical results up to the 8 th order in semi-major axis ratio, and it agrees slightly better than the analytical results. The dashed black lines represent the location of the mean motion resonances.
slightly better compared to the analytical expression. When α 0.3, the single averaged correction becomes similar in magnitude to the expansion at the n = 8 order (up to (a 1 /a 2 ) 8 ). Thus, one needs even higher order expansions to obtain more accurate semi-analytical results on the critical mass.

SURFACE OF SECTIONS
For systems with perturbing mass much lower than that of the central object(∼ 10 −3 to 10 −2 m), the secular results using the double average method provide a good approximation to the dynamics. Thus, to better understand the dynamics of mildly hierarchical triples, we analyze the secular results in this section in more detail. In the test particle limit, the secular dynamics can be reduced to two degrees of freedom, and the phase space is four dimensional. It is difficult to visualize the phase space directly and thus, we use surface of sections to characterize the dynamical properties, similar to the approach in e.g., Li et al. 2014a.
We first look at surface of sections in e cos ω − e sin ω space. We perform an ensemble of secular simulations using the method outlined in Section 2. From the trajectories of the test particles evolved in these simulations, we collect all points on the surface Ω = 0 for inner test particle configurations and Ω = π/2 for outer test particle configurations which satisfy the conditionΩ < 0. We choose Ω = π/2 for the outer test particle configurations because the librating trajectories are centered around Ω = π/2 due to the quadrupole resonances (Naoz et al. 2017).
To study how the surfaces change as the system become less hierarchical, we make surfaces for different values of the semi-major axes ratio (α = {0.1, 0.2, 0.3, 0.5, 2, 3, 5}), and we also consider perturbers with different eccentricities (e 2 = {0.2, 0.4, 0.6, 0.8}). We show a few representative panels in Figure 4. See the Appendix for the full set of surfaces. We take the argument of pericenter (ω 2 ) and longitude of ascending node(Ω 2 ) of the perturber to be zero.  (ω) space for e 2 = 0.6. We use the conditions Ω = 0 andΩ < 0 for inner test particle configurations and Ω = π/2 andΩ < 0 for outer test particle configurations to choose our points.
As the systems become less hierarchical, physically allowed regions change, particularly for those with an eccentric perturbers e 2 0.4. For instance, regions inside −π/2 < ω < π/2 can be allowed in the mildly hierarchical configurations with an eccentric perturber, as shown in the top right panel in Figure 4. However, the dynamical features typically have a weak dependence on the semi-major axis ratio in the same dynamically allowed regions. Similar to the hierarchical limit, where resonances have been identified at ω = π/2 and 3π/2 (Kozai 1962;Lidov 1962) in the quadrupole limit, and at ω = π in the octupole limit in surfaces (Li et al. 2014a), we identify resonances at ω = {π/2, π, 3π/2} in the inner test particle case. Overall, the system becomes more chaotic when it is less hierarchical.
For outer test particle configurations (bottom row of Figure 4), we can see resonances at ω = 3π/2 (note that resonances at ω = π/2 can also be seen in some surfaces, e.g., middle panel in e 2 = 0.6, α = 5 in the Appendix). This is also similar to the hierarchical limit. Specifically, ω has no resonances at the quadrupole limit, but there are resonances at ω = π/2 in the octupole order Hamiltonian (Naoz et al. 2017). Note that Naoz et al. (2017) identified resonances at π/2, because they useΩ > 0 as their condition to make the surfaces. To focus on the low inclination dynamics, we useΩ < 0 instead, since the inclination of points withΩ > 0 are typically over 90 • . Similar to the inner test particle cases, as the system become less hierarchical, there are generally more chaotic regions at higher values of e 2 . The chaotic regions result from overlapping of resonances, and we can see higher order resonances embedded in chaotic regions (See the panel for α = 2 , e 2 = 0.2).
The colors in the surfaces represent inclination of the particles, which oscillate as the eccentricities vary. For example, trajectories librating around ω = {π/2, 3π/2} undergo inclination variations that can flip over 90 • . Different from the hierarchical limit, orbits can flip from a near coplanar configuration in a near circular orbit under the influence of an eccentric outer perturber. We will discuss this in more detail in Section 5.
We now look at surface of sections in i cos(Ω)−i sin(Ω).
We plot surface of sections for systems with the same set of configurations as discussed above. In Figure 5 we show surfaces for α = {0.1, 0.3, 0.5, 2, 3, 5} and e 2 = {0.6} (see the Appendix for the rest of surfaces). We chose the surface ω = 0 with the conditionω > 0 to make these plots. Color in this case represents the eccentricities of the test particles. Similar to the surfaces in e cos(ω) − e sin(ω), physically allowed regions change for mildly hierarchical configurations with highly eccentric outer perturbers. In particular, additional resonant regions show up around Ω = 0 with high inclinations and moderate eccentricities (e.g., top right panel in Figure 5), under perturbations of an eccentric outer object (e 2 0.4) with α 0.3. At low values of α, resonances around Ω = π lead to high eccentricity excitation as the orbits flip cross 90 • . In addition, test particle's orbits circulate when the inclination is close to 0 • and 180 • due to an eccentric perturber (e.g., top left panel of Figure 5). The orbits librate for a wide range of inclinations in between the coplanar and counter orbiting configurations. To illustrate the flips of the orbit, we marked the circle with a radius of 90 • . Many librating orbits cross over 90 • at high eccentricities, similar to the hierarchical limit. When the systems become less hierarchical, eccentricity typically oscillates faster than inclination, and remains low when ω=0. Thus, we find particles with moderate or low eccentricities near 90 • on the surfaces (e.g., top right panel of Figure 5). However, the orbits still flip at relatively high eccentricities (e > 0.8), and the eccentricity when the orbit flips do not show a clear trend as a function of α.
The resonant regions in the plane of i cos Ω − i sin Ω are again similar to that in the hierarchical limit. We identify resonances at Ω = {0, π, π/2 and 3π/2} in inner test particle configurations. Note that the resonant regions around π/2 and 3π/2 exist in the surface of e.g., α = 0.3, e 2 = 0.2 near i ∼ 90 • (as shown in the the fifth column in Figure 17 in the Appendix). This is similar to the octupole order, where Li et al. (2014a) find resonances at Ω = {0, π}. We could also identify higher order resonances with high (α) values embedded in the chaotic sea. Similar to the surfaces in e cos ω − e sin ω, chaotic regions are more common at higher α values. For outer test particle configurations, we can identify resonances at Ω = {π/2, π, 3π/2}. This is consistent with the hierarchi-cal limit where at the quadrupole order, Ω is the resonant angle which librates around Ω = π/2 (Naoz et al. 2017). Similar to the surfaces shown in e cos ω − e sin ω, the size of the chaotic regions increases as the systems become less hierarchical.

ECCENTRICITY AND INCLINATION EXCITATION
The surface of section plots in Section 4 demonstrate that the eccentricity and inclination of the test particles can have large amplitude variations, and this may have important implications for astrophysical systems. In this section, we consider the eccentricity and inclination excitation in detail. First, to illustrate how the eccentricity excitation changes as the system become less hierarchical, we show the maximum eccentricity of the test particle as a function of the semi-major axes ratio (α) in Figure 6 for a circular perturber with a mass ratio m 2 /m = 3 × 10 −4 .  Figure 6. Maximum eccentricity as a function of the semi-major axis ratio (α) for the inner test particle configuration with a circular perturber. The test particle is initially on a circular orbit with i = π/3, ω = π/2 and Ω = 0. Results from N-body simulations are shown as blue dots whereas the results from secular simulations are shown as yellow dashed lines. It can be seen that results from secular simulations agree very well with N-body simulations for all values of α. The quadruple level result is shown for comparison. We also show results from the semi-analytical method which uses terms in the disturbing function up to the order α 8 . While it agrees with N-body simulation till α = 0.35, at higher values of α this method over predicts emax. Figure 6 shows that as the semi-major axis ratio increases so does the maximum eccentricity. While the quadruple order expression for the maximum eccentricity (e max = 1 − 5 3 cos 2 i 0 ) is accurate only for α < 0.05, the semi-analytical results at the 8 th order of expansion outlined in Section 3 agree more closely with N-body simulations for higher values of α. For α > 0.3, this approach also ceases to agree with N-body results and predicts a higher value for maximum eccentricity. This indicates the importance of even higher order terms of the disturbing function in this regime. The double average results agree well with the N-body results with the low mass perturber. With a circular perturber, the maximum eccentricity can be quite large, but it is mostly well below unity.
When the orbit of the perturber becomes eccentric, maximum eccentricity of the test particle becomes more complicated. Previous studies on hierarchical systems at the octupole order showed that the eccentricity of the test particle can be excited to values close to unity, and inclination can be excited to values beyond 90 • (Li et al. 2014b). Orbits can be flipped (inclination crosses over 90 • ) in three different scenarios 1) low initial eccentricity and high initial inclination(i > 40 • and i < 140 • ) a scenario similar to the standard Kozai-Lidov mechanism. 2) high initial ecentricity and low inclination and 3) medium eccentricity and high inclination. When the perturber's orbit is circular, the dynamics is mainly dominated by quadrupole resonances (as discussed in §3), while when the orbit of the perturber is eccentric, both octupole and quadrupole resonances become important (See e.g., Li et al. 2014a). We use a numerical approach to study the maximum eccentricity excitation under the influence of an eccentric perturber here. In Figure 7 we show the maximum eccentricity as a function of semi-major axis ratio(α) for test particles initially aligned (upper panel) and anti-aligned (lower panel) with respect to an eccentric perturber. The perturber (with m 2 /m = 3 × 10 −4 ) is on an orbit with an eccentricity of 0.6, and test particles start on an orbit with e = 0.4 and i = 5 • .
We plot the maximum eccentricity for both inner (left panels) and outer (right panels) test particle configurations. In addition, we consider initially aligned configurations (top panels), where the test particle pericenter starts aligned with that of the perturber, as well as the anti-aligned configurations (bottom panels). The maximum eccentricity increases as the system becomes less hierarchical for both inner and outer test particles. Except for when the test particles are close to mean-motion resonances, double averaged simulations agree with Nbody simulations across all values of α shown here. In contrast, the octupole level hamiltonian provides a good approximation only for α < 0.25 and α > 2 in the initially aligned configuration and for α < 0.1 and α > 8 for the initially anti-aligned configuration.
To investigate the eccentricity and inclination variations for different initial eccentricities, we performed an ensemble of secular simulations. We choose the following initial conditions for the test particles: α ∈ [0.1, 0.5] for inner test particle configurations and α ∈ [2, 10] for outer test particle configurations, e 0 ∈ {0, 0.9}, i 0 = 5 • along with three different initial orientations 0 = {0 • , 180 • , random} where is longitude of pericenter of the test particle and e 0 is the initial eccentricity of the test particle. For the perturber, we choose e 2 = 0.6, ω 2 = 0, Ω 2 = 0 and i 2 = 0. We evolve the inner-test particle systems for a maximum of 10t k and outer test-particle systems for 10t p and record the maximum eccentricity and inclination reached in each run. t k is the Kozai timescale, and t p is the precession time scale given by: where m is the mass of the central object, m pert is the mass of the perturber, a 1 is the semi major axis of the test particle, a pert is the semi-major axis of the perturber and P is orbital period of the test particle.
In Figure 8, we show the maximum eccentricity in the plane of the initial eccentricity(e 0 ) and semi-major axis ratio(α). The upper panels study the case of an inner test particle. The left panel starts with the aligned configuration, and the middle panel starts with the anti-aligned one. Similar to the hierarchical limit, it is generally easier to excite the eccentricities as compared to initially aligned ones. In the hierarchical limit at the octupole order, an analytical expression can be derived to predict the flip of the inner orbit (Li et al. 2014b), and we include the criterion using dashed lines in the figures. It shows that the criterion can still be valid up to α ∼ 0.1 for the anti-aligned configuration.
The criterion does not predict the flips when the system is initially aligned, since the orbit cannot be flipped starting in the aligned configuration in the hierarchical limit (α 0.1) at the octupole order (Li et al. 2014b). The right panel corresponds to randomly aligned initial configurations ( = random), and it shows the impor- . Maximum eccentricity excitation for mildly hierarchical triples. x axis represents the ratio of semi-major axes and y-axis represents the initial eccentricity of the test particle. The dashed lines in the upper panels illustrate the orbit-flip criteria from (Li et al. 2014b). Test particles start at an inclination of 5 • . Perturber is on an orbit with e = 0.6, i = 0, ω = 0 and Ω = 0. Color represents the maximum eccentricity reached in 10 Kozai timescales. The upper panels have an inner test particle, the lower panels have an outer test particle.
tance of the initial orientation. Specifically, the eccentricity can remain low for specific values of even at high values of α. This can be attributed to the resonant regions which can be seen in the surface of sections for high α configurations.
The lower panel corresponds to the outer test particle configurations. In contrast to the inner test particle cases, the eccentricity is rarely excited close to unity, which would lead to the cross of orbits for the outer test particles. Similar to the inner test particles, the eccentricity is more likely to be excited when the orbits are anti-aligned, and when the system is less hierarchical. Only in anti-aligned configurations can the near circular orbit excite to high eccentricities, even in the mildly hierarchical configurations.
It is also apparent that large inclination variations are accompanied by large eccentricity excitations. To illustrate the variations in inclination, we color code the maximum inclination as a function of the initial eccentricity (e 0 ) and semi-major axes ratio (α) in Figure 9. The flip criterion in the octupole limit by Li et al. (2014b) is included as black dashed lines, and it is consistent with the double averaged results for α < 0.1 and = 180. At higher values of α, the orbits can flip with lower initial eccentricities as compared to the flip criteria. This is consistent with the results shown in Figure 7. In particular, a near circular inner test particle orbit could flip due to the perturbation of an eccentric outer object starting anti-aligned with the test particle with α ∼ 0.4.
For outer test particles (lower panel), orbit flips are rare starting from a near co-planar configuration unless the orbit starts mildly hierarchical (α 10) and in an anti-aligned configuration. Eccentricities are excited very close to unity during the flips, and frequently causes the cross of orbits. Similar to the eccentricity excitation, the flips are only seen when the orbits are anti-aligned.

APPLICATIONS TO OUTER SOLAR SYSTEM
Mildly hierarchical triples are stable when the perturbing mass is much less massive than the central object, such as our own Solar System's extreme trans-Neptunian objects (eTNOs, with large semi-major axes and high eccentricity) under the perturbation of a possible undetected planet residing far from the Sun (Planet-9). In this section, we apply our understanding of the dynamics of the mildly hierarchical triples to the Solar System, and constrain the properties of the possible undetected Planet-9.
Planet-9 is proposed to explain mysterious features in the orbital distribution of extreme trans-neptunian ob- . Maximum inclination excitation for mildly hierarchical triples. x axis represents the ratio of semi-major axes and y-axis represents the initial eccentricity of the test particle. Test particles start at an inclination of 5 • . Perturber is on an orbit with e = 0.6, i = 0, ω = 0and Ω = 0. Color represents the maximum inclination reached in 10 Kozai timescales. The upper panels have an inner test particle, the lower panels have an outer test particle.
ject (objects that lie outside the orbit of Neptune with large semi-major axes 150 AU) (Sheppard & Trujillo 2016;Batygin & Brown 2016). In particular, they show clustering in the longitude of the ascending node, the argument of pericenter and their orbital planes. Many studies have shown that the alignment of the orbits is not due to selection biases (de La Fuente Marcos & de La Fuente Marcos 2014; Gomes et al. 2015;Batygin & Brown 2016), although Shankman et al. (2017) demonstrated that the Outer Solar System Origins Survey (OS-SOS; Bannister et al. 2016) contains nonintuitive biases for the detection of TNOs that lead to apparent clustering of orbital angles in their data, and the angular elements of the distant TNOs are consistent with uniform distribution (Bannister et al. 2018). Ongoing observational TNO surveys will provide a better understanding of the architecture of the outer solar system and the details of the TNO clustering, and it's important to obtain a better understanding on the dynamical origin of possible clusterings of eTNOs.
The proposed Planet-9 is on an eccentric orbit with a large semi-major axis which can perturb the orbits of the eTNOs and explain the clustering of the eTNOs. The underlying dynamics of the interactions between the eT-NOs and Planet-9 is rich, and it is found that secular dynamics plays an important role leading to the cluster-ing of the eTNO orbits Saillenfest & Lari 2017;Li et al. 2018). To illustrate the dominance of the secular resonances, we compare our secular results with the those we got from N-body simulations below.
In Figure 10, we show the projections for particles orbiting around the Sun at a = 250AU, under the perturbation of a Planet-9 of mass m 9 = 10M ⊕ , with semimajor axis a 9 = 500AU and eccentricity of e 9 = 0.2. To illustrate the inclination variations and the orbital plane clustering, we show the projections in the plane of i sin(Ω) and i cos(Ω). The projections are obtained in the same way as the surface of sections in Section 4. We choose the initial condition of the test particle based the results of the surface of section, so that they start with the same energy. The color represents the eccentricity of the particles. The top-left panel shows the results from N-body simulations, which do not include the effects of the giant planets. Without the effects of the giant planets, the results of the three-body secular results largely resembles that of the N-body simulations (compare with the last panel for α = 0.5 in Figure 17), except that parts of the libration region around Ω ∼ 180 • become chaotic. This is due to orbital crossing of eTNOs with Planet-9, which is more often when the orbits' pericenters are anti-aligned with repect to each other and lead to close scattering of eTNOs. The close scattering of eT- Bottom Left panel: Projection with the J 2 potential in mutual coordinates. In mutual coordinates, the location of resonance is restored to Ω = 0and π. Bottom Right panel: Projection evaluated in mutual coordinates using N-body code with Neptune as an active object and a J 2 potential for other planets of the solar system.
NOs by Planet-9 cannot be reproduced by the secular results, and it often leads to chaotic behaviors.
In the top-right panel, we show N-body results which include the effects of the giant planets as a J 2 potential. The J 2 potential causes precession in Ω, which leads to more frequent orbital crossing. This suppresses the dominance of resonances in Ω and leads to larger chaotic regions. However, eccentricity excitation near Ω ∼ 0 • due to the secular resonances can still be found in the N-body results, and this can produce the clustering of the eTNO orbits with high eccentricity. In addition, with the J 2 potential, the orbit of the Planet Nine is no longer constant. It undergoes precession and this makes a slight shift in the resonant region at Ω = 180 • .
In the secular results, the orbital elements are obtained relative to the orbit of the outer planet, however the orbital elements are relative to the ecliptic plane in the N-body results. While the orbital elements of the pertuber in non-hierarchical triples are constant, Planet-9's orbital elements are not. Neptune and other inner planets exert torques which causes Planet-9's orbit to precess. Thus, to better compare our secular results for mildly hierarchical triples with N-body simulations of eTNOs with Planet-9 and the J 2 potential, we use mutual coordinates, which is set up as the following. We choose our reference plane and reference direction such that inclination(i 9 ), argument of pericenter(ω 9 ) and longitude of ascending node(Ω 9 ) of the Planet-9 are set to 0, following the example of Gronchi (2002). The results based on the mutual coordinates are shown in in the bottom-left panel, and we can see that the resonance at Ω = 0 is restored in mutual coordinates. The precession due to the J 2 potential allows the orbit of the perturber and that of the test particles to cross more frequently. This causes the system to become significantly more chaotic as shown in a 9 ∈ {300, 700, 1000, 1400}AU e 9 ∈ {0.2, 0.4, 0.6, 0.8} m 9 ∈ {5, 10, 20, 30}M ⊕ i 9 ∈ {0, 30 • , 60 • , 90 • , 120 • , 150 • , 180 • } ω 9 = 0, Ω 9 = 0 Table 1 Planet-9 orbital elements. We note that some of the parameters above could be ruled out by observational constraints, but we included them here in order to analyze the dynamics in a broad range of parameter space.   Finally, we include a point mass Neptune in our Nbody simulation (as an active body instead of a J 2 term) and show the results in the bottom right panel of Figure  10. Including a point mass Neptune, the secular dynamical features get further suppressed. However, we could still see some clustering of the orbit in Ω with moderate eccentricity at high inclination around ∼ 180 • , due to the secular resonances.
We next ran an ensemble of N-body simulations and constrain the properties of Planet-9. Multiple studies have attempted to constrain the orbital parameters of a potential Planet-9 using N-body simulations. In their original paper Batygin & Brown (2016) reported that Planet-9 could be a 10M ⊕ (m 9 ) planet on an orbit with a semi-major axis(a 9 ) of 700 AU, an eccentricity of 0.6(e 9 ) and an inclination(i 9 ) of 30 • . These values were updated in a recent paper, where Batygin et al. (2019) find that m 9 ∼ (5, 10)M⊕, a 9 ∼ (400, 800)AU, e 9 ∼ (0.2, 0.5) and i 9 ∼ (15 − 25) • . The previous studies have focused on configurations of Planet-9 with inclination less that 35 • , but it is possible that Planet-9 lies largely out of the plane. Thus, we sample the inclination of Planet-9 over a wide range of inclinations here.
In Table 1 we list the configurations of Planet-9 we use for the N-body simulations. We use 1000 test particles to model the Kuiper belt. Initial conditions for these particles are listed in Table 2. We also include Neptune as a point particle in our simulations. Other giant planets in the solar system are modeled using a J 2 potential. We use the hybrid symplectic integration method in MER-CURY package (Chambers 1999) to do our simulations.
We first calculate the critical semi-major axis a c , beyond which we expect to see clustering. We do this by dividing the simulation data into bins in semi-major axis and mark the bin where the distribution of stops being uniform. In our analysis we look for clustering in the longitude of pericenter by fitting the distribution of ∆ (= − 9 ) for a > a c using the following probability density function: (17) where κ is a measure of the strength of the clustering, where a larger value of κ corresponds to less clustering. µ and σ are the mean and the standard deviation of the normal distribution, which denote the location and the strength of the cluster. We calculate the values of κ, µ and σ for all the simulations we do, and we use the Kolmogorov-Smirnov test to check the goodness of fit. We select the simulations which pass the test in the following analysis.
In Figure 11 we show an example. In the left panel we plot ∆ as a function of the semi-major axis ratio for test particles. In this simulation, Planet-9 has a semimajor axis of 700 AU with an eccentricity of 0.6 and mass of 10M ⊕ . We show the critical semi-major axis ratio in black at 260 AU for this simulation. The points are sampled from the trajectories of the test-particles. Points after 2 billion years of evolution are chosen at an interval of 1 million years. On the right we show the show the distribution of ∆ for a > a c (where a c = 260 AU), as well as the fit. For this fit we find µ = 177.7 • , κ = 0.321 and σ = 23.70 • .
We use association rule analysis to determine which configurations of Planet-9 leads to clustering of the test particles (lower κ). More specifically, association rule analysis is a machine learning method used to find relationships between variables in a data set. In the context of this study, each N-body simulation corresponds to a transaction which comprises of the orbital parameters of Planet-9 and the value of the fitting parameter κ, calculated using the method described above. To use the association rule analysis we discretize the parameter κ as the following: For example, for the simulation with Planet-9 orbital elements {a 9 = 300AU, e 9 = 0.6, i 9 = 0, m 9 = 5} we got κ = 0.55. This translates to a transaction t = {"a 9 = 300AU", "e 9 = 0.6", "i 9 = 0", "m 9 = 5" ,"medium cluster"}. Association rules are denoted by the notation: A → B, where A and B are sets of items. Two quantities are used to specify which rules are interesting: support (s) and confidence (c). Support indicates how frequently a set appears in a data set, and confidence indicates how often the rule is valid. For a given rule A → B, they are: where σ(A) is the count of the set A (i.e. number of times set A appears in the data set) and |D| is the number of transactions in the data set. In this case, |D| is the number of simulations which satisfy the fitting criteria for κ. Rules with higher support are more widely applicable. Confidence measures the reliability of a rule, and rules with higher confidence are more reliable.
This rule tells us that irrespective of what the mass of Planet-9 is, configurations with a 9 = 1000AU, e 9 = 0.8 and i 9 = 90 • lead to high clustering in .
Other rules with high confidences are listed below: The specific values of the support indicate the size of the parameter space that lead to the clustering based on the rules. Thus, rules with fewer Planet-9 parameters (e.g., the third and the fourth rule above) tend to have higher support. While rules with higher support can be applied widely, they tend to have lower confidence, which limits their validity. It should be noted that confidence and support are calculated only for simulations in which ∆ distribution was a good fit for the test distribution we defined above. Of the 448 simulations we ran, only 276 fit this criteria. It is intriguing that the high inclination Planet-9 with large semi-major axis (a 9 ∼1000-1400 AU) on an eccentric orbit(e 9 ∼ 0.6 − 0.8) could lead to strong clustering in . This is due to the long secular oscillation timescale of . During the 4 Gyr simulation, starts to converge to ∼ 170 • and leads to the clustering.
With high inclinations, the cluster in may not correspond to clustering in the orbital directions geometrically. To study the geometrical clustering, we find the spread in the pericenter directions directly, following the approach by Millholland & Laughlin 2019. We select particles with pericenter distance between 30 and 100AU and semi-major axis above 250AU for this analysis, since the scattering with the giant planets becomes important when r p < 30AU, and it is challenging to detect the eT-NOs with r p > 100AU. Below a = 250AU, the clustering Color represents the Kozai timescale for a particle with a = 250AU. It shows that near coplanar, polar (∼ 90 • ) and retrograde Planet-9 with short t K ∼ 1Gyr could lead to the clustering of the pericenter directions, and both near coplanar and counter-orbiting Planet-9 with short t K ∼ 1Gyr could produce low average eTNO inclination.
is weak due to the fast precession in the J 2 potential of the giant planets. We first estimate the average pericenter direction as the following: whereê i is the eccentricity unit vector of i th test particle, and N is the number of test particles which survived 4 billion years of evolution residing in the selected parameter region (30 < r p < 100AU, a > 250AU). Then, we calculate β i , which is the angle between the pericenter orientation of the i th test particle and the average pericenter direction. As a measure of the pericenter clustering, we calculate the average values of β i denoted by β : Since the particles start in the same plane, where β i is uniformly distributed between 0 • and 180 • , the initial β is ∼ 90 • . Using the same approach for the observed sample using data from the Minor Planet Center 2 , there are currently 19 eTNOs with β obs = 60 • , lying with low inclination i obs = 18 • .
In the upper panel of Figure 12, we plot β against inclination of Planet-9. To ensure the results to be statistically robust, we only select simulations with more than 50 particles left in the selected region (30 < r p < 100AU, a > 250AU). Color denotes the Kozai timescale due to the perturbations induced by Planet-9 (Eqn. 15) for an eTNO with semi-major axis of 250AU. It shows that near coplanar, polar and counter orbiting Planet-9 all could lead to the clustering in the pericenter directions, with short dynamical timecales t K 1Gyr. In addition, Planet-9 with shorter Kozai timescales could also lead to stronger clustering. The retrograde orbiting Planet-9 with the same dynamical timescales lead to the same results as that of the prograde cases. This is because the secular dynamics is independence on the direction of the motion of Planet-9.
When Planet-9 is highly inclined and retrograde (i 9 ∼ 150 • ), the massive Planet-9 (∼ 30M ⊕ ) with very short dynamical timescale t K ∼ 10Myr (a 9 ∼ 300AU), could also lead to clustering of the pericenters. However, this feature is missing in the prograde scenario, since not enough particles are left in the selected region as the prograde Planet-9 is more disruptive due to the non-secular effects. Nevertheless, Planet-9 in this parameter space are unlikely since the high mass Planet-9 with small and eccentric orbits can mostly be ruled out observationally.
On the lower panel, we plot the average inclination of the selected test particles (30 < r p < 100AU, a > 250AU) versus the inclination of Planet-9. We find the the average inclination can further constrain the inclination of Planet-9. Although both the polar and the planar Planet-9 could produce the clustering, the average inclination produced by the polar Planet-9 is much higher. The current observation of the eTNOs with low average inclination orbits (∼ 18 • ) favors the near co-planar and counter orbiting Planet-9. We note that the observational surveys are not complete for the high inclination region. The spread of eTNO inclinations due to a polar Planet-9 is large (∼ 40 • ), and selecting only those below 40 • still shows strong clustering. Thus, current observations cannot yet rule out a polar Planet-9, and future observations on the high inclination eTNOs are needed to better constrain the inclination of Planet-9.

CONCLUSIONS
In this paper we study the secular dynamics of mildlyhierarchical triple systems, and apply it to the dynamics of eTNOs in the outer Solar System. In contrast to the hierarchical limit, where one could use a perturbative approach to analyze the dynamics analytically by expanding in semi-major axes ratio, the dynamics can only be investigated numerically for mildly hierarchical systems. Thus, we developed code to numerically calculate the double averaged Hamiltonian, and to evolve the system to study its secular dynamics. The code is now publicly available (https://github.com/bhareeshg/gda3bd).
In mildly hierarchical systems, the perturbation of a massive object can lead to oscillations that render the secular double averaged results unreliable. However, the oscillation amplitude is small and secular results still provide a good approximation when the perturber mass is low. To illustrate this, we calculate the maximum eccentricity induced by the secular interactions and compare that with the N-body simulation in Section 3. We obtained an analytical expression for the critical mass, below which the secular results are reliable (Eqn. 14) for systems with a circular outer perturber. When the semi-major axes ratio is small α ∼ 0.1, the critical mass is m 2 /m 1 ∼ 0.05. With higher values of semi-major axes ratio (α ∼ 0.5), secular results hold for m 2 /m 1 < 10 −3.5 . The critical mass decays with α −3/2 .
In addition, we make surface of section plots to study the dynamics of mildly-hierarchical triples and compare these with the hierarchical limit. We find that the secular dynamics have a weak dependence on the semi-major axis ratio (α) if the perturber does not have a high eccentricity (e 2 0.8). For inner test particles, we find secular resonances at Ω = {0, π, π/2, 3π/2} and ω = {π, π/2, 3π/2}, and for outer test particle we find secular resonances at Ω = {π, π/2, 3π/2} and ω = {π/2, 3π/2}. This is similar to the hierarchical limit. The chaotic regions increase as the systems become less hierarchical. In contrast to the hierarchical limit, test particle orbits can flip from a near co-planar configuration starting with low eccentricities.
We then study the eccentricity and the inclination excitation of the test particles inside the mildlyhierarchical triples. For circular perturbers, we obtain semi-analytical results for the maximum eccentricity with an 8 th order expansion in the semi-major axes ratio, which agrees with the N-body results up to semi-major axes ratio of α ∼ 0.3 (section 3). For an eccentric perturber, double averaged results and N-body results agree well with each other, while octupole results typically give higher e max when α 0.1. We then perform an ensemble of secular simulations, and find that the eccentricity is easily excited for initially anti-aligned systems and for less hierarchical systems.
Finally, we apply our results to objects in the outer solar system. The secular dynamical features resemble that of the N-body results for the eTNOs under the perturbation of the undetected outer planet (Planet-9). Using surface of section, we find that the secular resonances can lead to clustering in eTNO orbits. Next, we perform an ensemble of N-body simulations to model the interactions between Planet-9 and eTNOs. In addition to a low inclination Planet-9 as identified in the literature (e.g., Batygin et al. 2019), we find that a polar (i 9 ∼ 90 • ) and a counter orbiting (i 9 ∼ 180 • ) Planet-9 with Kozai timescales of 1Gyr (for an eTNO with a = 250) could also produce strong clustering. The high inclination or counter orbiting Planet-9 could post challenges on the formation of such a highly inclined wide orbit planet, which may indicate that Planet-9 could be a captured planet from another star (Li & Adams 2016). The nearcoplanar and the counter orbiting Planet-9 could lead to low inclination eTNOs (average inclination of ∼ 20 • ), while the polar one leads to higher inclination ∼ 60 • eT-NOs with a large spread σ i ∼ 40 • . Future observations of the eTNOs with high inclinations could better constrain the Planet-9 inclination.
ACKNOWLEDGEMENT This work was supported by NASA grant 80NSSC20K0641 and 80NSSC20K0522.
For a given configuration of the system, only a finite range of Hamiltonian is physically allowed. To select the energy levels used to make the surfaces in e − ω space, we sample the values of H on a grid in e cos(ω) − e sin(ω) with Ω = 0. Then, we collect the H values at the grid points, and choose the mean ( H grids ) as well as one and two standard deviations away from the mean ( H grids ± σ Hgrids and H grids ± 2σ Hgrids ) to be the energy levels in the surfaces. The different columns in Figure 13-20 correspond to the different energy levels. We follow a similar procedure to select energy levels for surfaces in i − Ω space.
By looking at different surface of sections in a given row (with the same semi-major axes ratio), it is clear that the dynamics depends sensitively on the value of the Hamiltonian. Eccentricities are generally excited to higher values when the energy is low for inner test particles, while it's the opposite when the test particle is on the outer orbit. Various features of these surfaces are presented in detail in Section 4.