# The disturbing function for polar Centaurs and transneptunian objects

F. Namouni, M. Morais

► To cite this version:

F. Namouni, M. Morais. The disturbing function for polar Centaurs and transneptunian objects. Monthly Notices of the Royal Astronomical Society, 2017, 471 (2), pp.2097-2110. 10.1093/mn-ras/stx1714. hal-02398843

HAL Id: hal-02398843

https://hal.science/hal-02398843v1

Submitted on 13 Aug 2025

**HAL** is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire **HAL**, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# The disturbing function for polar Centaurs and transneptunian objects

F. Namouni<fnref n="1" />* and M. H. M. Morais<fnref n="2" />*

Accepted 2017 July 6. Received 2017 July 6; in original form 2017 April 12

## ABSTRACT

The classical disturbing function of the three-body problem is based on an expansion of the gravitational interaction in the vicinity of nearly coplanar orbits. Consequently, it is not suitable for the identification and study of resonances of the Centaurs and transneptunian objects on nearly polar orbits with the Solar system planets. Here, we provide a series expansion algorithm of the gravitational interaction in the vicinity of polar orbits and produce explicitly the disturbing function to fourth order in eccentricity and inclination cosine. The properties of the polar series differ significantly from those of the classical disturbing function: the polar series can model any resonance, as the expansion order is not related to the resonance order. The powers of eccentricity and inclination of the force amplitude of a *p*:*q* resonance do not depend on the value of the resonance order |*p* − *q*| but only on its parity. Thus, all even resonance order eccentricity amplitudes are ∝*e*<fnref n="2" /> and odd ones ∝*e* to lowest order in eccentricity *e*. With the new findings on the structure of the polar disturbing function and the possible resonant critical arguments, we illustrate the dynamics of the polar resonances 1:3, 3:1, 2:9 and 7:9 where transneptunian object 471325 could currently be locked.

**Key words:** celestial mechanics – comets: general – Kuiper belt: general – minor planets, asteroids: general – Oort Cloud.

## 1 INTRODUCTION

The increasing detections of Centaurs and transneptunian objects (TNOs) on nearly polar orbits (Gladman et al. 2009; Chen et al. 2016) raises the question of their origin and relationship to the Solar system planets. Amongst the dynamical processes that govern the evolution of such objects are mean motion resonances. In this context, it was shown recently, through intensive numerical simulations, that mean motion resonances are efficient at polar orbit capture (Namouni & Morais 2015, 2017). It is, therefore, important to have a thorough understanding of the processes of resonance crossing and capture for nearly polar Centaurs and TNOs so that we can identify the pathways that led to such orbits and ultimately uncover their origin.

Identifying a mean motion resonance for a Centaur or a TNO with a Solar system planet is a fundamentally simple task. One has to search for angle combinations of the form $φ = qλ - pλ' - kω + (p - q + k)Ω$ that can be stationary or oscillating around the equilibrium defined by the condition $φ = 0$. In the previous expression, $λ$ and $λ'$ are, respectively, the mean longitudes of the object and the planet, $ω$ and $Ω$ are, respectively, the object's longitudes of perihelion and ascending node, and *p*, *q* and *k* are integer coefficients. In the angle combination $φ$, we ignored

the planet's perihelion and node because the solar system planets' eccentricities and inclinations with respect to the invariable plane are small. As the number of integer combinations is infinite, one usually seeks and checks only the strongest resonances: those with a force amplitude that implies a sizable resonance width within which to capture the Centaur or TNO. This choice is as reasonable as it is rewarding provided that one remembers that the force amplitudes associated with a candidate resonance $φ$ are obtained from the classical disturbing function that is an expansion in powers of eccentricity and inclination of the planet–object gravitational interaction for nearly circular and coplanar orbits. Thus, our intuition regarding the angle combination $φ$ and its dynamical suitability for resonance is based on the assumption of near-coplanarity. It is the object of this work to remedy this shortcoming in the dynamical analysis of polar Centaurs and TNOs by deriving a disturbing function for nearly polar orbits and studying the properties of its force amplitudes.

The history of the classical disturbing function is intertwined with that of celestial mechanics. For a historical perspective, we refer the reader to Brouwer & Clemence (1961). For the purposes of this work, we note that the disturbing function of the three-body problem takes two different forms. The first form is a power series in terms of eccentricity *e* and $sin<fnref n="2" />(I/2)$, where *I* is the inclination. This form therefore assumes that the object's orbit does not depart significantly from prograde coplanar motion. It is used widely to study the formation and dynamics of planetary systems,

<footnote n="1">¹Université Côte d'Azur, CNRS, Observatoire de la Côte d'Azur, CS 34229, F-06304 Nice, France</footnote>

<footnote n="2">²Universidade Estadual Paulista (UNESP), Instituto de Geociências e Ciências Exatas, Av. 24-A, 1515 13506-900 Rio Claro, SP, Brazil</footnote>

<footnote>* E-mail: namouni@oca.eu (FN); helena.morais@rc.unesp.br (MHMM)</footnote>
the formation and evolution of planetary rings and the formation and resonance capture of planetary satellite systems (Murray & Dermott 1999; Ellis & Murray 2000). The second form is a power expansion in terms of the ratio $\alpha = a/a'$ where a and a' are the semimajor axes of the object and planet, respectively. This form is used mainly to study the dynamics of artificial and natural planetary satellites that have large inclinations as they are influenced by the distant Sun or the Moon. It was recently revisited for the study of the secular evolution of hierarchical planetary systems (Laskar & Boué 2010). In the case of Centaurs and TNOs, the second form is not particularly useful as such objects can be quite close to the planets' orbits but unlike satellites they revolve around the Sun not the planet. Reasonable order expansions with respect to $\alpha$ cannot model the dynamics when the semimajor axes' ratio does not satisfy $\alpha \ll 1$. We will therefore seek a disturbing function that is not expanded with respect to $\alpha$ but is written as a power series of eccentricity and some function of the inclination that vanishes if the object's orbit is exactly polar. We find in Section 2 that the natural function is simply $\cos I$. The classical disturbing function and its zero reference inclination can also be transformed into a disturbing function for nearly coplanar retrograde orbits, that is with $180^\circ$ reference inclination, to study the dynamics of retrograde resonances (Morais & Namouni 2013b). The retrograde disturbing function helped to identify the first Centaurs and Damocloids in retrograde resonance with Jupiter and Saturn (Morais & Namouni 2013a).

The plan of the paper is as follows. In Section 2, we write down the explicit steps of the literal expansion of the gravitational interaction of a planet with a particle on a nearly polar orbit in powers of eccentricity and inclination cosine. The reader who is not interested in the details of the expansion algorithm can skip this part and find the resulting disturbing function in Section 3. The properties of the polar disturbing function are compared to those of the classical disturbing function of nearly coplanar prograde orbits as well as that of nearly coplanar retrograde orbits in Section 4. The validity domain of the polar disturbing function is linked to secular evolution and discussed in Section 5. Examples of polar resonances are found in Section 6. Section 7 contains concluding remarks.

## 2 LITERAL EXPANSION FOR NEARLY POLAR ORBITS

Consider a test particle of negligible mass that moves under the gravitational influence of the sun of mass $M_*$ and a planet of mass $m' \ll M_*$. The motion of $m'$ with respect to $M_*$ is a circular orbit of radius $a'$ and longitude angle $\lambda'$. The reference plane is defined by the sun–planet orbit. The test particle’s osculating Keplerian orbit with respect to $M_*$ has semimajor axis $a$, eccentricity $e$, inclination $I$, true anomaly $f$, argument of pericentre $\omega$, and longitude of ascending node $\Omega$. After normalizing all distances to $a'$, the disturbing function reads:

$$R = G m' a'^{-1} (\Delta^{-1} - r \cos \psi) \equiv G m' a'^{-1} \bar{R}, \quad (1)$$

where $r = \alpha(1 - e^2)/(1 + e \cos f)$ is the orbital radius of the particle and $\alpha = a/a'$ is the particle’s normalized semimajor axis that can be larger or smaller than unity, $\psi$ is the angle between the radius vectors of the planet and the test particle, $\Delta^2 = 1 + r^2 - 2r \cos \psi$ is the planet–particle relative distance and

$$\cos \psi = \cos(\Omega - \lambda') \cos(f + \omega) - \sin(\Omega - \lambda') \sin(f + \omega) \cos I \quad (2)$$

The first term of $\bar{R}$ is the gravitational force from the mass $m'$ on the test particle also known as the direct perturbation that we

shall denote $\bar{R}_d$. The second term, that we denote $\bar{R}_i$, is the indirect perturbation that comes from the reflex motion of the star under the influence of the mass $m'$ as the standard coordinate system is chosen to be centred on the star. In the following, we use the notations and steps of the literal expansions for nearly coplanar prograde orbits (Murray & Dermott 1999) and nearly coplanar retrograde orbits by Morais & Namouni (2013a) so that the reader is able to see the similarities and differences of the three expansions.

The classical series of the disturbing function is expanded in powers of $e$ and $\sin^2(I/2)$, which is adequate for nearly coplanar prograde motion since $\sin^2(I/2)$ vanishes for $I = 0$. The classical series can also be used for nearly coplanar retrograde orbits after having applied the procedure devised by Morais & Namouni (2013a) that allows one to get retrograde resonant terms from prograde ones. In essence, the retrograde series is an expansion in terms of $e$ and $\cos^2(I/2)$ where the latter vanishes for $I = 180^\circ$. However, neither the prograde series nor its retrograde counterpart can be used for polar orbits. Instead, inspection of the expressions of $\cos \psi$ (2) and $\Delta$ reveals that a polar expansion has to be done with respect to $e$ and $\cos I$ that vanishes for $I = 90^\circ$. We therefore write

$$\Delta^2 = 1 + r^2 - 2r \cos(\Omega - \lambda') \cos(f + \omega) - 2r \psi, \quad (3)$$

where $\Psi$ is

$$\Psi = - \sin(\Omega - \lambda') \sin(f + \omega) \cos I. \quad (4)$$

Expanding the direct perturbation term $\Delta^{-1}$ in the vicinity of $\Psi = 0$, we have

$$\Delta^{-1} = \sum_{i=0}^{\infty} \frac{(2i)!}{2^i (i!)^2} (r\Psi)^i \Delta_0^{-(2i+1)}, \quad (5)$$

where $\Delta_0^2 = 1 + r^2 - 2r \cos(\Omega - \lambda') \cos(f + \omega)$. Defining $\epsilon = r/\alpha - 1 = \mathcal{O}(e)$ and expanding $\Delta_0^{-(2i+1)}$ around $\epsilon = 0$, we get

$$\Delta_0^{-(2i+1)} = \left( 1 + \sum_{k=1}^{\infty} \frac{\epsilon^k \alpha^k}{k!} \frac{d^k}{d\alpha^k} \right) \rho^{-(2i+1)}, \quad (6)$$

$$\rho^{-(2i+1)} = [1 + \alpha^2 - 2\alpha \cos(\Omega - \lambda') \cos(f + \omega)]^{-(i+1/2)}. \quad (7)$$

The validity of the expansion with respect to zero eccentricity will be discussed in Section 5, where we examine the secular coupling of eccentricity and inclination. In particular, a maximum value of eccentricity will be determined for the polar disturbing function of fourth order.

The next step in the literal expansion is to develop the function $\rho^{-(2i+1)}$ into a two-dimensional Fourier series with respect to the angles $f + \omega$ and $\Omega - \lambda'$ as follows:

$$\rho^{-(2i+1)} = \sum_{\substack{-\infty < j, k < \infty \\ j + k \text{ even}}} \frac{1}{4} b_{i+1/2}^{jk}(\alpha) \cos[k(f + \omega) + j(\Omega - \lambda')], \quad (8)$$

$$b_s^{jk}(\alpha) = \frac{1}{\pi^2} \int_0^{2\pi} \int_0^{2\pi} \frac{\cos(ju + kv) \, du \, dv}{(1 + \alpha^2 - 2\alpha \cos u \cos v)^s}, \quad (9)$$

where $b_s^{jk}(\alpha)$ are two-dimensional Laplace coefficients. The series (8) is summed over even $j + k$ owing to the invariance of the function $\rho^{2i+1}$ (7) with respect to the variable change $(f + \omega + \pi, \Omega - \lambda' + \pi)$ that makes $b_s^{jk} = 0$ if $j + k$ is odd. This invariance can be interpreted as a combined change of reference for $\Omega$ and $\omega$ to the descending node. We remark that the presence of the two-dimensional Laplace coefficients is related to the fact that the two angles $f + \omega$ and $\Omega - \lambda'$ enter the expression of $\rho^{-(2i+1)}$ independently in contrast to the expansions of nearly coplanar orbits, where these angles enter as the sum $f + \omega \pm (\Omega - \lambda')$ where
the ± signs are for prograde and retrograde orbits, respectively (Morais & Namouni 2013a). The two-dimensional Laplace coefficients also satisfy the following relations:

$$b_s^{jk} = b_s^{kj} = b_s^{(-j)k} = b_s^{(-k)j} = b_s^{(-j)(-k)}, \quad (10)$$

$$D b_s^{jk} = \frac{s}{2} \left( b_{s+1}^{(j+1)(k+1)} + b_{s+1}^{(j+1)(k-1)} + b_{s+1}^{(j-1)(k+1)} + b_{s+1}^{(j-1)(k-1)} \right) - 2\alpha s b_{s+1}^{jk}, \quad (11)$$

$$D^n b_s^{jk} = \frac{s}{2} \left( D^{n-1} b_{s+1}^{(j+1)(k+1)} + D^{n-1} b_{s+1}^{(j+1)(k-1)} + D^{n-1} b_{s+1}^{(j-1)(k+1)} + D^{n-1} b_{s+1}^{(j-1)(k-1)} \right) - 2\alpha s D^{n-1} b_{s+1}^{jk} - 2(n-1)s D^{n-2} b_{s+1}^{jk}, \quad (12)$$

where the operator $D = d/d\alpha$. Substituting the series (8) into the expression (6) and the latter into the expansion (5), the direct part of the perturbation is written as the following series:

$$\bar{R}_d = \sum_{\substack{0 \le i, j < \infty \\ -\infty < j, k < \infty \\ j+k \text{ even}}} \frac{1}{i! 2^{i+2} (i!)^2} (r \Psi)^i \epsilon^l A_{i, j, k, l} \cos[k(f + \omega) + j(\Omega - \lambda)], \quad (13)$$

where

$$A_{i, j, k, l} = \alpha^l D^l b_{i+1/2}^{jk}, \quad (14)$$

satisfy the same symmetry of the two-dimensional Laplace coefficients in that $A_{i, j, k, l} = 0$ if $j+k$ is odd. It now remains that we express the terms $\Psi$ and $r$ as functions of the mean longitude and the longitude of pericentre. This can be achieved using the classical elliptic expansions:

$$\begin{align}
\sin f &= 2(1 - e^2)^{\frac{1}{2}} \sum_{s=1}^{\infty} \frac{d}{sde} J_s(se) \sin[s(\lambda - \varpi)], \nonumber \\
\cos f &= -e + 2(1 - e^2)e^{-1} \sum_{s=1}^{\infty} J_s(se) \cos[s(\lambda - \varpi)], \nonumber \\
\frac{r}{\alpha} &= 1 + \frac{e^2}{2} - 2e \sum_{s=1}^{\infty} \frac{d}{s^2de} J_s(se) \cos[s(\lambda - \varpi)]. \tag{15}
\end{align}$$

Upon substituting the expressions (15) into the direct part of the perturbation (13) and truncating it to order $N$ in eccentricity and inclination cosine, we obtain a series that is not quite ready for use for a general $p:q$ resonance. Indeed, the series contains cosine terms whose seemingly unrelated arguments must be transformed to make them represent the same $p:q$ resonance. For example, amongst the various terms that appear to second order in eccentricity and zero order in inclination cosine, there are

$$T_1 = -\frac{e^2}{4} k^2 A_{0,j,k,0} \cos[k\lambda - j\lambda' + (j-k)\Omega], \quad (16)$$

$$T_2 = \frac{e^2}{32} A_{0,j,k,2} \cos[(k-2)\lambda - j\lambda' + (j-k)\Omega + 2\varpi]. \quad (17)$$

Both these terms can be made to correspond to the same resonance $p:q$ by choosing $j=p, k=q$ for $T_1$ and $j=p$ and $k=q+2$ for $T_2$. The transformed terms become

$$T_1 = -\frac{e^2}{4} q^2 A_{0,p,q,0} \cos[q\lambda - p\lambda' + (p-q)\Omega], \quad (18)$$

$$T_2 = \frac{e^2}{32} A_{0,p,q+2,2} \cos[q\lambda - p\lambda' + (p-q-2)\Omega + 2\varpi]. \quad (19)$$

Furthermore for $p \neq 0$ and $q \neq 0$, one needs to account for the resonant terms that are generated by $T_1$ and $T_2$ under the change $p \to -p$ and $q \to -q$ as the series (13) is summed over positive as well as negative $k$ and $j$. This transformation produces two new terms, $T_3$ and $T_4$, that correspond to the same resonance:

$$T_3 = -\frac{e^2}{4} q^2 A_{0,p,q,0} \cos[q\lambda - p\lambda' + (p-q)\Omega] = T_1, \quad (20)$$

$$T_4 = \frac{e^2}{32} A_{0,p,q-2,2} \cos[q\lambda - p\lambda' + (p-q+2)\Omega - 2\varpi]. \quad (21)$$

In the indices of $A_{0,p,q,0}$ and $A_{0,p,q-2,2}$, we use the properties (10) of the two-dimensional Laplace coefficients. The secular terms can be obtained by setting $p=0$ and $q=0$ in $T_1$ and $T_2$ but not in $T_3$ and $T_4$ because the same term would be counted twice. Another way of seeing this is that $(0,0)$ is a fixed point of the transformation $(p \to -p, q \to -q)$.

For the indirect part of the disturbing function, $\bar{R}_i$, one requires only the use of the elliptic expansions (15) to transform true anomalies into mean anomalies and perform the eccentricity expansion. The resulting expressions of the direct and indirect parts of the polar disturbing function are given in the next section. The secular part of the disturbing function is given in Section 5.

# 3 DISTURBING FUNCTION OF NEARLY POLAR ORBITS

## 3.1 Direct part

For a general resonance $p:q$ and an expansion of order $N$ in eccentricity and inclination cosine, the steps described in the previous section show that the direct part of the disturbing function is given as

$$\begin{align}
\bar{R}_d = \sum_{\substack{-N \le k \le N \\ |k| \le m \le N \\ 0 \le n \le N \\ m+n=N}} c_{mn}^k (p, q, \alpha) e^m \cos^n I \cos(\phi - k\omega), \nonumber \\
\phi = q\lambda - p\lambda' + (p-q)\Omega, \tag{22}
\end{align}$$

and $\omega = \varpi - \Omega$ is the argument of pericentre. The force coefficients $c_{mn}^k(p, q, \alpha)$ are given explicitly for the fourth-order series $N=4$ and $k=0, 1, 2, 3$ and $4$ in Tables 1–5. For negative $k$, the force coefficients can be obtained from the identity $c_{mn}^{-k}(p, q, \alpha) = c_{mn}^k(-p, -q, \alpha)$ – see for instance the example of $T_4$ above. Examination of the force coefficients shows the additional relationship $c_{mn}^{-k}(p, q, \alpha) = c_{mn}^k(p, -q, \alpha)$. We note that the force coefficients are applicable to inner and outer perturbers as the semi-major axis ratio can be smaller or larger than unity (Williams 1969).

The force coefficients $c_{mn}^k(p, q, \alpha)$ have a further important property related to the resonance order $o_r = |p-q|$. Examination of $c_{mn}^k(p, q, \alpha)$'s dependence on $A_{i, j, k, l}$ and recalling that $A_{i, j, k, l} = 0$ when $j+k$ is odd (or equivalently $j-k$ is odd<fnref n="1" />) show that for an odd resonance order $o_r$, $c_{mn}^k(p, q, \alpha) = 0$ when $k$ is even. Similarly for an even resonance order $o_r$, $c_{mn}^k(p, q, \alpha) = 0$ when $k$ is odd. This property guarantees that the integer coefficient of the longitude of ascending node, $\Omega$, that reads $p-q+k$ is always even.

<footnote n="1">¹ The integers $j+k$ and $j-k$ have the same parity.</footnote>
<table><caption>**Table 1.** Force coefficients $c_{mn}^0(p, q, \alpha)$ of the term $e^m \cos^n I \cos \phi$.</caption><tbody><tr><td>c<sup>0</sup><sub>00</sub></td><td>$\frac{1}{2}A_{0,p,q,0}$</td></tr><tr><td>c<sup>0</sup><sub>01</sub></td><td>$\frac{\alpha}{8}(A_{1,p-1,q-1,0} - A_{1,p-1,q+1,0} - A_{1,p+1,q-1,0} + A_{1,p+1,q+1,0})$</td></tr><tr><td>c<sup>0</sup><sub>20</sub></td><td>$\frac{1}{8}(-4q^2 A_{0,p,q,0} + 2A_{0,p,q,1} + A_{0,p,q,2})$</td></tr><tr><td>c<sup>0</sup><sub>02</sub></td><td>$\frac{3\alpha^2}{64}(A_{2,p-2,q-2,0} - 2A_{2,p-2,q,0} + A_{2,p-2,q+2,0} - 2A_{2,p,q-2,0} + 4A_{2,p,q,0} - 2A_{2,p,q+2,0} + A_{2,p+2,q-2,0} - 2A_{2,p+2,q,0} + A_{2,p+2,q+2,0})$</td></tr><tr><td>c<sup>0</sup><sub>21</sub></td><td>$-\frac{\alpha}{32}[(4q^2 - 2)(A_{1,p-1,q-1,0} - A_{1,p-1,q+1,0} - A_{1,p+1,q-1,0} + A_{1,p+1,q+1,0}) - 4(A_{1,p-1,q-1,0} + A_{1,p+1,q-1,0}) + A_{1,p-1,q+1,0} + A_{1,p+1,q-1,0} - A_{1,p+1,q+1,0}) - A_{1,p-1,q-1,0} + A_{1,p-1,q+1,0} + A_{1,p+1,q-1,0} + A_{1,p+1,q+1,0}]$</td></tr><tr><td>c<sup>0</sup><sub>03</sub></td><td>$-\frac{5\alpha^3}{256}[3(A_{3,p-3,q+1,0} - A_{3,p-3,q-1,0} - A_{3,p-1,q-3,0} + A_{3,p-1,q+3,0} + A_{3,p+1,q-3,0} + A_{3,p+1,q-3,0} - A_{3,p-3,q+1,0} + A_{3,p-3,q+1,0}) + 9(A_{3,p-1,q-1,0} - A_{3,p-1,q+1,0} - A_{3,p+1,q-1,0} + A_{3,p+1,q+1,0}) - A_{3,p+3,q-3,0} - A_{3,p-3,q+3,0} + A_{3,p-3,q-3,0} + A_{3,p+3,q+3,0}]$</td></tr><tr><td>c<sup>0</sup><sub>40</sub></td><td>$\frac{1}{128}[q^2(16q^2 - 9)A_{0,p,q,0} - 8q^2(A_{0,p,q,1} + A_{0,p,q,2}) + A_{0,p,q,3} + A_{0,p,q,4}]$</td></tr><tr><td>c<sup>0</sup><sub>04</sub></td><td>$\frac{35\alpha^4}{4096}[A_{4,p-4,q+4,0} + A_{4,p-4,q-4,0} + A_{4,p+4,q-4,0} + A_{4,p+4,q+4,0} - 4(A_{4,p-4,q-2,0} + A_{4,p-4,q+2,0} + A_{4,p-2,q-4,0} + A_{4,p-2,q+4,0} + A_{4,p+2,q-4,0} + A_{4,p+2,q+4,0} + A_{4,p+4,q-2,0} + A_{4,p+4,q+2,0}) + 6(A_{4,p-4,q,0} + A_{4,p,q-4,0} + A_{4,p,q+4,0} + A_{4,p+4,q,0}) + 16(A_{4,p-2,q-2,0} + A_{4,p-2,q+2,0} + A_{4,p+2,q+2,0} + A_{4,p+2,q-2,0}) - 24(A_{4,p-2,q,0} - A_{4,p,q-2,0} - A_{4,p,q+2,0} - A_{4,p+2,q,0}) + 36A_{4,p,q,0}]$</td></tr><tr><td>c<sup>0</sup><sub>22</sub></td><td>$-\frac{3\alpha^2}{256}[2(2q^2 - 3)(A_{2,p-2,q-2,0} - 2A_{2,p-2,q,0} + A_{2,p-2,q+2,0} - 2A_{2,p-2,q,0} - 2A_{2,p+2,q-2,0} + A_{2,p+2,q-2,0} + A_{2,p+2,q+2,0}) - 6A_{2,p-2,q-2,1} - A_{2,p-2,q-2,1} + 12A_{2,p-2,q,1} + 2A_{2,p-2,q,2} - 6A_{2,p-2,q+2,1} - A_{2,p-2,q+2,2} + 12A_{2,p,q-2,1} + 2A_{2,p,q-2,2} - 24A_{2,p,q,1} - 4A_{2,p,q,2} + 12A_{2,p,q,2,1} + 2A_{2,p,q+2,q-2,0} - 6A_{2,p+2,q-2,2} + 12A_{2,p+2,q,1} + 2A_{2,p+2,q,2} - 6A_{2,p+2,q+2,1} - A_{2,p+2,q+2,2}]$</td></tr></tbody></table>

To illustrate this property, we shall consider the even order inner 5:1 resonance and the odd order outer 2:9 resonance and write down the corresponding series to second order in eccentricity and inclination cosine. Using Tables 1 and 3 for the 5:1 resonance, we get

$$ \bar{R}_d^{5:1} = [c_{00}^0(5, 1, \alpha) + c_{01}^0(5, 1, \alpha) \cos I + c_{20}^0(5, 1, \alpha)e^2 + \\ + c_{02}^0(5, 1, \alpha) \cos^2 I] \cos(\lambda - 5\lambda' + 4\Omega) \\ + c_{20}^0(5, 1, \alpha)e^2 \cos(\lambda - 5\lambda' + 6\Omega - 2\sigma) \\ + c_{20}^0(5, -1, \alpha)e^2 \cos(\lambda - 5\lambda' + 2\Omega + 2\sigma). \quad (23) $$

The resonant terms $\cos(\lambda - 5\lambda' + 5\Omega - \omega)$ and $\cos(5\lambda - \lambda' + 3\Omega + \omega)$ whose force amplitudes are proportional to $e$ and to cannot appear as the corresponding force coefficients all vanish because $o_r = 4$ is even. They are written as

$$ c_{10}^1(5, 1, \alpha) = -\frac{1}{4}A_{0,5,0,0} = 0, \quad (24) $$

$$ c_{10}^1(5, -1, \alpha) = -\frac{1}{4}(4A_{0,5,2,0} + A_{0,5,2,1}) = 0, \quad (25) $$

<table><caption>**Table 2.** Force coefficients $c_{mn}^1(p, q, \alpha)$ of the term $e^m \cos^n I \cos(\phi - \omega)$.</caption><tbody><tr><td>c<sup>1</sup><sub>10</sub></td><td>$-\frac{1}{4}[2(1-q)A_{0,p,q-1,0} + A_{0,p,q-1,1}]$</td></tr><tr><td>c<sup>1</sup><sub>11</sub></td><td>$\frac{\alpha}{16}[(2q-3)(A_{1,p-1,q-2,0} - A_{1,p-1,q,0} + A_{1,p+1,q,0} - A_{1,p+1,q-2,0}) - A_{1,p-1,q-2,1} + A_{1,p-1,q,1} + A_{1,p+1,q-2,1} - A_{1,p+1,q,1}]$</td></tr><tr><td>c<sup>1</sup><sub>12</sub></td><td>$-\frac{3\alpha^2}{128}[2(2-q)(A_{2,p-2,q-3,0} - 2A_{2,p-2,q-1,0} + A_{2,p-2,q+1,0} - 2A_{2,p,q-3,0} + 4A_{2,p,q-1,0} + A_{2,p+2,q-3,0} - 2A_{2,p+2,q-1,0} + A_{2,p+2,q+1,0}) + A_{2,p-2,q-3,1} - 2A_{2,p-2,q-1,1} + A_{2,p-2,q+1,1} - 2A_{2,p,q-3,1} + 4A_{2,p,q-1,1} - 2A_{2,p,q+1,1} + A_{2,p+2,q-3,1} - 2A_{2,p+2,q-1,1} + A_{2,p+2,q+1,1}]$</td></tr><tr><td>c<sup>1</sup><sub>30</sub></td><td>$\frac{1}{32}[2q(7q - 4q^2 - 3)A_{0,p,q-1,0} + q(4q-1)A_{0,p,q-1,1} - 4A_{0,p,q-1,2} + 2qA_{0,p,q-1,2} - A_{0,p,q-1,3}]$</td></tr><tr><td>c<sup>1</sup><sub>13</sub></td><td>$\frac{5\alpha^3}{512}[(2q-5)(A_{3,p-3,q-4,0} - 3A_{3,p-3,q-2,0} + 3A_{3,p-3,q,0} - A_{3,p-3,q+2,0} - 3A_{3,p-1,q-4,0} + 9A_{3,p-1,q-2,0} - 9A_{3,p-1,q,0} + 3A_{3,p-1,q+2,0} + 3A_{3,p+1,q-4,0} - 9A_{3,p+1,q-2,0} + 9A_{3,p+1,q,0} - 3A_{3,p+1,q+2,0} - A_{3,p+3,q-4,0} + 3A_{3,p+3,q-2,0} - 3A_{3,p+3,q,0} + A_{3,p+3,q+2,0}) - A_{3,p-3,q-4,1} + A_{3,p+3,q-4,1} - A_{3,p+3,q+2,1} + A_{3,p-3,q+2,1} + 3(A_{3,p-3,q-2,1} - A_{3,p-3,q,1} + A_{3,p-1,q-4,1} - A_{3,p-1,q+2,1} - A_{3,p-1,q+2,1} + A_{3,p+3,q,1}) + 9(A_{3,p-1,q,1} - A_{3,p-1,q-2,1} + A_{3,p+1,q-2,1} - A_{3,p+1,q-2,1} + A_{3,p+1,q-2,1} - A_{3,p+1,q,1})]$</td></tr><tr><td>c<sup>1</sup><sub>31</sub></td><td>$-\frac{\alpha}{128}[q(7-18q+8q^2)(A_{1,p-1,q-2,0} - A_{1,p-1,q,0} - A_{1,p+1,q-2,0} + A_{1,p+1,q,0}) + (8-3q-4q^2)(A_{1,p-1,q-2,1} - A_{1,p-1,q,1} - A_{1,p+1,q-2,1} + A_{1,p+1,q,1}) + (7-2q)(A_{1,p-1,q-2,2} - A_{1,p-1,q,2} - A_{1,p+1,q-2,2} + A_{1,p+1,q,2}) + A_{1,p-1,q-2,3} - A_{1,p-1,q,3} - A_{1,p+1,q-2,3} + A_{1,p+1,q,3}]$</td></tr></tbody></table>

<table><caption>**Table 3.** Force coefficients $c_{mn}^2(p, q, \alpha)$ of the term $e^m \cos^n I \cos(\phi - 2\omega)$.</caption><tbody><tr><td>c<sup>2</sup><sub>20</sub></td><td>$\frac{1}{16}[(6-11q+4q^2)A_{0,p,q-2,0} + (6-4q)A_{0,p,q-2,1} + A_{0,p,q-2,2}]$</td></tr><tr><td>c<sup>2</sup><sub>21</sub></td><td>$\frac{\alpha}{64}[(4q^2 - 15q + 12)(A_{1,p-1,q-3,0} - A_{1,p-1,q-1,0} - A_{1,p+1,q-3,0} + A_{1,p+1,q-1,0}) - 4(q-2)(A_{1,p-1,q-3,1} + A_{1,p-1,q-1,1} + A_{1,p+1,q-3,1} - A_{1,p+1,q-1,1}) + A_{1,p-1,q-3,2} - A_{1,p-1,q-1,2} - A_{1,p+1,q-3,2} + A_{1,p+1,q-1,2}]$</td></tr><tr><td>c<sup>2</sup><sub>40</sub></td><td>$\frac{1}{192}[(12+26q-88q^2+68q^3 - 16q^4)A_{0,p,q-2,0} - 2(6-23q+24q^2 - 8q^3)A_{0,p,q-2,1} + (6-9q)A_{0,p,q-2,2} + 4(2-q)A_{0,p,q-2,3} + A_{0,p,q-2,4}]$</td></tr><tr><td>c<sup>2</sup><sub>22</sub></td><td>$\frac{3\alpha^2}{512}[(20-19q+4q^2)(A_{2,p-2,q-4,0} - 2A_{2,p-2,q-2,0} + A_{2,p-2,q,0} - 2A_{2,p,q-4,0} + 4A_{2,p,q-2,0} - 2A_{2,p,q,0} + A_{2,p+2,q-4,0} - 2A_{2,p+2,q-2,0} + A_{2,p+2,q,0}) + (10-4q)(A_{2,p-2,q-4,1} - 2A_{2,p-2,q-2,1} + A_{2,p-2,q,1} - 2A_{2,p,q-4,1} + 4A_{2,p,q-2,1} - 2A_{2,p,q,1} + A_{2,p+2,q-4,1} - 2A_{2,p+2,q-2,1} + A_{2,p+2,q,1}) + A_{2,p-2,q-4,2} - 2A_{2,p-2,q-2,2} + A_{2,p-2,q,2} - 2A_{2,p,q-4,2} + 4A_{2,p,q-2,2} - 2A_{2,p,q,2} + A_{2,p+2,q-4,2} - 2A_{2,p+2,q-2,2} + A_{2,p+2,q,2}]$</td></tr></tbody></table>

$$ c_{11}^1(5, 1, \alpha) = \frac{\alpha}{16}[-(A_{1,4,1,0} - A_{1,4,1,0} + A_{1,6,1,0} - A_{1,6,1,0}) - A_{1,4,1,1} + A_{1,4,1,1} + A_{1,6,1,1} - A_{1,6,1,1}] = 0, \quad (26) $$
<table><caption>**Table 4.** Force coefficients $c_{mn}^3(p, q, \alpha)$ of the term $e^m \cos^n I \cos(\phi - 3\omega)$.</caption><tr><td>c<sup>3</sup><sub>30</sub></td><td><sup>1</sup>&frasl;<sub>96</sub>[(8q<sup>3</sup> - 42q<sup>2</sup> + 62q - 24)A<sub>0,p,q-3,0</sub> - 3(12 - 15q + 4q<sup>2</sup>)A<sub>0,p,q-3,1</sub> + (6q - 12)A<sub>0,p,q-3,2</sub> - A<sub>0,p,q-3,3</sub>]</td></tr><tr><td>c<sup>3</sup><sub>31</sub></td><td><sup>α</sup>&frasl;<sub>384</sub>[(60 - 107q + 54q<sup>2</sup> - 8q<sup>3</sup>)(A<sub>1,p-1,q-2,0</sub> - A<sub>1,p-1,q-4,0</sub> + A<sub>1,p+1,q-4,0</sub> - A<sub>1,p+1,q-2,0</sub>) - 3(20 - 19q + 4q<sup>2</sup>)(A<sub>1,p-1,q-4,1</sub> - A<sub>1,p-1,q-2,1</sub> - A<sub>1,p+1,q-4,1</sub> + A<sub>1,p+1,q-2,1</sub>) + (15 - 6q)(A<sub>1,p-1,q-2,2</sub> - A<sub>1,p-1,q-4,2</sub> + A<sub>1,p+1,q-4,2</sub> - A<sub>1,p+1,q-2,2</sub>) - A<sub>1,p-1,q-4,3</sub> + A<sub>1,p-1,q-2,3</sub> + A<sub>1,p+1,q-4,3</sub> - A<sub>1,p+1,q-2,3</sub>]</td></tr></table>

<table><caption>**Table 5.** Force coefficients $c_{mn}^4(p, q, \alpha)$ of the term $e^m \cos^n I \cos(\phi - 4\omega)$.</caption><tr><td>c<sup>4</sup><sub>40</sub></td><td><sup>1</sup>&frasl;<sub>768</sub>[(120 - 394q + 379q<sup>2</sup> - 136q<sup>3</sup> + 16q<sup>4</sup>)A<sub>0,p,q-4,0</sub> - 4(-60 + 107q - 54q<sup>2</sup> + 8q<sup>3</sup>)A<sub>0,p,q-4,1</sub> + (120 - 114q + 24q<sup>2</sup>)A<sub>0,p,q-4,2</sub> + (20 - 8q)A<sub>0,p,q-4,3</sub> + A<sub>0,p,q-4,4</sub>]</td></tr></table>

$$c_{11}^{1}(5, -1, \alpha) = \frac{\alpha}{16} \left[ \begin{aligned} & -5(A_{1,4,3,0} - A_{1,4,1,0} + A_{1,6,1,0} - A_{1,6,3,0}) - A_{1,4,3,1} + A_{1,4,1,1} \\ & + (20 - 114q + 24q^2)A_{0,p,q-4,2} \\ & + (20 - 8q)A_{0,p,q-4,3} + A_{0,p,q-4,4} \end{aligned} \right] = 0, \quad (27)$$

where we have used the two-dimensional Laplace coefficient relations (10) and the property $A_{r,j,k,l} = 0$ when $j+k$ is odd. We remark that unlike the classical disturbing function for nearly coplanar orbits, the resonance order does not appear in the powers of eccentricity and inclination cosine of the force amplitudes. Moreover, such terms as (23) would not exist in the classical disturbing function as the lowest order pure eccentricity term would be proportional to $e^4$. To dispel doubt on the existence of an unknown symmetry that would make the force coefficients of (23) vanish, we list their non-zero numerical values for nominal resonance $\alpha = 0.341995: c_{00}^0(5, 1) = 0.00069676$, $c_{01}^0(5, 1) = 0.000586162$, $c_{20}^0(5, 1) = 0.00432225$, $c_{02}^0(5, 1) = -0.00111703$, $c_{20}^0(5, 1) = 0.00242241$ and $c_{20}^2(5, -1) = 0.00397014$. Furthermore, in Section 6 we show examples of capture in high-order resonances for low values of the integer $k$ that are smaller than the resonance order $o_r$ (see also the next paragraph). This more general fundamental difference between the two disturbing functions of nearly coplanar and nearly polar orbits will be discussed in Section 4.

The next example is given by the second order in eccentricity and inclination cosine series of the 2:9 resonance that is free of second-order eccentricity terms because $o_r = 7$ is odd. The corresponding expressions are obtained from Table 2 and read:

$$ \bar{R}_d^{2.9} = e[c_{10}^1(2, 9, \alpha) + c_{11}^1(2, 9, \alpha) \cos I] \\ \times \cos(9\lambda - 2\lambda' - 6\Omega - \omega) \\ + e[c_{10}^1(2, -9, \alpha) + c_{11}^1(2, -9, \alpha) \cos I] \\ \times \cos(9\lambda - 2\lambda' - 8\Omega + \omega). \quad (28) $$

The values of the various force coefficients evaluated at nominal resonance, $\alpha = 2.72568$, are $c_{10}^1(2, 9) = 0.000227527$, $c_{11}^1(2, 9) = 0.000149595$, $c_{10}^0(2, -9) = 8.68079 \times 10^{-6}$ and $c_{11}^0(2, -9) = 15.6653 \times 10^{-6}$. Similarly to the previous example, the terms in (28) would not exist in the classical disturbing function as the lowest order pure eccentricity term would be proportional to $e^7$. Using the numerical integration of the full equations of motion, examples of the 2:9 resonance are given in Section 6 where libration is shown

<table><caption>**Table 6.** Force amplitudes and cosine arguments of the indirect part.</caption><thead><tr><th>Cosine argument</th><th>Force amplitude</th></tr></thead><tbody><tr><td>λ' - ω</td><td><sup>3α</sup>&frasl;<sub>4</sub>(1 + cos I)e</td></tr><tr><td>λ - λ'</td><td>-<sup>α</sup>&frasl;<sub>2</sub>(1 + cos I - <sup>1</sup>&frasl;<sub>2</sub>e<sup>2</sup> - <sup>1</sup>&frasl;<sub>2</sub>e<sup>2</sup>cos I - <sup>1</sup>&frasl;<sub>64</sub>e<sup>4</sup>)</td></tr><tr><td>λ - λ' - 2ω + 2Ω</td><td>-<sup>α</sup>&frasl;<sub>48</sub>e<sup>2</sup>(3 + e<sup>2</sup> - 3cos I)</td></tr><tr><td>2λ - λ' - ω</td><td><sup>α</sup>&frasl;<sub>16</sub>(1 + cos I)(3e<sup>3</sup> - 4e)</td></tr><tr><td>2λ - λ' - 3ω + 2Ω</td><td><sup>α</sup>&frasl;<sub>48</sub>(cos I - 1)e<sup>3</sup></td></tr><tr><td>3λ - λ' - 2ω</td><td><sup>3α</sup>&frasl;<sub>16</sub>e<sup>2</sup>(e<sup>2</sup> - cos I - 1)</td></tr><tr><td>3λ - λ' - 4ω + 2Ω</td><td>-<sup>3α</sup>&frasl;<sub>256</sub>e<sup>4</sup></td></tr><tr><td>4λ - λ' - 3ω</td><td>-<sup>α</sup>&frasl;<sub>6</sub>(1 + cos I)e<sup>3</sup></td></tr><tr><td>5λ - λ' - 4ω</td><td>-<sup>125α</sup>&frasl;<sub>768</sub>e<sup>4</sup></td></tr><tr><td>λ' + ω - 2Ω</td><td>-<sup>3α</sup>&frasl;<sub>4</sub>(cos I - 1)e</td></tr><tr><td>λ + λ' - 2Ω</td><td>-<sup>α</sup>&frasl;<sub>128</sub>(64 - 32e<sup>2</sup> - e<sup>4</sup> + 32(e<sup>2</sup> - 2)cos I)</td></tr><tr><td>λ + λ' - 2ω</td><td>-<sup>α</sup>&frasl;<sub>48</sub>e<sup>2</sup>(3 + e<sup>2</sup> + 3cos I)</td></tr><tr><td>2λ + λ' - 3ω</td><td>-<sup>α</sup>&frasl;<sub>48</sub>(1 + cos I)e<sup>3</sup></td></tr><tr><td>2λ + λ' - ω - 2Ω</td><td>-<sup>α</sup>&frasl;<sub>16</sub>(cos I - 1)e(3e<sup>2</sup> - 4)</td></tr><tr><td>3λ + λ' - 4ω</td><td>-<sup>3α</sup>&frasl;<sub>256</sub>e<sup>4</sup></td></tr><tr><td>3λ + λ' - 2ω - 2Ω</td><td>-<sup>3α</sup>&frasl;<sub>16</sub>e<sup>2</sup>(1 - e<sup>2</sup> - cos I)</td></tr><tr><td>4λ + λ' - 3ω - 2Ω</td><td><sup>α</sup>&frasl;<sub>6</sub>(cos I - 1)e<sup>3</sup></td></tr><tr><td>5λ + λ' - 4ω - 2Ω</td><td>-<sup>125α</sup>&frasl;<sub>768</sub>e<sup>4</sup></td></tr></tbody></table>

to occur in the e-resonant term ($k=1$) but also in higher order terms such as $e^3$ ($k=3$) and $e^5$ ($k=5$) that like (28) would not exist in the classical disturbing function of nearly coplanar orbits for the seventh order 2:9 resonance.

## 3.2 Indirect part

The arguments and force amplitudes up to and including fourth order of the disturbing function's indirect part for nearly polar orbits are given in Table 6. The terms present in the expansion concern only resonances of the type $1:q$ with $0 \le q \le 5$. These terms therefore concern only perturbers located inside the object's orbit.

## 4 COMPARING THE DISTURBING FUNCTIONS OF NEARLY COPLANAR ORBITS AND NEARLY POLAR ORBITS

The first main difference between the disturbing functions of nearly coplanar orbits and nearly polar orbits is the fact that the expansion order is not related to the resonance orders. To understand this difference, recall that a literal expansion of the classical disturbing function (of nearly coplanar orbits) to order $N$ in eccentricity and inclination produces cosine terms that represent at most resonances of order $N$. For instance, the fourth order series of Murray & Dermott (1999) applied to a particle perturbed by a planet on a nearly coplanar prograde circular orbit produces only the cosine terms: $\cos[j(\lambda - \lambda') + f_0(\omega, \Omega)]$, $\cos[j\lambda - (j-1)\lambda' + f_1(\omega, \Omega)]$, $\cos[j\lambda - (j-2)\lambda' + f_2(\omega, \Omega)]$, $\cos[j\lambda - (j-3)\lambda' + f_3(\omega, \Omega)]$ and $\cos[j\lambda - (j-4)\lambda' + f_4(\omega, \Omega)]$, where the functions $f_i$ represent the correct combinations of the longitudes of pericentre and ascending node that we do not reproduce explicitly to avoid cumbersome notation. This shows that the possible resonances are of order zero to four but no more. The literal expansion of the disturbing function of nearly polar orbits produces cosine terms for any type of resonance $p:q$ with no restriction on the resonance order $o_r = |p - q|$. A close inspection of the expansion shows that this property is related to the presence of the two independent angles
$f + \omega$ and $\Omega - \lambda'$ that require the use of two-dimensional Laplace coefficients unlike the classical disturbing function that makes use of one-dimensional Laplace coefficients. Therefore, whereas the nearly polar disturbing function of order 4 can be used for the 2:9 resonance discussed in the previous section to study the e-terms associated with $k = 1$ and $k = 3$, a literal expansion of order 7 of the classical disturbing function is required to get the first possible resonant term. This property motivated Ellis & Murray (2000) to come up with an algorithm that produces the force amplitude of a given cosine term of any resonance order without having to expand the classical disturbing function literally.

The second main difference between the two disturbing functions is the fact that the powers of eccentricity and inclination cosine in the force amplitudes of the polar disturbing function are independent of the value of the resonance order $o_r$ (except its parity discussed in the next paragraph). In the classical disturbing function, the lowest order eccentricity and inclination power of the force amplitude of a given cosine term is $o_r$. Indeed for any resonance of order $o_r$, the force amplitude of the cosine term is proportional to $e^{o_r-2k} \sin(I/2)^{2k}$ to the lowest order in $e$ and $I$ where the integer $k$ satisfies $0 \le 2k \le o_r$. The examples of the 5:1 and 2:9 resonance in the previous section showed that their force amplitudes to lowest order in $e$ and $\cos I$ were affine with respect to $\cos I$ as $c_{00}^0 + c_{01}^0 \cos I$ and quadratic in $e$ through the terms, $c_{20}^0 e^2$, $c_{20}^2 e^2$ and $c_{20}^{-2} e^2$ for 5:1 (equation 23). To lowest order in $e$ and $\cos I$, the 2:9 force amplitude is linear in $e$ through the terms $c_{10}^1 e$ and $c_{10}^{-1} e$ (equation 28). In the classical disturbing function, if we seek a linear dependence in eccentricity for 2:9, we must carry along the inclination to the sixth power. The only inclination-free force amplitude is proportional to $e^7$.

When the algorithm of Morais & Namouni (2013a) is applied to the classical disturbing function (of prograde orbits) to produce a series for nearly coplanar retrograde orbits, the force amplitude of a p:q retrograde resonance to lowest order in $e$ and $\cos(I/2)$ is proportional to $e^{(p+q)-2k}\cos(I/2)^{2k}$, where $0 \le 2k \le |p+q|$. This gives even to a first-order resonance (i.e. $|p-q|=1$), the dynamical structure of a high-order resonance. For instance, the planar 1:2 resonance is equivalent to the third-order 1:4 resonance (Morais & Giuppone 2012). Therefore, unlike the polar case, retrograde force amplitudes involve a retrograde resonance order defined as $\bar{o}_r = |p+q|$.

The third main difference between the classical disturbing function and that of nearly polar orbits is the dependence on the parity of the resonance order and the corresponding universal binarity of the force amplitudes of resonant terms. To lowest order in $e$ and $\cos I$, all resonances p:q with an even $o_r = |p-q|$ have force amplitudes that are quadratic with respect to eccentricity and constant with respect to inclination whereas all resonances p:q with an odd $o_r = |p-q|$ have force amplitudes that are linear with respect to eccentricity. As was mentioned in the previous section, this curious behaviour stems from the presence of the two independent angles $f + \omega$ and $\Omega - \lambda'$ in the relative distance $\Delta$. We suspect, but cannot prove, that the property is related to the fact that for a given resonance, one must have prograde as well as retrograde arguments in the same polar series unlike the classical disturbing functions of nearly coplanar prograde or retrograde orbits.

We gain further insight into the structure of the disturbing function for nearly polar orbits by seeking the natural variables with which polar motion can be studied. To do this we recall that instead of using the standard orbital elements, Poincaré devised canonical action-angle variables for studying the three-body problem that have the property of including two actions related to eccentricity and inclination that vanish when motion is exactly circular,

coplanar and prograde. The Poincaré canonical variables are given as the three pairs:

$$ 
\begin{aligned}
\Lambda &= mna^2, \quad \lambda = M + \omega + \Omega, \\
\Gamma &= mna^2 [1 - (1 - e^2)^{1/2}], \quad \gamma = -\omega - \Omega, \\
Z &= mna^2 (1 - e^2)^{1/2} (1 - \cos I), \quad z = -\Omega,
\end{aligned} 
\quad (29) 
$$

where $m$ is the particle's mass, $n$ its mean motion and $M$ its mean longitude. The appropriate variables for polar motion must have the property that the actions related to eccentricity and inclination vanish when motion is exactly polar and circular. The Poincaré action $\Gamma$ satisfies this condition but not $Z$. The latter can be replaced by $Z^* = \Lambda - \Gamma - Z = mna^2(1 - e^2)^{1/2} \cos I$, the normal component of angular momentum. The remaining variables are obtained from the following generating function:

$$ F = (\lambda + z)\Lambda^* + (\gamma - z)\Gamma^* - zZ^*. \quad (30) $$

Using $Z = \partial_z F$, $z^* = \partial_z^* F$ etc., we find

$$ 
\begin{aligned}
\Lambda^* &= \Lambda, \quad \lambda^* = \lambda - \Omega = M + \omega, \\
\Gamma^* &= \Gamma, \quad \gamma^* = -\omega, \\
Z^* &= mna^2 (1 - e^2)^{1/2} \cos I, \quad z^* = \Omega.
\end{aligned} 
\quad (31) 
$$

It can be seen that the choice of the correct variable $Z^*$ modifies the mean longitude $\lambda^*$, the longitude of pericentre $\gamma^*$ and the angle $z^*$ associated with the longitude of ascending node showing that the argument of pericentre is one of the natural angles that should describe polar motion. For comparison, when we modified the Poincaré canonical variables in Namouni & Morais (2015) to study retrograde resonances by choosing the new inclination action as $Z' = 2(\Lambda - \Gamma) - Z$ so that $Z'$ vanishes for exactly coplanar retrograde motion, the mean longitude and longitude of pericentre were modified to $\lambda' = M + \omega - \Omega$ and $\gamma' = \Omega - \omega$, thus producing the natural angles with which retrograde resonances can be studied. We note that our choice of the third canonical action is not unique. For instance instead of $mna^2(1 - e^2)^{1/2} \cos I$, one could employ $mna^2(1 - e^2)^{1/2} (1 - \sin I)$ that has the added advantage of being positive regardless of inclination. The corresponding new angles, however, are no longer function of the old angles; they will also depend on $\Lambda^*, \Gamma^*$ and $Z^*$.

Using the new polar canonical angles, the argument of the cosine terms in the disturbing function (22) is transformed as $\phi_k^{p,q} = \phi - k\omega = p\lambda^* - q(\lambda' - z^*) + k\gamma^*$ implying a simple physical meaning that polar mean longitude need only be measured as if motion were two-dimensional. The particle's longitude of ascending node must be used as a reference line to measure the mean longitude of the planet as the latter two angles lie in the same plane. The remaining term $k\gamma^*$ gives the kth-harmonic that could be excited by the planet. Lastly, we also mention that the new polar canonical variables are related to the classical Delauney variables given as $(L = \Lambda^*, l = \lambda^* + \gamma^*)$, $(G = \Lambda^* - \Gamma^*, g = -\gamma^*)$ and $(H = Z^*, h = Z^*)$. These variables were used by Kozai (1962) to study the secular evolution at large eccentricity and inclination in the three-body problem that will be discussed in the next section.

## 5 SECULAR POTENTIAL AND VALIDITY DOMAIN OF THE DISTURBING FUNCTION

The validity domain of the disturbing function of nearly polar orbits is related to the secular potential that governs the long-term dynamics of the particle. The reason is the large inclination of the particle's orbital plane relative to the planet's that could lead
to large eccentricity and inclination oscillations. In the three-body problem with a planet on a circular orbit, secular evolution of a particle with a non-resonant orbit is given by the Kozai–Lidov potential (Kozai 1962; Lidov 1962). Its integral expression (Quinn, Tremaine & Duncan 1990) is written with our notations as

$$ 
\begin{aligned}
R_{KL} &= \frac{Gm'}{\pi a' \alpha^2 (1 - e^2)^{1/2}} \int_0^{2\pi} R^{-1/2} k K(k) r^2 \, df, \\
k^2 &= \frac{4R}{(R+1)^2 + z^2}, \\
R^2 &= r^2 - z^2, \\
z &= r \sin I \sin(f + \omega),
\end{aligned} 
\quad (32) $$

where $K$ is the complete elliptic integral of the first kind and $r = \alpha(1 - e^2)/(1 + e \cos f)$ is the particle's orbital radius defined in Section 2. The Kozai–Lidov potential is the doubly averaged gravitational potential with respect to the particle's and planet's mean longitudes and does not involve any expansion with respect to eccentricity and inclination. As it depends on the sole angle $\omega$, one can use the Delauney variables and find that both the semimajor axis $\alpha$ and the normal component of angular momentum $(1 - e^2)^{1/2} \cos I$ are constants of secular evolution. Motion generated by the Kozai–Lidov potential thus occurs in the eccentricity-argument of pericentre plane. In writing the expression (32), no assumption was made on the inclination; therefore, it can equally be $\le$ or $\ge$90°. Close examination of (32) and noting the manner in which the normal coordinate, $z$, enters its expression reveal that the secular structures in the $\omega$-plane for prograde and retrograde orbits are identical, and one can study the former and deduce the latter because of the potential's reflection symmetry with respect to the polar plane.

We study the validity domain of the disturbing function by comparing the secular potential it produces with the Kozai–Lidov potential as in essence the former is an expansion of order $N$ of the latter with respect to eccentricity and inclination cosine for nearly polar orbits. This comparison is valid only in the absence of resonance libration but should illustrate the typical values of eccentricity and inclination cosine where the fourth-order series can be used. The literal expansion of Section 2 shows that the secular potential to order $N$ in eccentricity and inclination cosine is given as

$$ \bar{R}_s = \frac{1}{2} b_{\frac{3}{2}}^{00}(\alpha) + \sum_{\substack{0 \le k, n \le N \\ k \le m \le N \\ k, m, n \text{ even} \\ m+n=N}} s_{mn}^k(\alpha) e^m \cos^n I \cos(k\omega), \quad (33) $$

where the corresponding force coefficients $s_{mn}^k(\alpha)$ are given in Table 7 for $N=4$.

We assess the possible large variations of eccentricity and inclination produced by the secular potentials for nearly polar orbits by plotting in the $\omega$-plane the level curves of $R_{KL}a'/Gm'$ and $\bar{R}_s$ for two initially circular orbits located at $\alpha = 2$ and $\alpha = 0.5$, so as to illustrate the effects of internal and external perturbers, respectively. Since these locations are near the 1:3 and 3:1 resonances, respectively, both secular potentials reflect the dynamics of circulating orbits to first order in the perturber's mass.<fnref n="2" /> We will see in Section 6 the evolution of eccentricity and inclination of resonant orbits is somewhat modified by the critical arguments' libration.

<table><caption>Table 7. Force coefficients $s_{mn}^k(\alpha)$ of the secular term $e^m \cos^n I \cos(k\omega)$.</caption><tbody><tr><td rowspan="4">k = 0</td><td>s<sup>0</sup><sub>20</sub></td><td><sup>1</sup>&frasl;<sub>16</sub>(2A<sub>0,0,0,1</sub> + A<sub>0,0,0,2</sub>)</td></tr><tr><td>s<sup>0</sup><sub>02</sub></td><td><sup>3α<sup>2</sup></sup>&frasl;<sub>32</sub>(A<sub>2,0,0,0</sub> - 2A<sub>2,2,0,0</sub> + A<sub>2,2,2,0</sub>)</td></tr><tr><td>s<sup>0</sup><sub>22</sub></td><td><sup>3α<sup>2</sup></sup>&frasl;<sub>128</sub>(6A<sub>2,0,0,0</sub> + 6A<sub>2,0,0,1</sub> + A<sub>2,0,0,2</sub> - 12A<sub>2, 2, 0, 0</sub> - 12A<sub>2, 2, 0, 1</sub> - 2A<sub>2, 2, 0, 2</sub> + 6A<sub>2, 2, 2, 0</sub> + 6A<sub>2, 2, 2, 1</sub> + A<sub>2, 2, 2, 2</sub>)</td></tr><tr><td>s<sup>0</sup><sub>40</sub></td><td><sup>1</sup>&frasl;<sub>256</sub>(4A<sub>0,0,0,3</sub> + A<sub>0,0,0,4</sub>)</td></tr><tr><td rowspan="2">k = 2</td><td>s<sup>0</sup><sub>04</sub></td><td><sup>35α<sup>4</sup></sup>&frasl;<sub>2048</sub>(9A<sub>4,0,0,0</sub> - 24A<sub>4,0,2,0</sub> + 3A<sub>4,0,4,0</sub> + 16A<sub>4, 2, 2, 0</sub> + 3A<sub>4, 4, 0, 0</sub> - 8A<sub>4, 4, 2, 0</sub> + A<sub>4, 4, 4, 0</sub>)</td></tr><tr><td>s<sup>2</sup><sub>20</sub></td><td><sup>1</sup>&frasl;<sub>16</sub>(6A<sub>0,0,2,0</sub> + 6A<sub>0,0,2,1</sub> + A<sub>0,0,2,2</sub>)</td></tr><tr><td rowspan="5">k = 4</td><td>s<sup>2</sup><sub>22</sub></td><td>- <sup>3α<sup>2</sup></sup>&frasl;<sub>256</sub>(20A<sub>2,0,0,0</sub> + 10A<sub>2,0,0,1</sub> + A<sub>2,0,0,2</sub> - 60A<sub>2, 2, 0, 0</sub> - 30A<sub>2, 2, 0, 1</sub> - 3A<sub>2, 2, 0, 2</sub> + 40A<sub>2, 2, 2, 0</sub> + 20A<sub>2, 2, 2, 1</sub> + 2A<sub>2, 2, 2, 2</sub> - 20A<sub>2, 2, 4, 0</sub> - 10A<sub>2, 2, 4, 1</sub> - A<sub>2, 2, 4, 2</sub> + 20A<sub>2, 4, 0, 0</sub> + 10A<sub>2, 4, 0, 1</sub> + A<sub>2, 4, 0, 2</sub>)</td></tr><tr><td>s<sup>2</sup><sub>40</sub></td><td><sup>1</sup>&frasl;<sub>192</sub>(12A<sub>0,0,2,0</sub> - 12A<sub>0,0,2,1</sub> + 6A<sub>0,0,2,2</sub> + 8A<sub>0, 0, 2, 3</sub> + A<sub>0, 0, 2, 4</sub>)</td></tr><tr><td>s<sup>4</sup><sub>40</sub></td><td><sup>1</sup>&frasl;<sub>768</sub>(120A<sub>0,0,4,0</sub> + 240A<sub>0,0,4,1</sub> + 120A<sub>0, 0, 4, 2</sub> + 20A<sub>0, 0, 4, 3</sub> + A<sub>0, 0, 4, 4</sub>)</td></tr></tbody></table>

The initial inclinations are taken as $I(e=0) = 85°, 75°, 65°$ and $55°$ (Fig. 1). Owing to the symmetry with respect to the polar plane, these values produce exactly the same $\omega$-portraits as $I(e=0) = 95°, 105°, 115°$ and $125°$, respectively. The only difference is the inclination range that reads $[I(e=0):180°]$ for retrograde orbits instead of $[I(e=0):0°]$ for prograde orbits. The various level curves in each panel correspond to additional orbits with a normal angular momentum $Z^*$ equal to that of the reference particle namely $(1 - e^2)^{1/2} \cos I = \cos [I(e=0)]$. It is seen on the full Kozai–Lidov potential that particles located inside the planet's orbit ($\alpha = 0.5$) are unstable in the sense that they inevitably reach a near-unit eccentricity corresponding to an orbit that is nearly coplanar with the planet. In the absence of mean motion resonance libration, the time it takes for a particle on a nearly polar orbit with a moderate eccentricity to reach a nearly coplanar orbit is long (see e.g. in Section 6). The fourth-order secular potential is found to reproduce the dynamics quite well for $e \le 0.5$ and $I(e=0) \ge 65°$ in the case of an external perturber (i.e. $\alpha = 0.5$). The potential $\bar{R}_s$ can be used to follow the dynamics on time-scales shorter than the libration around the Kozai–Lidov resonance at $\omega = 90°$. For time-scales comparable to the libration time, one needs to push the expansion order to larger values so as to improve the dynamics' rendition. Particles outside the planet's orbit (i.e. $\alpha = 2$) fall into two types of motion. Up to an eccentricity $e = 0.5$, the argument of pericentre circulates and eccentricity and therefore also inclination have small amplitude variations. The validity domain of the fourth-order polar disturbing function is $e \le 0.5$ for $I(e=0) \ge 75°$ and $e \le 0.2$ for $55° \le I(e=0) \le 65°$. Above $e = 0.5$, eccentricity can be made close to unity by the two Kozai–Lidov resonances at $\omega = 0$ and $90°$ (with the obvious exception of the resonance centres vicinity) and the use of the disturbing function would require a larger (than 4) expansion order like in the case of the external perturber.

In order to understand the secular evolution of the line of nodes that was absent from the previous analysis, we use the secular potential (33) and apply it to the motion of a massless particle perturbed by an internal planet (i.e. $\alpha > 1$) when the particle's orbit is far from the Kozai–Lidov resonances. The corresponding $\omega$-curve is therefore located in the bottom part of the second-row plots

<footnote n="2">² Near-mean motion resonances introduce an additional secular potential effect for circulating orbits whose amplitude is of the second order in the perturber's mass (Hagihara 1972).</footnote>
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Level curves of the secular potentials for an initially circular orbit with semimajor axis ratio $\alpha = 0.5$ (top panels) and $\alpha = 2$ (bottom panels) for four inclination values $I = 85^{\circ}, 75^{\circ}, 65^{\circ}, 55^{\circ}$. The solid blue lines represent the Kozai–Lidov potential (32) and the dashed red lines the fourth-order polar secular potential (33). The thin lines represent the collision singularity.</figcaption></figure>

of Fig. 1. The secular potential (33) restricted to second order in eccentricity and inclination cosine will suffice to follow the motion of the longitude of ascending node. It is written as

$$ \bar{R}_s = s_{20}^0 e^2 + s_{20}^2 e^2 \cos 2\omega + s_{02}^0 \cos^2 I, \quad (34) $$

where we removed the first term of (33) as it does not influence eccentricity and inclination. The Lagrange planetary equations can be written in terms of $e, \omega, I$ and $\Omega$ (Brouwer & Clemence 1961, page 289) and truncated for nearly polar orbits to lowest order in eccentricity and inclination cosine to give

$$ \dot{e} = -(na^2e)^{-1} \partial_{\omega} R_s, \quad \dot{\omega} = (na^2e)^{-1} \partial_e R_s, \quad (35) $$

$$ \dot{I} = -(na^2)^{-1} \partial_{\Omega} R_s, \quad \dot{\Omega} = (na^2)^{-1} \partial_I R_s, \quad (36) $$

where $n = (GM_*a^{-3})^{1/2}$ is the particle's mean motion and $R_s = Gm'a^{-1}\bar{R}_s$ the fully dimensional secular potential. We note how the $I\Omega$-equations (36) differ from those of nearly coplanar orbits in that the variation rates are not inversely proportional to inclination unlike the $e\omega$-equations (35) that keep their classical form. Furthermore, by using the potential (34) in the Lagrange equations, it is found that the inclination $I$ is a constant of motion implying far smaller variations for $I$ than for $e$ when resonant terms are included. The longitude of ascending node's variation rate for nearly polar orbits is also constant, as it depends only on $I$ and given as

$$ \dot{\Omega} = -\frac{n\alpha m' s_{02}^0 \sin 2I}{M_*}. \quad (37) $$

The secular force coefficient that enters the expression of $\dot{\Omega}$ can be approximated by $s_{02}^0 \sim 11(\alpha - 1)(6.5 + 370(\alpha - 1)^{2.1})^{-1}$ for $1 < \alpha \le 5$. Therefore, as $s_{02}^0 > 0$, the line of nodes of prograde nearly polar orbits regresses, whereas that of retrograde nearly polar orbits precesses. The nearer to the polar motion, the smaller the variation rate of the longitude of ascending node. These results are confirmed in the next section.

## 6 EXAMPLES OF POLAR RESONANCE

In this section, we provide an illustration of how the disturbing function helps us to identify the correct resonant arguments associated with the polar motion of a particle that interacts with a Neptune-mass planet on a circular orbit ($m'/M_* = 5.12 \times 10^{-5}$). We will not develop a comparison of the analytical polar disturbing function using the Lagrange equations with numerical integrations, as it is beyond the scope of this work. Instead, we integrate the full equations of motion only to follow the evolution of the particle's orbit and show a variety of polar resonances. We shall consider the following resonances 1:3, 3:1, 2:9 and 7:9.

We learned in Section 3 that the arguments that enter the disturbing function are of the form $\phi_k^{p,q} = \phi - k\omega$, where $\phi = q\lambda - p\lambda' + (p - q)\Omega$. The fundamental mode $k = 0$ occurs only for even-order resonances. It is a pure-inclination term that in principle could librate regardless of eccentricity, as its first-order force amplitude is $c_{00}^0 + c_{01}^0 \cos I$ giving the resonance a pendulum-like dynamical structure almost independent of inclination for polar-like orbits $I \sim 90^{\circ}$. However, nearly polar orbits have a large relative inclination with respect to the planet's orbit and that in turn can force a coupling of eccentricity and inclination variations similar to what was seen in the previous section with the Kozai–Lidov resonance. We therefore illustrate the fundamental mode $k = 0$ by placing the particle directly in the Kozai–Lidov resonance that is coupled to the outer 1:3 mean motion resonance, thus ensuring that the argument of perihelion $\omega$ is stationary and allowing the $k = 0$ mode to librate. Fig. 2 shows as functions of time, normalized to the planet's orbital period $T$, the evolution of the orbital elements along with the resonant angles $\phi_k^{1:3} = 3\lambda - \lambda' - 2\Omega - k\omega$ for $k = 0$ and $-2$. The particle's initial orbital elements are $e = 0.2, \omega = 90^{\circ}, I = 120^{\circ}, \Omega = 0^{\circ}$ and $\alpha = 3^{2/3}$ is the nominal resonance location. The argument $\phi_0^{1:3}$ librates with a variable period starting from a maximum of 4000 $T$ and evolving to a minimum of 1000 $T$. A similar behaviour is forced on the $\phi_{-2}^{1:3}$ argument because of secular resonance. The 1:3 resonance also modifies the $e\omega$-secular structure in that it allows the Kozai–Lidov resonance at $\omega = 90^{\circ}$ to occur at a much lower
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Orbital elements $α, e, ω, I, Ω$ and resonant arguments $φ_0^{1:3}, φ_{-2}^{1:3}$ as a function of time for a particle at the outer 1:3 resonance with a Neptune-mass planet. Initial parameters are eccentricity $e = 0.2$, inclination $I = 120°$, longitude of ascending node $Ω = 0°$, argument of pericentre $ω = 90°$ and relative mean longitude $λ - λ' = 0°$. The bottom right panel is a zoom of $φ_0^{1:3}$ for $0 ≤ t/T' ≤ 10^4$.</figcaption></figure>

eccentricity than that of non-resonant orbits considered in Section 5. However, whereas the eccentricity's variations are moderate, the inclination's are quite small and the line of nodes precesses linearly as the orbit is retrograde confirming the results of Section 5.

On the subject of the effect of the Kozai–Lidov secular potential on polar orbits, we also examine the evolution of a particle that librates in the inner 3:1 resonance located at $α = 3^{-2/3}$ to give an example with an external perturber and ascertain the similarities and differences with the secular evolution of the non-resonant orbits discussed in Section 5. The possible resonant critical arguments now read $φ_k^{3:1} = λ - 3λ' + 2Ω - kω$, where $k$ is an even integer. Placing a particle at the bottom of the $eω$-plane of the first top panel of Fig. 1 with $e = 0.1, ω = 0, Ω = 0$ and choosing an inclination $I = 95°$ whose secular structure away from resonance is identical to that of $I = 85°$ (as explained at the beginning of Section 5), we show in Fig. 3 the evolution of orbital elements as a function of time. As expected from the effects of the secular potential of an external perturber, the argument of pericentre circulates periodically in the narrow strip adjacent to the Kozai–Lidov resonances where the eccentricity reaches a maximum value $e_m = 0.995$. However, the minimum inclination that should have been 180° at maximum eccentricity if the particle were circulating instead of librating in mean motion resonance is now reduced to 150°, indicating that the conserved quantity involving inclination is no longer the normal component of angular momentum like that of the Kozai–Lidov potential of non-resonant orbits. In effect, of the three critical arguments $φ_0^{3:1}, φ_2^{3:1}, φ_4^{3:1}$ the particle is found to stably librate in the $k = 4$ resonance with a 120°-amplitude and $10^5 T'$-period explaining why its secular evolution is not completely described by the Kozai–Lidov potential. It is also interesting to note how other arguments display quasi-librations. For instance, $φ_2^{3:1}$ follows the libration of $φ_4^{3:1}$ over half the libration period only to circulate

rapidly at maximum eccentricity. In a complementary way, the argument $φ_0^{3:1}$ circulates during $φ_4^{3:1}$'s libration and briefly librates at maximum eccentricity.

With the next example, we illustrate how a resonance of order $o_r \gg 1$ can display librations with $|k| < o_r$ and consequently a force amplitude $αe^{|k|}$, a property inherent to the polar disturbing function first encountered in Section 3 and discussed in Section 4. We therefore examine the outer 2:9 resonance example (28) located at $α = (9/2)^{2/3}$ that has an odd resonance order, $o_r = 7$, and possible resonant arguments $φ_k^{2:9} = 9λ - 2λ' - 7Ω - kω$, where $|k| \ge 1$ is an odd integer, and force amplitudes $c_{[k]0}^k(2, 9) e^{|k|}$ to lowest order in eccentricity.

Figs 4–6 show the evolution of the orbital elements as well as the arguments $φ_k^{2:9}$ for $k = 1, 3$ and 5 for three different initial conditions. With the initial parameters $I = 70°, e = 0.05, ω = 0°$ and $Ω = 0°$, the particle librates with the critical argument $φ_1^{2:9}$ while the other arguments circulate (Fig. 4). More precisely, the libration involves two periods: a short one 9800 $T'$ and a longer modulation $2 \times 10^6 T'$. The first period is from the fundamental mode $φ_0^{2:9}$ that circulates on the second time-scale. We know this because the fundamental mode is not related to eccentricity and consequently is not influenced by its secular evolution. The longer period corresponds to the time necessary to counter the fundamental mode's long-term slow circulation with the slow drift of the argument of pericentre $ω$, giving rise to the resonant angle $φ_1^{2:9}$. We note that the resonant arguments $φ_3^{2:9}$ and $φ_5^{2:9}$ display the fast libration of the fundamental mode's short period, yet each is circulating with a slow rate for $φ_3^{2:9}$ and slightly faster for $φ_5^{2:9}$. This shows that when identifying Centaurs and TNOs in polar resonance, one must integrate their orbits over long time spans so as not to misinterpret evolutions such as those described for $φ_3^{2:9}$ and $φ_5^{2:9}$ as true librations in resonance.
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Orbital elements $\alpha, e, \omega, I, \Omega$ and resonant arguments $\phi_0^{3:1}, \phi_1^{3:1}$ and $\phi_4^{3:1}$ as a function of time for a particle at the inner 3:1 resonance with a Neptune-mass planet. Initial parameters are eccentricity $e = 0.1$, inclination $I = 95^\circ$, longitude of ascending node $\Omega = 0^\circ$, argument of pericentre $\omega = 0^\circ$ and relative mean longitude $\lambda - \lambda' = 180^\circ$.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 4.** Orbital elements' and resonant angles' evolution of particle at the outer 2:9 resonance with a Neptune-mass planet. Libration occurs for $\phi_1^{2:9}$. Initial parameters are eccentricity $e = 0.05$, inclination $I = 70^\circ$, longitude of ascending node $\Omega = 0^\circ$, argument of pericentre $\omega = 0^\circ$ and relative mean longitude $\lambda - \lambda' = 90^\circ$.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Orbital elements' and resonant angles' evolution of particle at the outer 2:9 resonance with a Neptune-mass planet. Libration occurs for $\phi_3^{2.9}$. Initial parameters are eccentricity *e* = 0.1, inclination *I* = 70°, longitude of ascending node *Ω* = 0°, argument of pericentre *ω* = 0° and relative mean longitude $\lambda - \lambda' = 90°$.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 6.** Orbital elements' and resonant angles' evolution of particle at the outer 2:9 resonance with a Neptune-mass planet. Libration occurs for $\phi_5^{2.9}$. Initial parameters are eccentricity *e* = 0.2, inclination *I* = 85°, longitude of ascending node *Ω* = 0°, argument of pericentre *ω* = 180° and relative mean longitude $\lambda - \lambda' = 0°$.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Orbital elements and resonant angles' evolution of TNO 471325 at the outer 7:9 resonance with a Neptune-mass planet in the context of the three-body problem. Initial parameters are eccentricity *e* = 0.3, inclination *I* = 110°, longitude of ascending node *Ω* = 0°, argument of pericentre *ω* = 90° and relative mean longitude *λ* − *λ*′ = 180°.</figcaption></figure>

Increasing the particle's eccentricity to *e* = 0.1 has two effects: the argument of perihelion's regression is faster (see below) leading to resonance with the critical argument *φ*<sub>3</sub><sup>2:9</sup> (Fig. 5). The two libration periods are decreased: the short one to 4500 *T* and the long one to 1.66 × 10<sup>6</sup> *T*. The cautionary note that we pointed out in the previous example is still valid for the new eccentricity as the arguments *φ*<sub>1</sub><sup>2:9</sup> and *φ*<sub>5</sub><sup>2:9</sup> circulate but display the librations associated with the short period. The situation is even more misleading for *φ*<sub>5</sub><sup>2:9</sup> because if the integration time span were restricted to the interval [7,10] × 10<sup>5</sup> *T*, the argument would seem to be genuinely librating.

For the example in Fig. 6, we increased the inclination to $I = 85^{\circ}$ and eccentricity to $e = 0.2$, resulting in the libration of the critical argument $\phi_2^{2:9}$ with a small amplitude and the two periods $1400 T'$ (fast one associated with the fundamental mode $k=0$) and $5 \times 10^5 T'$ (the slower one associated with $k=5$). With this eccentricity, the behaviour of the arguments $\phi_1^{2:9}$, $\phi_3^{2:9}$ is less misleading regarding the importance of the short period librations discussed with the previous examples and their evolution is more clearly circulating. The reason is the small amplitude of the short period libration seen in $\phi_5^{2:9}$.

The 2:9 resonance examples also illustrate how the inclination's relative variation is small and how the line of nodes of subpo- lar orbits regresses linearly at a rate that decreases as the incli- nation approaches 90°. We also note that when the eccentricity is increased, the librating critical argument's integer *k* increases. This trend however is not generally valid, as the next example will show.

With the last example we examine the outer 7:9 resonance in which TNO 471325 (Chen et al. 2016) could currently be captured. Since our interest is the polar disturbing function in the context of

the three-body problem, these simulations do not reflect the actual evolution of the object but will give us an idea on the possible resonant arguments that might be involved. With this in mind, the initial orbital elements are: semimajor axis $\alpha = 1.182$, eccentricity is $e = 0.3$ and inclination $I = 110^{\circ}$. We choose $\Omega = 0^{\circ}$ and $\omega = 90^{\circ}$. As an even order resonance, the permissible resonant arguments are $\phi_k^{7:9} = 9\lambda - 7\lambda' - 2\Omega - k\omega$, where *k* is an even integer. Fig. 7 shows the evolution of the orbital elements as well as the arguments $\phi_k^{7:9}$ for $k = 2, 4$ and $6$. It is seen that the argument $\phi_4^{7:9}$ for TNO 471325 in the three-body problem librates around 180° with an amplitude and a period of 68° and 15550 $T'$, respectively. The argument of pericentre $\omega$ circulates rapidly, the longitude of ascending node $\Omega$ precesses linearly and the inclination has moderate variations as predicted in Section 5. The both resonant arguments $\phi_2^{7:9}$ and $\phi_6^{7:9}$ circulate but only the latter displays temporary librations similar to those of the previous resonance.

Increasing the eccentricity to $e = 0.4$ and initializing $\omega$ at 180° makes the orbit librate with the critical argument $\phi_2^{7:9}$ with an amplitude of 120° and a period of 22 000 $T'$ (Fig. 8). Thus, the trend noted in the previous example, about how a larger *k* could be associated with a larger *e*, is not confirmed. The remaining arguments circulate but now it is $\phi_4^{7:9}$ that displays temporary librations whereas $\phi_6^{7:9}$ circulates with the fastest rate.

When the eccentricity is increased further to $e = 0.64$ ($\omega = 180^{\circ}$), the particle librates with the critical argument $\phi_6^{7:9}$ of amplitude 70° and period of 16 000 $T'$ (Fig. 9). The arguments $\phi_2^{7:9}$ and $\phi_4^{7:9}$ circulate rapidly without temporary librations.

We conclude that for the same location and inclination as TNO 471325, the resonant argument depends strongly on the observed eccentricity. It is likely that $\phi_4^{7:9}$ is the correct librating critical
<figure><img src="image_9.png" /><figcaption>**Figure 8.** Orbital elements and resonant angles' evolution of particle at the outer 7:9 resonance with a Neptune-mass planet. Initial parameters are eccentricity *e* = 0.4, inclination *I* = 110°, longitude of ascending node *Ω* = 0°, argument of pericentre *ω* = 180° and relative mean longitude *λ* − *λ*′ = 180°.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 9.** Orbital elements and resonant angles' evolution of particle at the outer 7:9 resonance with a Neptune-mass planet. Initial parameters are: eccentricity *e* = 0.64, inclination *I* = 110°, longitude of ascending node *Ω* = 0°, argument of pericentre *ω* = 180° and relative mean longitude *λ* − *λ*′ = 180°.</figcaption></figure>

argument as the effect of the other planets on resonant polar asteroids must be reduced because of their peculiar orbital geometry unless there is strong interaction between the planets such as a mean motion resonance that carries over to the motion of the polar asteroid.

# 7 CONCLUDING REMARKS

The main technical results of this work consist of (i) an explicit algorithm for generating a literal expansion of the disturbing
function for nearly polar orbits of general order *N* in eccentricity and inclination cosine (ii) the explicit form of the fourth-order polar disturbing function through the direct part (22) with its force coefficients given in Tables 1–5, the indirect part written explicitly in Table 6, and the secular potential (33) whose force coefficients are given in Table 7. Beyond the technical results, our original motivation for deriving a literal expansion of the disturbing function for nearly polar orbits is the realization that general attitude regarding resonance identification for polar Centaurs and TNOs is based on decades, and for some aspects more than a century, of use in planetary dynamics of the classical disturbing function derived for nearly coplanar prograde orbits. It is therefore not surprising that we have revealed new features unseen in the classical disturbing function, especially the structure of the force amplitudes that define resonance strength. In particular, the fact that regardless of resonance order, a particle can librate in the lowest harmonics (small *k* in equation 22) of the disturbing function is interesting, as it explains an important observation that was made in our numerical studies of resonance capture at arbitrary inclination (Namouni & Morais 2015, 2017) that resonance order is not a good indicator of resonance strength nor capture efficiency. This observation was particularly striking for the outer 1:5 resonance that exceeded 80 per cent capture efficiency for the most eccentric nearly polar orbits (fig. 6 of Namouni & Morais 2017). TNO 471325 provides a good example of a near polar asteroid locked in resonance. Our three-body simulations suggest that the resonant critical argument is $\phi_4^{7:9}$. Further simulations including all Solar system planets are required to confirm this possibility and discover yet more Centaurs and TNOs in polar resonance with the giant planets.

## ACKNOWLEDGEMENTS

The authors thank Zoran Knežević for his thorough and helpful review of the manuscript. FN thanks the University of Rio Claro (UNESP) for their hospitality during his stay when part of this work was performed. This work was supported by grant no. 2015/17962-5 of São Paulo Research Foundation (FAPESP).

## REFERENCES

Brouwer D., Clemence G. M., 1961, Celestial Mechanics. Academic Press, New York

Chen Y.-T. et al., 2016, ApJ, 827, L24

Ellis K. M., Murray C. D., 2000, Icarus, 147, 129

Gladman B. et al., 2009, ApJ, 697, L91

Hagihara Y., 1972, Celestial mechanics. Vol. 2, Perturbation Theory. Massachusetts Institute of Technology (MIT), Cambridge, MA

Kozai Y., 1962, AJ, 67, 591

Laskar J., Boué G., 2010, A&A, 522, A60

Lidov M. L., 1962, Planet. Space Sci., 9, 719

Morais M. H. M., Giuppone C. A., 2012, MNRAS, 424, 52

Morais M. H. M., Namouni F., 2013a, Celest. Mech. Dyn. Astron., 117, 405

Morais M. H. M., Namouni F., 2013b, MNRAS, 436, L30

Murray C. D., Dermott S. F., 1999, Solar System Dynamics. Cambridge Univ. Press, Cambridge

Namouni F., Morais M. H. M., 2015, MNRAS, 446, 1998

Namouni F., Morais M. H. M., 2017, MNRAS, 467, 2673

Quinn T., Tremaine S., Duncan M., 1990, ApJ, 355, 667

Williams J. G., 1969, PhD thesis, University of California at Los Angeles

This paper has been typeset from a *TEX/LATEX* file prepared by the author.