# The disturbing function for asteroids with arbitrary inclinations

F. Namouni, M. Morais

► To cite this version:

F. Namouni, M. Morais. The disturbing function for asteroids with arbitrary inclinations. Monthly Notices of the Royal Astronomical Society, 2018, 474 (1), pp.157-176. 10.1093/mnras/stx2636. hal-02398839

HAL Id: hal-02398839

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

Submitted on 22 Nov 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished 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
<figure><img src="image_2.png" /></figure>

# The disturbing function for asteroids with arbitrary inclinations

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

Accepted 2017 October 6. Received 2017 October 6; in original form 2017 July 7

## ABSTRACT

The classical disturbing function of the three-body problem widely used in planetary dynamics studies is an expansion of the gravitational interaction of the three-body problem with respect to zero eccentricity and zero inclination. This restricts its validity to nearly coplanar orbits. Motivated by the dynamical study of asteroids, Centaurs and transneptunian objects with arbitrary inclinations, we derive a series expansion of the gravitational interaction with respect to an arbitrary reference inclination that generalizes our work on the polar and retrograde disturbing functions. The new disturbing function is similar to the polar one and may model any resonance as expansion order is unrelated to resonance order. The powers of eccentricity and inclination of the force amplitude of a *p*:*q* resonance depend only on the parity of the resonance order |*p* − *q*|. Disturbing functions with non-zero reference inclinations are thus physically different from the classical disturbing function as the former are based on the three-dimensional three-body problem and the latter on the two-dimensional one. We illustrate the use of the new disturbing function by showing that what is known as pure eccentricity resonances are intrinsically dependent on inclination contrary to the prediction of the classical disturbing function. We determine the inclination dependence of the resonance widths of the 2:1 and 3:1 prograde and retrograde inner resonances with Jupiter as well as those of the asymmetric librations of the 1:2 and 1:3 prograde outer resonances with Neptune.

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

## 1 INTRODUCTION

The recent discoveries that Centaurs and transneptunian objects (TNOs) with polar and retrograde inclinations may be captured in mean motion resonance by the Solar system's planets (Morais & Namouni 2013b, 2017; Wiegert, Connors & Veillet 2017) as well as the increasing number of TNOs on polar orbits (Gladman et al. 2009; Chen et al. 2016) opened the new field of study on the origins and evolution of asteroids with arbitrary inclinations. Hopes for an analytical understanding of the intricacies of asteroid dynamics and capture at arbitrary inclinations have been hampered by the inadequacy of the classical analytical tools. Chief among them is the most basic one: the disturbing function of the three-body problem, the series expansion in the vicinity of circular and coplanar orbits of the gravitational interaction of two bodies that revolve around the sun. Lacking such a basic tool for arbitrary inclinations, analytical studies could concentrate only on nearly coplanar retrograde orbits (Morais & Namouni 2013a) whereas intensive numerical simulations probed capture mechanisms and dynamical

stability at arbitrary inclination (Namouni & Morais 2015, 2017c; Morais & Namouni 2016). Of the most interesting and unexpected numerical results is the high efficiency of resonance capture for polar and retrograde inclinations. The basic physical origin of such efficiency is known to be the encounter geometry as the retrograde asteroid encounters the planet for a shorter duration and larger velocity than if it were on a prograde orbit making retrograde resonances harder to destabilize (Namouni & Morais 2017a). A more precise and quantitative analytical understanding of how capture operates would be possible if a disturbing function were derived to cover all inclinations and not only nearly coplanar prograde motion (Ellis & Murray 2000) and nearly coplanar retrograde motion (Morais & Namouni 2013a).

A step in the right direction was recently achieved with the disturbing function for polar Centaurs and TNOs (Namouni & Morais 2017b, hereafter Paper I). A series expansion was carried out literally for nearly polar motion in powers of eccentricity and inclination cosine. The resulting disturbing function was quite different from the classical one for nearly circular coplanar prograde motion (Ellis & Murray 2000) and challenged long held beliefs on the dependence of the force amplitude of resonant terms on the asteroids's orbital elements. For instance, it was found that any

<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">²Instituto de Geociências e Ciências Exatas, Universidade Estadual Paulista (UNESP), 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>
$p:q$ mean motion resonance has, to lowest order in eccentricity, a force amplitude proportional to eccentricity $e$ if the resonance order $|p-q|$ is odd and proportional to $e^2$ if the order is even. Such behaviour originates from the fact that the expansion is carried out with respect to a three-dimensional unperturbed configuration of reference inclination $I_r = 90^\circ$ unlike the classical expansion that is built on the two-dimensional three-body problem of reference inclination $I_r = 0^\circ$. This difference adds a new degree of freedom to the disturbing function embodied by new two-dimensional Laplace coefficients that bear the semimajor axis dependence of the force amplitude.<fnref n="1" /> The polar disturbing function helped identify TNO (471325) as the first TNO in polar resonance with Neptune (Morais & Namouni 2017).

The three available disturbing functions: the classical one for nearly circular coplanar prograde motion (Ellis & Murray 2000), its companion, the nearly circular coplanar retrograde expansion (Morais & Namouni 2013a) and the polar expansion seem to cover large portions of inclination space for the study of resonance dynamics. However, this situation is not satisfactory as it poses the problem of in-between inclinations. If an asteroid's orbit has an inclination of 45° or its polar symmetric, 135°, which one of the three expansions must be used? It was shown in Paper I that below an inclination of 55° (and beyond the polar symmetric 125°), the polar disturbing function to fourth order in eccentricity and inclination cosine is imprecise and higher expansion orders are required. The classical disturbing function and its retrograde companion suffer the same limitations at different inclination thresholds (Knežević et al. 1995). Therefore, the range of intermediate inclinations is substantial and may not be modelled by any of the known disturbing functions. In order to produce an analytical theory of mean motion resonance capture and crossing at arbitrary inclination, the opportune action is to develop an expansion with respect to eccentricity and an arbitrary and unspecified reference inclination $I_r$. If this is possible, then even for inclinations such as 10° or 100°, the corresponding disturbing function with an equal reference inclination will be far more accurate than the classical expansion or the polar one even if they can be used, and will require fewer terms in the series to model their dynamics. What is said more precisely in mathematical terms is that approximating a function $f(x)$ by Taylor-expanding it in the vicinity of $x_0 = 0$ in order to apply it to the vicinity of an $x_1 \gg 0$ is neither precise nor acceptable when $f(x)$ can be Taylor-expanded in the vicinity of $x_1$ directly.

In this paper, we expand the gravitational potential of an asteroid with an eccentric and inclined orbit, which interacts with a planet on a circular orbit, around zero eccentricity and a reference orbital plane of unspecified inclination $I_r$. As noted in Paper I, disturbing functions come in two types: the first is that of expansions with respect to eccentricity and some function of inclination. Detailed expositions and historical perspectives on the classical disturbing function for nearly coplanar orbits are found in Brouwer & Clemence (1961). The polar series expansion of Paper I is another such example. The second type is that of expansions with respect to the semimajor axis ratio that were recently reformulated in modern language by Laskar & Boué (2010). The first type is the most widely used in nearly all subjects in planetary dynamics because of a significant shortcoming in the second type. Indeed, series expansions with respect to semimajor axis ratio $\alpha$ are notorious for their slow convergence if the condition $\alpha \ll 1$ is not satisfied, requiring an exceed-

ingly large expansion order (see the discussion of the $\alpha$-expansion of the two-dimensional Laplace coefficients in Appendix A). This feature confines the second type disturbing functions to the study of the perturbation of distant bodies such as the Sun on natural and artificial satellites. For our application needs, namely the dynamics of asteroids, Centaurs and TNOs that usually do not satisfy $\alpha \ll 1$, only series expansions of the first type are useful.

We therefore develop a series expansion in powers of eccentricity $e$ and inclination variable $s = \sin(I - I_r)$. Surprisingly the resulting disturbing function provided explicitly to fourth order in $e$ and $s$ is not particularly more involved than the polar expansion of Paper I itself a special case thereof for $I_r = 90^\circ$. Most of the new properties of the polar disturbing function are present in the general disturbing function and show that the latter two expansions belong to a class of disturbing functions based on the three-dimensional three-body problem unlike the classical disturbing function that is based on a perturbation of the two-dimensional problem through the assumption of nearly coplanar orbits for the asteroid and the planet. In Section 2, we provide the literal expansion of the disturbing function with an arbitrary reference inclination and show how the new two-dimensional Laplace coefficients that are commonly functions of the semimajor axis ratio of the asteroid and the planet, become functions of the reference inclination $I_r$ as well. The properties of the new Laplace coefficients are discussed in Appendix A. The reader who is not interested in the derivation details of the new disturbing function can skip Section 2 and find the expression of the series expansion in Section 3. The properties of the new disturbing function and how it reduces to the classical one whose $I_r = 0^\circ$, its companion expansion for retrograde motion whose $I_r = 180^\circ$ and the polar expansion of Paper I whose $I_r = 90^\circ$ are discussed in Section 4. In Section 5, we illustrate the importance of the new disturbing function by deriving the resonance width of pure eccentricity resonances whose resonant arguments do not include inclination through the absence of the longitude of ascending node in their expressions. According to the classical disturbing function, these resonances are independent of inclination. We show that this conclusion is erroneous and pure eccentricity resonances depend strongly on inclination. In order to perform the resonance width calculations and apply them to prograde and retrograde inner resonances with Jupiter and outer resonances with Neptune, we develop a new pendulum model of resonance based on two resonant harmonics that is able to reproduce the known asymmetric librations of outer resonances (Beaugé 1994; Bruno 1994; Malhotra 1996; Winter & Murray 1997). Section 6 contains our concluding remarks.

## 2 LITERAL EXPANSION WITH ARBITRARY REFERENCE INCLINATION

We perform the literal expansion for arbitrary inclination using the same notations and steps of the literal expansion for nearly polar orbits (Paper I). The notations are also consistent with the development of the literal expansions for nearly coplanar prograde orbits (Murray & Dermott 1999) and nearly coplanar retrograde orbits (Morais & Namouni 2013a).

We consider an asteroid that moves under the gravitational influence of the sun of mass $M_*$ and a planet of mass $m' \ll M_*$. The planet's orbit is circular of radius $a'$ and longitude angle $\lambda'$. The reference plane is defined by the sun–planet orbit. The asteroid's osculating Keplerian orbit with respect to the star has semimajor axis $a$, eccentricity $e$, inclination $I$, true anomaly $f$, argument of pericentre $\omega$ and longitude of ascending node $\Omega$. After normalizing

<footnote n="1">¹ These aspects are discussed in detail in Section 4.</footnote>
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 \psi)$ is the orbital radius of the asteroid and $\alpha = a/a'$ is the asteroid's normalized semimajor axis that may be larger or smaller than unity, $\psi$ is the angle between the radius vectors of the planet and the asteroid, $\Delta^2 = 1 + r^2 - 2r \cos \psi$ is the planet-asteroid relative distance and

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

The first term of $\bar{R}$ is the direct perturbation that we denote $\bar{R}_d$ and the second term, that we denote $\bar{R}_t$, is the indirect perturbation that comes from the reflex motion of the star under the influence of the planet as the standard coordinate system is chosen to be centred on the star.

The classical disturbing function is expanded in powers of $e$ and $\sin^2(I/2)$ with respect to a reference inclination $I_r = 0^\circ$. It is therefore valid for nearly coplanar prograde motion as $\sin^2(I/2)$ vanishes for $I = I_r$. The classical disturbing function may be transformed into an expansion in terms of $e$ and $\cos^2(I/2)$ valid for nearly coplanar retrograde orbits as $\cos^2(I/2)$ vanishes for the reference inclination $I_r = 180^\circ$ (Morais & Namouni 2013a). For nearly polar orbits and hence $I_r = 90^\circ$, the disturbing function may be expanded in powers of $e$ and $\cos I$ (Paper I). For an arbitrary reference inclination $I_r$, we write $I = I_r + \delta I$ in the expression of $\cos \psi$ and reduce it as follows:

$$\cos \psi = \cos(\Omega - \lambda') \cos(f + \omega) - \sin(\Omega - \lambda') \sin(f + \omega) \cos I_r + \Psi \quad (3)$$

$$\Psi = \left( [1 - (1 - s^2)^{1/2}] \cos I_r + s \sin I_r \right) \times \sin(f + \omega) \sin(\Omega - \lambda'), \quad (4)$$

where $s = \sin \delta I = \sin(I - I_r)$. At this point, we note that the interpretation of the new disturbing function is twofold because of the presence of the reference inclination. We can choose and expand the gravitational interaction with respect to $e$ and $s$ and obtain a series that generalizes the classical nearly coplanar and nearly polar disturbing functions. This will be done shortly. The second interpretation arises if we set $s \equiv 0$ at this level of the expansion development or in the full series expansion given in Section 3, then the new disturbing function is only an expansion with respect to eccentricity and keeps a full (unexpanded) dependence on inclination as the reference inclination is actually the asteroid's inclination for $s \equiv 0$. This aspect will be discussed further in Section 4.3. With this in mind, we may write the asteroid-planet relative distance as follows:

$$\Delta^2 = 1 + r^2 - 2r \cos(\Omega - \lambda') \cos(f + \omega) + 2r \sin(\Omega - \lambda') \sin(f + \omega) \cos I_r - 2r \Psi. \quad (5)$$

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

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

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

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

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

The next step 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, I_r) \times \cos[k(f + \omega) + j(\Omega - \lambda')], \quad (9)$$

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

where $b_s^{jk}(\alpha, I_r)$ are two-dimensional Laplace coefficients that reduce to those of the polar expansion for $I_r = 90^\circ$. The index $s$ is a half-integer and should not be confused with the inclination variable. We chose these notation conventions to be in accord with the significant body of literature on the classical disturbing function both historical and contemporary. The series (9) is summed over even $j + k$ owing to the invariance of the function $\rho^{2i+1}$ (8) with respect to the variable change $(f + \omega + \pi, \Omega - \lambda' + \pi)$ that makes $b_s^{jk} = 0$ if $j + k$ is odd. The appearance of the two-dimensional Laplace coefficients is related to the presence of the two independent angles $f + \omega$ and $\Omega - \lambda'$ in the expression of $\rho^{-(2i+1)}$ just like the case of the polar expansion (Paper I) but unlike the case of the expansions of nearly coplanar orbits where those angles enter only as the sum $f + \omega \pm (\Omega - \lambda')$ where the $\pm$ signs are for prograde and retrograde orbits, respectively (Morais & Namouni 2013a). The properties of the two-dimensional Laplace coefficients with an arbitrary reference inclination are discussed in Appendix A.

Substituting the series (9) into the expression (7) and the latter into the expansion (6), the direct part of the perturbation is written as the following series:

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

where

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

satisfy the same restriction on the two-dimensional Laplace coefficients in that $A_{i,j,k,l} = 0$ if $j + k$ is odd, as well as additional properties of the Laplace coefficients discussed in Appendix A. Next, we express $\Psi$ and $r$ in terms of the mean longitude and the longitude of pericentre using the classical elliptic expansions:

$$\sin f = 2(1 - e^2)^{1/2} \sum_{\sigma=1}^{\infty} \frac{d}{s \, d\sigma} J_{\sigma}(\sigma e) \sin[\sigma(\lambda - \sigma)], \\ \cos f = -e + 2(1 - e^2)e^{-1} \sum_{\sigma=1}^{\infty} J_{\sigma}(\sigma e) \cos[\sigma(\lambda - \sigma)], \\ \frac{r}{\alpha} = 1 + \frac{e^2}{2} - 2e \sum_{\sigma=1}^{\infty} \frac{d}{\sigma^2 \, d\sigma} J_{\sigma}(\sigma e) \cos[\sigma(\lambda - \sigma)]. \quad (13)$$
Upon substituting the expressions (13) into the direct part of the perturbation (11) and truncating it to order N in eccentricity e and inclination s, we obtain a series that needs to be transformed further to model a general p:q resonance. This classical reduction step for literal expansions arises because seemingly unrelated arguments actually pertain to the same p:q resonance. For example, among the various terms that appear to second order in s = sin δI and zero order in eccentricity, there are

$$
\begin{align}
T_1 &= \frac{\alpha^2 s^2}{64} A_{2,j,k,0} \cos[k\lambda - j\lambda' + (j-k)\Omega], \nonumber \\
T_2 &= -\frac{3\alpha^2 s^2 \cos(2I_r)}{512} A_{2,j,k,0} \nonumber \\
    &\qquad \times \cos[(k-2)\lambda - (j+2)\lambda' + (j-k+4)\Omega]. \tag{15}
\end{align}
$$

As according to equation (11), the indices j and k can take any positive and negative integer values independently of each other; both these terms may be made to correspond to the same resonance p:q by choosing j = p, k = q for T₁ and j = p - 2 and k = q + 2 for T₂ to find

$$
\begin{align}
T_1 &= \frac{\alpha^2 s^2}{64} A_{2,p,q,0} \cos[q\lambda - p\lambda' + (p-q)\Omega], \tag{16} \\
T_2 &= -\frac{3\alpha^2 s^2 \cos(2I_r)}{512} A_{2,p-2,q+2,0} \nonumber \\
    &\qquad \times \cos[q\lambda - p\lambda' + (p-q)\Omega]. \tag{17}
\end{align}
$$

Furthermore, for $p \neq 0$ and $q \neq 0$, $(-p):( -q)$, cosine terms obviously correspond to the p:q resonance as the series (11) is summed over positive as well as negative k and j as noted previously. Changing p and q to their opposites in T₁ and T₂ produces two new terms, T₃ and T₄, which correspond to the p:q resonance:

$$
\begin{align}
T_3 &= \frac{\alpha^2 s^2}{64} A_{2,p,q,0} \cos[q\lambda - p\lambda' + (p-q)\Omega], \tag{18} \\
T_4 &= -\frac{3\alpha^2 s^2 \cos(2I_r)}{512} A_{2,p+2,q-2,0} \nonumber \\
    &\qquad \times \cos[q\lambda - p\lambda' + (p-q)\Omega]. \tag{19}
\end{align}
$$

In the indices of A<sub>0,p,q,0</sub> and A<sub>0,p+2,q-2,2</sub>, we use the properties (A7) and (A8) of the two-dimensional Laplace coefficients. The secular terms may be obtained by setting p = 0 and q = 0 in T<sub>1</sub> and T<sub>2</sub> but not in T<sub>3</sub> and T<sub>4</sub> because the same term would be counted twice.<sup>2</sup>

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

3 DISTURBING FUNCTION WITH ARBITRARY
REFERENCE INCLINATION

For a *p*:q cosine term and an expansion order *N* in eccentricity *e* and inclination *s*, 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{k \text{ m } \\ -N \le k \le N \\ |k| \le m \le N \\ 0 \le m \le N \\ m+n=N}} c_{mn}^k (p, q, \alpha, I_r) e^m s^n \cos \phi_k^{p,q}, \nonumber \\
\phi_k^{p,q} &= q\lambda - p\lambda' + (p-q)\Omega - k\omega. \tag{20}
\end{align}
$$

The force coefficients $c_{mn}^k(p, q, \alpha, I_r)$ are given explicitly for the fourth-order series $N=4$ and $k=0, 1, 2, 3$ and 4 in Tables B1–B5. For negative $k$, the force coefficients may be obtained from the identity $c_{mn}^{-k}(p, q, \alpha, I_r) = c_{mn}^k(-p, -q, \alpha, I_r + 180^\circ)$ where $I_r + 180^\circ$ is defined modulo $180^\circ$. Examination of the force coefficients shows the additional relationship:

$$
c_{mn}^{-k}(p, q, \alpha, I_r) = c_{mn}^{k}(p, -q, \alpha, I_r + 180^{\circ}). \quad (21)
$$

As no assumption was made regarding the ratio α of the semi-
major axis of the asteroid's orbit and that of the perturber's circular
orbit, the force coefficients are valid for α > 1 as well as α < 1
(Williams 1969).

At this point, we caution the user of the arbitrary inclination disturbing function that the fourth-order expansion provided in Appendix B does not model accurately the highest eccentricity orbits. However, it is a simple matter to code the literal disturbing function for arbitrary inclinations with the steps we provided in the previous section using formal algebra software and easily obtain high-eccentricity power terms in a relatively short runtime on a personal computer.

The force coefficients $c_{mn}^k(p, q, \alpha, I_r)$ have an important property related to the resonance order $o_r = |p - q|$. Examination of $c_{mn}^k(p, q, \alpha, I_r)$'s dependence on $A_{i,j,k,l}$ and recalling that $A_{i,j,k,l} = 0$ when $j \pm k$ is odd show that for an odd resonance order $o_r$, $c_{mn}^k(p, q, \alpha, I_r) = 0$ when $k$ is even. Similarly, for an even resonance order $o_r$, $c_{mn}^k(p, q, \alpha, I_r) = 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.

These properties were encountered for the first time in the direct part of the disturbing function for nearly polar orbits (Paper I); they now appear as a general property of the disturbing function for arbitrary inclination and suggest that the classical disturbing function obtained with $I_r = 0^\circ$ is a degenerate two-dimensional limit of the general expansion. More on this aspect in Section 4.2.

The arguments and force amplitudes up to and including fourth order of the disturbing function's indirect part for an arbitrary reference inclination $I_r$ are given in Table B6 to fourth order of eccentricity $e$ and inclination $s$. As in the polar expansion, the indirect part concerns only cosines of the type 1:$q$ with $0 \le q \le 5$. For $q \ne 0$, arguments that include $q\lambda - \lambda'$ may induce resonance whereas those that start out as $q\lambda + \lambda'$ do not.

The secular potential to order N in eccentricity and inclination is obtained by setting p = q = 0 in the literal expansion while taking care not to count the same terms twice in the expansion as explained in the previous section. Its expression is found as

$$
\begin{equation}
\begin{split}
\bar{R}_s = \frac{1}{2} b_{\frac{1}{2}}^{00}(\alpha, I_r) &+ \sum_{\substack{k \le m \le N \\ 0 \le k, n \le N}} s_{mn}^k(\alpha, I_r) e^m s^n \cos(k\omega), \\
& \text{k,m even} \\
& \text{m+n=N}
\end{split}
\tag{22}
\end{equation}
$$

where in contrast to that of the polar disturbing function of Paper I,
the powers of inclination s are no longer limited to even integers.
This property is related to the symmetry of the expansion coeffi-
cients (and the Laplace coefficients) with respect to the polar plane
and is discussed further in Section 4 and Appendix A. The expres-
sions of the secular coefficients $s_{mn}^k(\alpha, I_r)$ are given in Table B7 for
the $N=4$ series expansion.

<footnote n="2">² The pair (0, 0) is a fixed point of the transformation (p → −p, q → −q).</footnote>
# 4 PROPERTIES OF THE DISTURBING FUNCTION WITH ARBITRARY REFERENCE INCLINATION

## 4.1 Comparison with the disturbing function for nearly polar orbits

The polar disturbing function of Paper I ($I_r = 90^\circ$) was the first disturbing function to be expanded with respect to a non-coplanar motion.<fnref n="3" /> It taught us that its properties differ significantly from those of the classical disturbing function whose $I_r = 0^\circ$ (Ellis & Murray 2000) and its companion retrograde expansion whose $I_r = 180^\circ$ (Morais & Namouni 2013a). The main differences were three. First, the literal expansion of order $N$ was applicable to any resonance ratio $p:q$ in contrast to the nearly coplanar disturbing function that produces terms only for the ratios $p:p \pm N$ (Murray & Dermott 1999). For example, the arguments $\phi_1^{3:10} = 10\lambda - 3\lambda' - 7\Omega - \omega$ and $\phi_3^{3:10} = 10\lambda - 3\lambda' - 7\Omega - 3\omega$ of the 3:10 resonance may be studied by the $N=4$ order series whereas with the classical disturbing function, in order to study the same arguments, the expansion order must be at least $N=7$. This is not a mere technical detail as the fact that the integers $p$ and $q$ are not related to one another through $N$ in the polar expansion is rooted in the presence of the two independent angles $f + \omega$ and $\Omega - \lambda'$ in the expression (9) that in turn gave rise to the two-dimensional Laplace coefficients, thus adding an extra degree of freedom to the problem. As already stated in Section 2, these two angles show up in the classical disturbing function only as the combination $f + \omega \pm (\Omega - \lambda')$ where $\pm$ refers to the expansion with $I_r = 0^\circ$ and $I_r = 180^\circ$, respectively, giving rise to the one-dimensional Laplace coefficients. The independence of the two angles found in the polar expansion is shared by the disturbing function with an arbitrary reference inclination $I_r$, provided that $I_r \neq 0^\circ$ and $I_r \neq 180^\circ$. The new disturbing function is indeed applicable to any resonance ratio $p:q$. The expansion order $N$ limits only the integer coefficient $k$ of the argument of pericentre in the resonance angle $\phi_k^{p:q} = q\lambda - p\lambda' + (p-q)\Omega - k\omega$ to the values $|k| \le N$ as can be seen in the expressions of the direct and secular parts of the disturbing function (20) and (22).

The next two properties found in the polar disturbing function are also shared by the new series expansion and originate in the properties of the coefficients $c_{mn}^k(p, q, \alpha, I_r)$. One is the fact that the powers of eccentricity $e$ and inclination $s$ in the force amplitudes of the new disturbing function are independent of the value of the resonance order and the other property is that only the parity of the resonance order decides the powers of $e$ and $s$ in the force amplitude of a $p:q$ resonance. As explained in Paper I, the former property does not apply for $I_r = 0^\circ$ and $I_r = 180^\circ$ since the powers of say the 1:4 resonance with $I_r = 0^\circ$ are of the form $e^{3-2m}\sin 2m(I/2)$ to lowest order in eccentricity and inclination. If one seeks the force amplitude linear with respect to $e$, one has to accept the presence of the inclination factor $\sin^2(I/2)$ in the force amplitude because the sum of the powers of eccentricity and inclination must equal the resonance order $|p-q|=3$. With the new disturbing function, there is a term of power 1 in eccentricity and power zero in inclination $s$ of coefficient $c_{10}^k(1, 4, \alpha, I_r)$ (Table B2) as long as the chosen $I_r \neq 0^\circ$ and $I_r \neq 180^\circ$. We discuss below in more precise terms

how the coefficients $c_{mn}^k(p, q, \alpha, I_r)$ behave when $I_r = 0^\circ$ and $I_r = 180^\circ$.

The property regarding the force amplitude and the parity of resonance order tells us that all odd order $p:q$ resonances have force amplitudes proportional to eccentricity as $(c_{10}^1 + c_{11}^1 s)e$ to lowest order in eccentricity and inclination. All even order $p:q$ resonances have force amplitudes of the form $c_{00}^0 + c_{01}^0 s + c_{20}^2 e^2 + c_{02}^2 s^2$ to lowest order in eccentricity and inclination. This universal binarity is also a property of the new disturbing function with arbitrary inclination $I_r$, provided $I_r \neq 0^\circ$ and $I_r \neq 180^\circ$. The dependence on the smallest powers of $e$ and $s$ (linear and quadratic) is realized with the smallest possible integer $|k|$ in the argument $\phi_k^{p:q} = q\lambda - p\lambda' + (p-q)\Omega - k\omega$.

The polar and arbitrary inclination disturbing functions share the previous three general properties and even have some terms with identical expressions such as the pure eccentricity coefficients $c_{k0}^k$ that do not depend explicitly on $I_r$ but only implicitly through the Laplace coefficients (10). The similarities however end there as the polar disturbing function is a special case of the general one for $I_r = 90^\circ$. For instance, the two-dimensional Laplace coefficients for an arbitrary $I_r$ do not share all the symmetries of their polar counterparts. Whereas for $I_r = 90^\circ$, $b_s^{j(k)} = b_s^{jk}$, the general Laplace coefficients only have reflection symmetry with respect to the polar plane as $b_s^{j(k)}(I_r) = b_s^{jk}(I_r + 180^\circ)$ modulo $180^\circ$ (see Appendix A). One of the salient consequences of this symmetry is the presence of terms linear with respect to inclination $s$ in the secular potential (equation 22 and Table B7) that are not present in the polar disturbing function ($I_r = 90^\circ$). More generally, one can easily check that the coefficients of Tables B1–B6 agree with those of Paper I for $I_r = 90^\circ$.<fnref n="4" />

The three properties of the polar expansion, generalized to the arbitrary inclination disturbing function signal different physical behaviours between expansions performed with respect to the two-dimensional three-body problem, namely the classical disturbing function and its companion for retrograde orbits, and expansions performed with respect to the three-dimensional three-body problem such as the series we derive in this work.

## 4.2 Relationship to the classical and retrograde disturbing functions

To examine further the relationship between the nearly coplanar disturbing functions and the three-dimensional ones, we consider the disturbing function's coefficients for $I_r = 0^\circ$ and $I_r = 180^\circ$. We note that for $I_r = 0^\circ$ or $180^\circ$, the inclination variable is $s = \pm\sin I$ where the $\pm$ sign is for the classical prograde and retrograde disturbing functions, respectively. However, the classical coplanar expansions are performed with respect to the variables $\sin(I/2)$ and $\cos(I/2)$. This slight difference does not affect our comparison as consistency requires that we use the classical expansions near $I \sim 0^\circ$ and $180^\circ$, respectively, which in turn implies that for the same cosine arguments, the classical disturbing function coefficients have extra $2^n$ factors for inclination terms with respect to those of the new disturbing function where $n$ is the inclination power in $s^n$. An example of this situation will be given below using an inclination term of the $I_r = 0^\circ$ disturbing function.

<footnote n="3">³ We consider only expansions in eccentricity and inclination but not semi-major axis employed historically for the motion of planetary satellites. See Laskar & Boué (2010) for a historical perspective and a modern presentation of such expansions.</footnote>

<footnote n="4">⁴ Note that in this paper, the inclination variable $s = -\cos I$ for $I_r = 90^\circ$, the opposite of that of Paper I. This means that all odd powers of $s$ will have an additional negative sign with respect to the coefficients derived in Paper I.</footnote>
Inspection of Tables B1–B7 shows that choosing $I_t = 0^\circ$ or $180^\circ$ eliminates all terms with an odd power of the inclination $s$ in the disturbing function’s direct part, indirect part as well as the secular potential. This recovers the classical dependence of the nearly coplanar series on powers of $s^2$.

In Appendix A, we find that for $I_t = 0^\circ$, the coefficients $A_{i,j,k,l} = 0$ unless $j = k$ (A13), and that for $I_t = 180^\circ$, the coefficients $A_{i,j,k,l} = 0$ unless $j = -k$ (A14). Applying these two properties to the terms $c_{k0}^k(p, q, \alpha, I_t = 0^\circ \text{ or } 180^\circ) e^{ikl} \cos \phi_k^{p,q}$ of Tables B1–B6 (for which $k > 0$) reveals that the integer $q$ is no longer independent of $p$ but must satisfy $q - p = k$ for the prograde expansion and $q + p = k$ for the retrograde one. For negative $k$, the force coefficients are obtained from equation (21) as $c_{k0}^k(p, -q, \alpha, I_t = 180^\circ \text{ or } 0^\circ)$ for the prograde and retrograde expansions, respectively, implying that $-q + p = k$ for the prograde expansion (that now obeys the retrograde rule $A_{i,j,k,l} = 0$ unless $j = -k$) and $-(q + p) = k$ for the retrograde expansion. This leads to the general requirement $|p - q| = |k|$ for the prograde expansion and $|p + q| = |k|$ for the retrograde one. We thus recover the expansion order to resonance order relationship of the classical expansions discussed at the beginning of this section.

It is important to emphasize that this relationship between the integers $p$ and $q$ of a cosine term in the disturbing function is specific to the disturbing functions of nearly coplanar orbits ($I_t = 0^\circ$ and $I_t = 180^\circ$). There are no similar relationships for $I_t \neq 0^\circ$ and $I_t \neq 180^\circ$. This indicates that if the general disturbing function is defined as a functional of the spatial dimension of the three-body problem that is chosen as a reference system with respect to which the series expansion is carried out, then the functional is extremal for dimension 2, as the relationship $q = p \pm N$ removes a significant number of terms from the expansion – that is why we termed the coplanar disturbing functions degenerate. It is beyond the scope of this paper to provide the mathematical proof for this conjecture.

In order to verify the agreement of the new disturbing function with the classical one, we recall that the two-dimensional Laplace coefficients $b_s^j(\alpha, I_t = 0^\circ) = 2b_s^j(\alpha)$ where $b_s^j(\alpha)$ is the classical (one-dimensional) Laplace coefficient ($A_{s1}$) (Appendix A). The prograde force amplitudes of the eccentricity terms discussed previously are then reduced to those of the classical disturbing function of Ellis & Murray (2000) as follows: $ec_{10}^1(j - k, j, \alpha, I_t = 0^\circ) = ef_{27}$ (term 4D1.1), $e^2c_{20}^2(j - k, j, \alpha, I_t = 0^\circ) = e^2f_{45}$ (term 4D2.1), $e^3c_{30}^3(j - k, j, \alpha, I_t = 0^\circ) = e^3f_{82}$ (term 4D3.1) and $e^4c_{40}^4(j - k, j, \alpha, I_t = 0^\circ) = e^4f_{90}$ (term 4D4.1).

The case of inclination terms is treated similarly with the difference that for a given inclination term, the resonances concerned are diverse. For instance, the term $c_{02}^0(p, q, \alpha, I_t = 0^\circ) s^2 \cos \phi_0^{p,q}$ has a force coefficient $c_{02}^0(p, q, \alpha, I_t = 0^\circ) = \alpha(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})/16$ and applies to one of the three resonances $q = p, q = p - 2$ and $q = p + 2$ as the integers $j$ and $k$ in each of the $A_{1,j,k,0}$ are different. For the latter two resonances, there is only one term in the coefficient’s expression that does not vanish. Explicitly, with $p = j + 2$ and $q = j$, we have $s^2c_{02}^0(j + 2, j, \alpha, I_t = 0^\circ) = \alpha s^2 A_{1,j-1,j-1,0}/16 = \alpha s^2 b_{j/2}^{j-1}/8 = f_{57}s^2/4$ where the latter coefficient is that of the corresponding term 4D2.4 of Ellis & Murray (2000) – an additional factor $1/4$ is present because of the choice of $s = \sin I$ instead of $\sin(I/2)$ as the inclination variable. More generally, it can be checked that the new fourth-order disturbing function with $I_t = 0^\circ$ agrees completely with the classical one (Ellis & Murray 2000).

## 4.3 Two interpretations of the arbitrary inclination disturbing function

The last property of the new disturbing function that we discuss is its versatility. Unlike the classical disturbing function, its retrograde companion and the polar disturbing function, the new series expansion has two interpretations depending on whether the inclination variable $s$ is set to zero or kept as a free variable (Section 2). In the latter case, the new disturbing function is a double series expansion with respect to eccentricity $e$ and inclination $s$. However, if $s$ is identically set to zero, the new disturbing function becomes an expansion with respect to eccentricity $e$ only and remains unexpanded with respect to inclination. The identity $s \equiv 0$ then means that the reference inclination is just the asteroid’s inclination. The dependence on inclination is transferred exclusively to the new two-dimensional Laplace coefficients.

The two interpretations of the new disturbing function have their advantages. For instance, when calculations are required to follow the time evolution of the asteroid’s dynamics, the double es expansion version is preferable as applying the Lagrange planetary equations to the *e*-expansion may prove difficult. Indeed, using the *e*-expansion with the Lagrange equations requires the evaluation of the Laplace coefficients and their derivatives with respect to inclination which in practice means the use of the *α*-expansion of the Laplace coefficients derived in Appendix A. To achieve good accuracy the *α*-expansion in turn requires a large order (see Appendix A for more details). For this reason, an es-expansion is preferable for numerical integrations as the Laplace coefficients depend only on the constant reference inclination and the semimajor axis ratio *α*. The version of the new disturbing function as an *e*-expansion is useful for analytical developments as we show in the next section with the derivation of mean motion resonance widths as a function of the asteroid’s inclination.

## 5 RESONANCE WIDTH AT ARBITRARY INCLINATION

We illustrate the use of the disturbing function for arbitrary inclination by deriving the resonance width of a number of pure eccentricity mean motion resonances. Such resonances have arguments that do not involve the longitude of ascending node explicitly and consequently force amplitudes that depend exclusively on eccentricity. The classical disturbing function that has always been used for their study affirms that pure eccentricity resonances whose arguments for prograde motion may be written as $q\lambda - p\lambda' + (p-q)\omega$ do not depend on the orbital inclination. In effect, only when the secular potential is added to the perturbing force, does the inclination affect the resonance albeit indirectly. Here, we show using the new disturbing function that pure eccentricity resonances do intrinsically depend on the inclination for both prograde and retrograde motion. This helps us demonstrate the importance of the new disturbing function in modelling three-dimensional resonance configurations accurately and caution against the general use of resonant dynamics properties obtained from the classical disturbing function that is based on the two-dimensional three-body problem.

In order to derive resonance widths as functions of inclination, we upgrade the pendulum model of resonance and apply it to specific examples that include inner resonances with Jupiter and outer resonances with Neptune. To this end, in the next section, we review the pendulum model then extend it to the presence of two harmonics so as to model the asymmetric librations of Neptune’s resonances.
## 5.1 Two-harmonics pendulum model of resonance

The pendulum model of resonance is the simplest method that allows the determination of the semimajor axis width of a mean motion resonance. Its drawback is the fact that it does not apply to orbits with very small eccentricity. The method is therefore adequate for our illustration needs as most asteroids, Centaurs and TNOs that interest us do not have negligible eccentricities. In the following, we briefly review how the pendulum model works and its relationship to the Poincaré model of resonance (Poincaré 1902).

Consider an asteroid in a p:q resonance with the planet of resonant argument $\phi_k^{p:q} = q\lambda - p\lambda' - (q - p)\Omega - k\omega$. We may write down the Lagrange planetary equation that governs the time evolution of the mean motion $n = (GM_*a^{-3})^{1/2}$ (Brouwer & Clemence 1961) exactly as

$$ \dot{n} = -3a^{-2}\partial_{\sigma}R = -\frac{3n^2\alpha m'}{M_*}\partial_{\sigma}\bar{R}, \quad (23) $$

where $\sigma$ is related to the mean longitude at epoch $\tau = \lambda - nt$ by $\dot{\sigma} = t\dot{n} + \dot{t}$ and $\bar{R}$ is the disturbing function with its direct, indirect and secular parts written as $\bar{R} = (f_{k,d}^{p,q} + f_{k,i}^{p,q}) \cos \phi_k^{p,q} + \bar{R}_s$ for the p:q resonance. The force amplitudes $f_{k,d}^{p,q}$ and $f_{k,i}^{p,q}$ include the eccentricity and inclination dependence of the resonant term. Next, consider the second derivative of the resonant argument with respect to time:

$$ \ddot{\phi}_k^{p,q} = q(\dot{n} + \ddot{\sigma}) - (q - p)\ddot{\Omega} - k\ddot{\omega}. \quad (24) $$

The reference to the mean motion of the planet has disappeared as the asteroid is considered massless leaving the planet's orbit unperturbed.

The working assumption of the pendulum model of resonance is to neglect, in the evolution of the resonant argument (24), the second derivatives $\ddot{\sigma}$, $\ddot{\Omega}$ and $\ddot{\omega}$ in comparison to the time variation of the mean motion $n$. This is generally valid unless the eccentricity is small because such angles as $\sigma$, $\Omega$ and $\omega$ would have rapid variations (Murray & Dermott 1999). Following the same logic, the secular term $\bar{R}_s$ has little effect on the resonance's dynamics as it does not influence the mean motion but concerns the discarded rates $\ddot{\sigma}$, $\ddot{\Omega}$ and $\ddot{\omega}$. The resonant state that the pendulum model describes is that of the stable branch of critical points that extends to large eccentricity in the Poincaré model of resonance (Murray & Dermott 1999). Since our main interest is the dynamics of asteroids, Centaurs and TNOs whose orbits are not circular, we may use the pendulum model assumption to study the resonance's dynamics by substituting (23) into (24) to obtain the following:

$$ \ddot{\phi}_k^{p,q} = \frac{3n^2 q^2 m' \alpha}{M_*} (f_{k,d}^{p,q} + f_{k,i}^{p,q}) \sin \phi_k^{p,q}. \quad (25) $$

This equation describes the evolution of the pendulum of angle $\phi_k^{p,q}$ and of natural frequency $|3n^2q^2m'M_*^{-1}\alpha(f_{k,d}^{p,q} + f_{k,i}^{p,q})|^{1/2}$. The sign of $f_{k,d}^{p,q} + f_{k,i}^{p,q}$ determines whether stable libration occurs around $\phi_k^{p,q} = 0^\circ$ when it is negative or around $\phi_k^{p,q} = 180^\circ$ when it is positive.

So far, we have considered a single harmonic of the p:q resonance given by the angle $\phi_k^{p,q}$. Modelling resonances with a single harmonic whether through the pendulum model or the Poincaré Hamiltonian is known to be sufficient for inner resonances. However, for outer resonances, higher harmonics' amplitudes modify the topology of phase space giving rise to asymmetric librations that depend on the orbit's eccentricity (Malhotra 1996; Winter & Murray 1997) (we will show in Section 5.3 that they are intrinsically dependent on inclination as well). This behaviour affects particularly outer

resonances of the type 1:q (Bruno 1994). For such resonances, we upgrade the pendulum model by including the second harmonic $\ddot{\phi}_{2k}^{p:2q} = 2\phi_k^{p:q}$ in a way similar to Andoyer's Hamiltonian model (Andoyer 1902; Beaugé 1994). The same restriction on the one-harmonic model applies here in that very small eccentricities may not be studied by the two-harmonics pendulum model. For outer resonances, the two-harmonics pendulum equation is given as

$$ \ddot{\phi}_k^{p,q} = \frac{3n^2 q^2 m' \alpha}{M_*} \left[ (f_{k,d}^{p,q} + f_{k,i}^{p,q}) \sin \phi_k^{p,q} + 2(f_{2k,d}^{p:2q} + f_{2k,i}^{p:2q}) \sin 2\phi_k^{p,q} \right]. \quad (26) $$

In order to derive the properties of this dynamical system, we write it in the following way $\ddot{\phi} = u \sin \phi + v \sin 2\phi$. We may choose $u > 0$ as when it is negative, the variable change $\phi \to \phi + 180^\circ$ makes the coefficient of $\sin \phi$ positive. With this convention, stable librations of the first harmonic ($v=0$) occur around $180^\circ$. Next, we change the time variable to $t' = |u|^{1/2}t$. The reduced two-harmonics pendulum equation then becomes $\phi'' = \sin \phi(1 + \beta \cos \phi)$ where $\phi''$ denotes the second derivative with respect to $t'$ and

$$ \beta = \frac{2v}{|u|} = \frac{4(f_{2k,d}^{p:2q} + f_{2k,i}^{p:2q})}{|f_{k,d}^{p,q} + f_{k,i}^{p,q}|}. \quad (27) $$

The corresponding Hamiltonian is given as

$$ H = \frac{p^2}{2} + \cos \phi + \frac{\beta}{4} \cos(2\phi), \quad (28) $$

where the momentum $p = \phi'$. The critical points of this system are given by $p = 0$ and $\sin \phi = 0$ and when possible by the asymmetric centres:

$$ \phi_a = \pm \arccos(-\beta^{-1}). \quad (29) $$

The Hamiltonian's level curves are shown in Fig. 1 for three values of the parameter $\beta$.

If $|\beta| < 1$, the position and stability of the two equilibrium points of the one-harmonic pendulum model are unchanged by the addition of the second harmonic. The motion separatrix passes through the unstable point $\phi = 0^\circ$ (as $u > 0$) and has an energy $H = 1 + \beta/4$. Therefore, the corresponding resonance half-width is given as $p(\phi = 180^\circ) = 2$ which in fully dimensional units reads $\dot{\phi} = 2|u|^{1/2}$ the same as that of the one-harmonic pendulum model ($\beta = 0$). Hence, if the second harmonic is weaker than a quarter of the first, libration stability and resonance width are unchanged with respect to the one-harmonic pendulum model. The modified properties are the libration time and elongation of the angle $\phi$.

When the second harmonic is stronger than a quarter of the first, $|\beta| \ge 1$, two possibilities must be considered $\beta < 0$ and $\beta > 0$. For the former ($\beta \le -1$), librations around $\phi = 180^\circ$ remain stable as can be seen directly from the reduced two-harmonics pendulum equation but in addition to them, stable librations around $\phi = 0^\circ$ are now possible. Two unstable points form at $\phi_{a,u\pm}(29)$. As the separatrix now passes through these points, the corresponding energy is given as $H = -\beta/4 - 1/2\beta$. This enables us to determine the maximum half-widths of stable librations around $0^\circ$ and $180^\circ$ as $p(\phi = 0^\circ) = |\beta|^{-1/2}|\beta + 1|$ and $p(\phi = 180^\circ) = |\beta|^{-1/2}|\beta - 1|$.

When $\beta \ge 1$, both one-harmonic critical points at $\phi = 0^\circ$ and $180^\circ$ become unstable and stable asymmetric librations appear around $\phi_{a,s\pm}(29)$. There are two separatrices in this case. One that goes through $\phi = 180^\circ$ and defines the maximum limit for asymmetric librations with an energy $H = \beta/4 - 1$. The second separatrix passes through $\phi = 0^\circ$ and delineates the domain of symmetric libration around both points $\phi_{a,s\pm}(29)$ from the general circulation domain;
<figure><img src="image_3.png" /><figcaption>Figure 1. Level curves of the two-harmonics pendulum Hamiltonian $H$ (28). Separatrices are shown with bold lines.</figcaption></figure>

its energy is $H = \beta/4 + 1$. The asymmetric libration half-width is $p(\phi_{a,s\pm}) = |\beta|^{-1/2}|\beta - 1|$ whereas that around both equilibria is $p(\phi_{a,s+} + \phi_{a,s-}) = |\beta|^{-1/2}|\beta + 1|$.

In order to find the resonance widths in terms of the semimajor axis, we recall that in the pendulum approximation $\dot{n} = q^{-1}\ddot{\phi}$. This leads to $n' = q^{-1}|u|^{1/2}\phi''$. Eliminating the time variable $t'$ in favour of $\phi$ leads to $dn/d\phi = q^{-1}|u|^{1/2}\phi''/\phi'$. Using the energy integral $H$ and integrating<fnref n="5" /> with respect to $\phi$ yields the general expression of the mean motion as a function of $\phi$ as

$$n = q^{-1}|u|^{1/2}[2(H - \cos \phi - \beta \cos(2\phi)/4)]^{1/2}. \quad (30)$$

The general expression of the mean motion can be used to derive the resonance width by applying it to a given separatrix defined by its energy and the unstable points it passes through and combining it with Kepler's third law $\Delta a = 2(3n)^{-1}\Delta n a$. Equivalently, noting that $dn = q^{-1}|u|^{1/2}dp$, the mean motion widths are proportional to the momentum widths derived previously through $\Delta n = q^{-1}|u|^{1/2}\Delta p$. Applying Kepler's third law yields the resonance widths in terms of semimajor axis.

Thus, when the second harmonic is weaker than a quarter of the first, $|\beta| < 1$, stable libration occurs around $\phi = 180^{\circ}$ with $H = 1 + \beta/4$ and the semimajor axis resonance width is given as

$$\Delta_0 a_k^{p,q} = \frac{4|u|^{1/2} a_k^{p,q}}{3qn} = \left[ \frac{16 \alpha m' |f_{k,d}^{p,q} + f_{k,i}^{p,q}|}{3M_*} \right]^{1/2} a_k^{p,q}, \quad (31)$$

where $a_k^{p,q} = \alpha a'$ is the resonant semimajor axis. This is the classical formula of the one-harmonic pendulum model (Murray & Dermott 1999) which is used for the critical arguments obtained from the classical disturbing function ($I_r = 0^{\circ}$) to visualize the resonance width as a function of eccentricity, assess how it increases with eccentricity and determine the possible resonance overlapping that triggers chaotic motion.

When the second harmonic is stronger than a quarter of the first and its amplitude is negative, $\beta \le -1$, stable librations may occur around $\phi = 0^{\circ}$ and $180^{\circ}$, with a larger width for the latter. There is a single separatrix corresponding to the energy $H = -\beta/4 - 1/2\beta$ and the semimajor axis width at $\phi = 0^{\circ}$ is given as

$$\Delta_1 a_k^{p,q} = \frac{2|2v + |u||a_k^{p,q}|}{3qn|2v|^{1/2}},$$

$$\Delta_1 a_k^{p,q} = \left[ \frac{\alpha m'}{3M_*} \right]^{1/2} \frac{4(f_{2k,d}^{2p:2q} + f_{2k,i}^{2p:2q}) + |f_{k,d}^{p,q} + f_{k,i}^{p,q}|}{|f_{2k,d}^{2p:2q} + f_{2k,i}^{2p:2q}|^{1/2}} a_k^{p,q}. \quad (32)$$

Libration around $\phi = 180^{\circ}$ occurs with maximum semimajor axis width:

$$\begin{align}
\Delta_2 a_k^{p,q} &= \frac{2|2v - |u||a_k^{p,q}|}{3qn|2v|^{1/2}}, \nonumber \\
\Delta_2 a_k^{p,q} &= \left[ \frac{\alpha m'}{3M_*} \right]^{1/2} \frac{4(f_{2k,d}^{2p:2q} + f_{2k,i}^{2p:2q}) - |f_{k,d}^{p,q} + f_{k,i}^{p,q}|}{|f_{2k,d}^{2p:2q} + f_{2k,i}^{2p:2q}|^{1/2}} a_k^{p,q}. \tag{33}
\end{align}$$

When the second harmonic is stronger than a quarter of the first and its amplitude is positive, $\beta \ge 1$, asymmetric librations at either $\phi_{s\pm} = \pm \arccos(-\beta^{-1})$ are possible and are delineated by a separatrix of energy $H = \beta/4 - 1$, the resonance's semimajor axis width is given by $\Delta_2 a_k^{p,q}$. A second separatrix is associated with symmetric librations about both $\phi_{s\pm}$ centred around $\phi = 180^{\circ}$ and of energy $H = \beta/4 + 1$. The corresponding resonance width is given by $\Delta_1 a_k^{p,q}$.

All this analysis was done for a positive first harmonic amplitude ($f_{k,d}^{p,q} + f_{k,i}^{p,q} > 0$). If that amplitude is negative, then the only change that occurs is that all critical points have to be shifted by $180^{\circ}$ as the variable change $\phi = 180^{\circ} + \bar{\phi}$ produces a two-harmonics pendulum equation for $\bar{\phi}$ identical to that of $\phi$ except for a positive sign for the first harmonic coefficient. Thus, for instance, stable asymmetric librations for $f_{k,d}^{p,q} + f_{k,i}^{p,q} < 0$ occur at $\phi_{a,s\pm} = 180^{\circ} \pm \arccos(-\beta^{-1})$ when $\beta \ge 1$. Examples of the use of the two-harmonics pendulum model are given in Section 5.3.2 for Neptune's outer 1:2 and 1:3 resonances.

## 5.2 Prograde and retrograde resonant arguments

The disturbing function of an asteroid that interacts with a planet on a circular orbit involves the argument $\phi_k^{p,q} = q\lambda - p\lambda' - (q - p)\Omega - k\omega$ with the usual orbital angles defined for prograde motion. To emphasize the prograde nature of the chosen angles, we write the argument as a function of the longitude of pericentre, $\omega$, instead of the argument of pericentre, $\omega'$, as the former's definition changes when motion is retrograde:

$$\phi_+ = q\lambda - p\lambda' - k\omega + (p - q + k)\Omega. \quad (34)$$

For a given $(p, q, k)$ integer set, we can follow the dependence of the corresponding resonance width on the reference inclination even when $I_r$ exceeds $90^{\circ}$ as it is known that libration in prograde

<footnote n="5">⁵ Denoting $V(\phi)$ the perturbing potential $\cos \phi + \beta \cos(2\phi)/4$, we have $\phi' = [2(H - V)]^{1/2}$ and $\phi'' = -dV/d\phi$. The quantity $dn/d\phi = -|u|^{1/2}q^{-1}(dV/d\phi)/[2(H - V)]^{1/2}$.</footnote>
resonance exists at retrograde inclinations albeit with a modified dynamical behaviour. Such libration was observed when retrograde motion is about to be captured in mean motion resonance in a universal process termed as the three-stage capture process of high-order resonance (Namouni & Morais 2015). It may be summed up briefly as follows: capture first occurs in the prograde resonant argument that behaves like an inclination resonance. Motion enters the second stage when the Kozai–Lidov resonance is activated followed by the final capture in the retrograde resonance termed as the third stage. Examples of this mechanism are given in Namouni & Morais (2015) for the 1:2 resonance and an inclination $I = 170^{\circ}$ and in Morais & Namouni (2016) for the 1:1 resonance with an inclination $I = 179.999$.

We may also use the argument $\phi_k^{p,q}$ to study retrograde resonance by adopting the angles $\lambda^*$ and $\omega^*$ that are compatible with retrograde motion and related to the mean longitude through $\lambda = \lambda^* + 2\Omega$ and to the longitude of perihelion through $\omega = \omega^* + 2\Omega$ (Morais & Namouni 2013a). The retrograde resonant argument thus reads

$$ \phi_{-} = q\lambda^{*} - p\lambda' - k\omega^{*} + (q + p - k)\Omega. \quad (35) $$

It is evident from the expressions (34) and (35) that what defines the prograde or retrograde flavour of a $p:q$ resonance is the choice of the integer $k$. For instance, if we are interested in the prograde pure eccentricity resonant term (34), then $k = q - p$ whereas if we go for the equivalent retrograde pure eccentricity resonant term (35), $k = q + p$. In this regard, we remark that for prograde pure eccentricity arguments, the integer $k$ depends on the perturber's location: inside perturbers have positive $k$ whereas outside perturbers have negative $k$. Retrograde pure eccentricity arguments do not distinguish between outside and inside perturbers as $k$ is always positive. In the following, we estimate prograde as well as retrograde resonance widths.

## 5.3 Resonance widths

The resonance width is a function of eccentricity $e$, inclination $s$ and reference inclination $I_r$ through the force amplitudes $f_{k,d}^{p,q}$ and $f_{k,i}^{p,q}$ (equations 31–33). In this illustration of the use of the new disturbing function, we seek to understand how pure eccentricity resonance widths depend on the inclination of the resonant motion of an asteroid. The classical disturbing function being a series expansion with respect to circular coplanar orbits predicts that pure eccentricity terms do not depend on the asteroid’s orbital inclination. They may depend indirectly on the inclination variations only through the secular part of the disturbing function that couples the evolution of eccentricity to that of inclination but the inner workings of such resonances remain inclination free. The fact that the force amplitudes of the new disturbing function depend on the reference inclination $I_r$ through the two-dimensional Laplace coefficients will invalidate this conclusion.

As we need to select a meaningful reference inclination $I_r$ for a resonant asteroid with an inclination $I_r$, we can invoke the two interpretations of the new disturbing function to this effect. If we choose to cancel the free inclination variable through $s \equiv 0$ and retain only the eccentricity expansion, then the reference inclination in the disturbing function is just the asteroid’s inclination and is present in the two-dimensional Laplace coefficients. Since the resonance width formulae in the pendulum approximation assume exact resonance for $e$ and $I_r$, these elements are actually constant. With the interpretation of the new disturbing function as a double

es-expansion, the reference inclination can be set to the asteroid’s mean inclination $I_r = \langle I \rangle$, that is, its value after short period variations are averaged out. More precisely, the reference inclination can be set to the asteroid’s proper inclination (Knežević & Milani 2000). Choosing those orbits whose inclination amplitude $s$ is small compared to the mean or proper value $\langle I \rangle$ gives us the resonance widths of the previous choice with $s \equiv 0$. Therefore, in the context of the pendulum model of resonance, whether the reference inclination is interpreted as the asteroid’s inclination, or the proper value, the derived resonance widths are identical.

For the asteroid’s orbital eccentricity, we choose three values $e = 0.1, 0.3$ and $0.5$. The classical studies of resonance dynamics indicate that resonance width is an increasing function of eccentricity when orbits are not nearly circular (Murray & Dermott 1999).

### 5.3.1 Inner 2:1 and 3:1 Jupiter eccentricity resonances

The appropriate resonant argument for the pure eccentricity prograde inner 2:1 resonance is obtained with $k = -1$ as $\phi_{-1}^{2:1} = \lambda - 2\lambda' + \omega$. As the inclination’s oscillations are assumed negligible ($s \ll 1$), we neglect all inclination terms in the force amplitudes (proportional to powers of $s$). To fourth order in eccentricity, the amplitude thus involves the coefficients $c_{10}^{1}(2, -1, \alpha = 2^{-2/3}, I_r + 180^{\circ})$ and $c_{30}^{1}(2, -1, \alpha = 2^{-2/3}, I_r + 180^{\circ})$ in Table B2 where we used the relationship (21) given in Section 3 between coefficients $c_{mn}^k$ of positive and negative $k$. The force amplitude expression reads

$$ f_{-1,d}^{2:1} = -\frac{e}{4}(4\bar{A}_{0,2,-2,0} + \bar{A}_{0,2,-2,1}) + \frac{e^3}{32}(28\bar{A}_{0,2,-2,0} + 5\bar{A}_{0,2,-2,1} - 6\bar{A}_{0,2,-2,2} - \bar{A}_{0,2,-2,3}), \quad (36) $$

where $\bar{A}_{i,j,k,l}(I_r) = A_{i,j,k,l}(I_r + 180^{\circ})$. We may remove the negative signs and the bar from the previous expression as $\bar{A}_{i,j,-k,l}(I_r) = \bar{A}_{i,j,k,l}(I_r + 180^{\circ}) = A_{i,j,k,l}(I_r + 360^{\circ}) = A_{i,j,k,l}(I_r)$ (see equation A8). For inner resonances, there is no contribution from the indirect part of the disturbing function.

Similarly, the resonant argument for the pure eccentricity prograde inner 3:1 resonance is obtained with $k = -2$ as $\phi_{-2}^{3:1} = \lambda - 3\lambda' + 2\omega$. The corresponding force amplitudes to fourth order in eccentricity are obtained from the relationship (21) as $c_{20}^2(3, -1, \alpha = 3^{-2/3}, I_r + 180^{\circ})$ and $c_{40}^2(3, -1, \alpha = 3^{-2/3}, I_r + 180^{\circ})$ (Table B3):

$$ f_{-2,d}^{3:1} = \frac{e^2}{16}(21A_{0,3,3,0} + 10A_{0,3,3,1} + A_{0,3,3,2}) + \frac{e^4}{192}(15A_{0,3,3,2} - 186A_{0,3,3,0} - 122A_{0,3,3,1} + 12A_{0,3,3,3} + A_{0,3,3,4}), \quad (37) $$

where like the 2:1 expression we removed the bar and negative signs in $A_{i,j,k,l}$ and set their inclinations to $I_r$ instead of $I_r + 180^{\circ}$.

In order to estimate the resonances’ widths, we use the one-harmonic pendulum model as it is most appropriate for inner resonances. Throughout this section and the next, the values of the inclination-dependent two-dimensional Laplace coefficients are determined using the series derived in Appendix A (A29) with an $\alpha$-expansion order $N_\alpha = 20$. Fig. 2 shows the resonance widths obtained by inserting the expressions (36) and (37) into the resonance width formula (31) as a function of the asteroid’s inclination and the three possible eccentricities $0.1, 0.3$ and $0.5$. Jupiter’s parameters are semimajor axis $a' = 5.2$ au, mass $m'/M_* = 10^{-3}$. The widths of
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Inner 2:1 and 3:1 Jupiter pure eccentricity resonance widths (31) as functions of the asteroid's average inclination $I_r$ for three eccentricity values $e = 0.1, 0.3$ and $0.5$. Green solid lines indicate libration of the resonant argument around $180^\circ$ and red dashed lines around $0^\circ$. The dashed vertical line denotes the location of nominal resonance. The number next to each curve is the corresponding eccentricity value.</figcaption></figure>

prograde inner resonances are found to be monotonically decreasing functions of inclination that vanish for exactly retrograde coplanar motion. The eccentricity dependence agrees with the classical result that larger eccentricity implies larger width.

The appropriate resonant arguments for the pure eccentricity retrograde inner 2:1 and 3:1 resonances are obtained with $k = 3$ and $k = 4$, respectively, as $\phi_3^{3:1} = \lambda^* - 2\lambda' - 3\omega^*$ and $\phi_4^{3:1} = \lambda^* - 3\lambda' - 4\omega^*$. The corresponding force amplitudes to fourth order in eccentricity are found in Tables B4 and B5 for the direct part $c_{30}^3(2, 1, \alpha = 2^{-2/3}, I_r)$ and $c_{40}^4(3, 1, \alpha = 3^{-2/3}, I_r)$ as

$$f_{3,d}^{2:1} = \frac{e^3}{96}(4A_{0,2,-2,0} - 3A_{0,2,-2,1} - 6A_{0,2,-2,2} - A_{0,2,-2,3}), \quad (38)$$

$$f_{4,d}^{3:1} = \frac{e^4}{768}(30A_{0,3,-3,2} - 15A_{0,3,-3,0} - 4A_{0,3,-3,1} + 12A_{0,3,-3,3} + A_{0,3,-3,4}), \quad (39)$$

where unlike the prograde expressions (36) and (37), we may not remove the negative signs in $A_{i,j,k,l}$ unless we change the inclination from $I_r$ to $I_r + 180^\circ$ (see equation A8). It is interesting to note that the expressions (36)–(39) are identical to those obtained with nearly coplanar motion, that is, $j = k$ in $A_{i,j,k,l}$ for prograde motion and $j = -k$ for retrograde motion, except for the fact that the new Laplace coefficients in $A_{i,j,k,l}$ now depend on the reference inclination. The widths of the 2:1 and 3:1 retrograde Jupiter resonances obtained from expression (31) are shown in Fig. 2 where they are found to be naturally narrower than their
prograde counterparts. The retrograde 2:1 resonance width does not decrease monotonically and has a maximum near ~40° whereas both retrograde resonance widths vanish for exact prograde coplanar motion.

We therefore conclude that pure eccentricity inner resonances whether prograde or retrograde depend intrinsically on the asteroid's inclination and therefore will influence the resonance's capture efficiency. This conclusion is contrary to what may be inferred from the classical disturbing function whose reference inclination is $I_r = 0°$.

### 5.3.2 Outer 1:2 and 1:3 Neptune eccentricity resonances

In the two-dimensional three-body problem, outer resonances of the type 1:q are known to exhibit asymmetric libration whose resonant arguments' equilibria lie somewhere between 50° and 180° depending on the asteroid's orbital eccentricity (Beaugé 1994; Bruno 1994; Malhotra 1996; Winter & Murray 1997). The origin of such equilibria for a p:q resonance of argument $\phi_k^{p:q}$ is related to the importance of the second harmonic's $\phi_{2k}^{2p:2q}$ amplitude.<fnref n="6" /> For this reason, we developed the new two-harmonics pendulum model and proceed to apply it to Neptune's outer resonances.

The pure eccentricity prograde outer 1:2 resonant argument is obtained with $k=1$ as $\phi_1^{1:2} = 2\lambda - \lambda' - \omega$. The corresponding force amplitudes to fourth order in eccentricity are found in Table B2 for the direct part: $c_{10}^1(1, 2, \alpha = 2^{2/3}, I_r)$ and $c_{30}^1(1, 2, \alpha = 2^{2/3}, I_r)$, and in Table B6 for the indirect part. The force amplitude of the second harmonic $\phi_2^{2:4} = 4\lambda - 2\lambda' - 2\omega$ is read off Table B3 ($k=2$) with the coefficients $c_{20}^2(2, 4, \alpha = 2^{2/3}, I_r)$ and $c_{40}^2(2, 4, \alpha = 2^{2/3}, I_r)$ (no indirect terms are applicable). We thus have

$$f_{1,d}^{1:2} = \frac{e}{4}(2A_{0,1,1,0} - A_{0,1,1,1}) + \frac{e^3}{32}(14A_{0,1,1,1} - 20A_{0,1,1,0} - A_{0,1,1,3}), \quad (40)$$

$$f_{1,i}^{1:2} = \frac{e\alpha}{16}(3e^2 - 4)(1 + \cos I_r), \quad (41)$$

$$f_{2,d}^{2:4} = \frac{e^2}{16}(26A_{0,2,2,0} - 10A_{0,2,2,1} + A_{0,2,2,2}) + \frac{e^4}{192}(428A_{0,2,2,1} - 1036A_{0,2,2,0} - 30A_{0,2,2,2} - 8A_{0,2,2,3} + A_{0,2,2,4}). \quad (42)$$

Although there is no explicit inclination dependence in the direct force amplitude, we remind the reader that the coefficients $A_{i,j,k,l}$ (12) are functions of inclination through the two-dimensional Laplace coefficients $b_s^{jk}$ (10) like in the previous examples of Jupiter's inner resonances. The amplitudes are evaluated numerically for Neptune with a semimajor axis $a' = 30.11$ au and mass $m'/M_* = 5.12 \times 10^{-5}$.

When estimating the resonance width with the two-harmonics model, two quantities must be monitored, the ratio $\beta = 4f_{2,d}^{2:4}/|f_{1,d}^{1:2} + f_{1,i}^{1:2}|$ (27) of the second harmonic's amplitude to a quarter of the first harmonic's absolute amplitude, as well as the sign of the first harmonic's amplitude $\text{sgn}1^2 = \text{sign}(f_{1,d}^{1:2} + f_{1,i}^{1:2})$.

The 1:2 resonance's $\beta$ is found to be always positive. For $e=0.1$, the ratio $\beta \ge 1$ for $0° \le I_r \le 129°$ implying the presence of asymmetric libration up to that inclination and simple 180° libration beyond

it. The sign $\text{sgn}1^2 < 0$ in the inclination interval $[37°, 96°]$ and positive otherwise. The sign change allows the asymmetric libration to exist below the value $\phi_1^{1:2} = 90°$ that is reached at the edges of that interval where the separatrices of asymmetric libration and 180° libration intersect. The resonance widths are thus obtained from equations (32) and (33) for $I_r \le 129°$ and from equation (31) beyond that value. The resulting resonant widths as function of the inclination are shown in Fig. 3. It is found that the resonance width of asymmetric libration shows a slight increase starting from the planar configuration to reach a maximum at 36° whereas the 180° libration about both asymmetric centres increases and reaches a minimum at that inclination. A second extremum is reached at 96° after which asymmetric libration width decreases until it disappears in favour of simple 180° libration ($0° < \beta < 1°$). Fig. 4 shows the location of the positive asymmetric libration centre (29); the negative one is symmetric with respect to $\phi_1^{1:2} = 0°$. The asymmetric libration centre starts out at $\phi_1^{1:2} = 108°$ for $I_r = 0°$, reaches a minimum $\phi_1^{1:2} = 75°$ for $I_r = 69°$ then increases steadily to 180° and disappears at $I_r = 129°$.

Increasing the eccentricity value to $e=0.3$ increases the resonance widths and produces a similar behaviour except for the asymmetric libration width that decreases as the asteroid's orbit leaves the planet's orbital plane. The three critical inclinations are $I_r = 36°, 96°$ where the libration centres reach 90° and $I_r = 155°$ beyond which asymmetric libration no longer exist. The location of asymmetric libration at $I_r = 0°$ is $\phi_1^{1:2} = 98°$ and the minimum of 83° is reached at $I_r = 68°$. Regardless of eccentricity, all resonant libration disappear for exactly coplanar retrograde motion $I_r = 180°$.

We now examine the pure eccentricity prograde outer 1:3 resonance. Its argument is obtained with $k=2$ as $\phi_2^{1:3} = 3\lambda - \lambda' - 2\omega$ and the corresponding force amplitudes are found in Table B3 for the direct part: $c_{20}^2(1, 3, \alpha = 3^{2/3}, I_r)$, $c_{40}^2(1, 3, \alpha = 3^{2/3}, I_r)$, and in Table B6 for the indirect part. The second harmonic $\phi_4^{2:6} = 6\lambda - 2\lambda' - 4\omega$ has the $k=4$-amplitude coefficient $c_{40}^4(2, 6, \alpha = 3^{2/3}, I_r)$ found in Table B5 (no indirect terms are applicable). The force amplitudes are given explicitly as

$$f_{2,d}^{1:3} = \frac{e^2}{16}(9A_{0,1,1,0} - 6A_{0,1,1,1} + A_{0,1,1,2}) + \frac{e^4}{192}(126A_{0,1,1,1} - 162A_{0,1,1,0} - 21A_{0,1,1,2} - 4A_{0,1,1,3} + A_{0,1,1,4}), \quad (43)$$

$$f_{2,i}^{1:3} = \frac{3e^2\alpha}{16}(e^2 - 1)(1 + \cos I_r), \quad (44)$$

$$f_{4,d}^{2:6} = \frac{e^4}{768}(2760A_{0,2,2,0} - 1464A_{0,2,2,1} + 300A_{0,2,2,2} - 28A_{0,2,2,3} + A_{0,2,2,4}). \quad (45)$$

For the planar 1:3 resonance, asymmetric libration are triggered at a larger eccentricity value than the 1:2 resonance, a threshold estimated numerically at $e \sim 0.13$ (Beaugé 1994; Malhotra 1996). It is therefore not surprising that for an eccentricity $e=0.1$, the ratio $0° < \beta < 1°$ and asymmetric libration are absent for $I_r = 0°$ (Fig. 3). At the same eccentricity but for three-dimensional motion, two small asymmetric libration islands appear in the inclination intervals $[15°, 62°]$ and $[78°, 106°]$ (Fig. 3). The position of the positive libration centre decreases with inclination from 180° to 0° in the first interval and increases with inclination from 0° to 180°. In the same inclination intervals, 180° libration around both

<footnote n="6">⁶ Historically, it was thought that at least nine harmonics were required to establish asymmetric libration (Message 1958).</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Outer 1:2 and 1:3 Neptune pure eccentricity resonance widths as functions of the asteroid’s average inclination $I_r$ for two eccentricity values $e = 0.1$ and $0.3$. Solid green lines indicate the asymmetric resonance widths (33), dashed red lines indicate the librations around $180^\circ$ encompassing both asymmetric centres (32) and dotted blue lines indicate simple libration around $180^\circ$ when the second harmonic is weaker than a quarter of the first (31).</figcaption></figure>

asymmetric centres is naturally possible. Otherwise, all libration occurs around the centre defined by the first harmonic namely $180^\circ$ (as $0 < \beta \le 1$). For the larger eccentricity $e = 0.3$, asymmetric librations exist for coplanar prograde motion and the librations’ structure and width are similar to that of 1:2 resonance with the three critical inclination values $I_r = 38^\circ, 96^\circ$ and $138^\circ$ (Fig. 4). Regardless of eccentricity, all resonant librations disappear for exactly coplanar retrograde motion $I_r = 180^\circ$.

Lastly, the 1:2 and 1:3 retrograde Neptune resonances cannot be examined by the current fourth-order disturbing function as their second harmonics given by $\phi_6^{2:4}$ and $\phi_8^{2:6}$, respectively, require expansions of the order of $N = 6$ and $N = 8$.<fnref n="7" /> We remind

the reader that an expansion of order $N$ models all resonant angles $\phi_k^{p:q} = q\lambda - p\lambda' - (q-p)\Omega - k\omega$ with $|k| \le N$ regardless of the values of $p$ and $q$ – see, for instance, the expression of the direct part (20).

## 6 CONCLUSION

This work’s main motivation of understanding the processes of resonance passage and capture at arbitrary inclination was inspired by our earlier numerical work on the subject (Namouni & Morais 2015, 2017a,c), and encouraged by the realization that the polar disturbing function of Paper I challenged the standard view on the force amplitudes associated with a resonant argument. We had therefore wondered whether the new properties of the polar disturbing function signalled a new general behaviour of

<footnote n="7">⁷ We used the property that for pure eccentricity retrograde resonances $k = p + q$ in Section 5.2.</footnote>
<figure><img src="image_6.png" /></figure>

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

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

<figure><img src="image_9.png" /><figcaption>**Figure 4.** Outer 1:2 and 1:3 Neptune positive asymmetric resonance centres (29) as a function of the asteroid's average inclination $I_r$ for two eccentricity values $e = 0.1$ and 0.3.</figcaption></figure>

disturbing functions that are based on the three-dimensional three-body problem. This work answers that query with the affirmative. We did so by providing the algorithm for the literal expansion of the gravitational interaction in the three-body problem with respect to an arbitrary inclination $I_r$ and giving explicitly its terms to fourth order in eccentricity $e$ and inclination $\sin(I - I_r)$. We then showed how the new disturbing function reduces to the classical one by taking $I_r = 0°$. The properties of the classical disturbing function once considered standard henceforth seem those of a degenerate expansion embodied by (i) the presence of the one-dimensional Laplace coefficients, (ii) the constraint that an expansion of order $N$ can only model resonances of order $N$ and (iii) the fact that only a pure eccentricity resonance of order $k$ has, to the lowest order in eccentricity, a force amplitude $\propto e^{|k|}$ and no inclination dependence. These once-thought fundamental properties are not true in general. When expansions are not carried out with respect to the two-dimensional problem, i.e. $I_r \neq 0°$ and $I_r \neq 180°$ in the first

interpretation of the new disturbing function, or not carried out at all in the second interpretation of the new disturbing function with $s \equiv 0$, then the expansion of order $N$ may model any resonance order. More importantly, the force amplitudes of a resonant term are mostly independent of the resonance order giving, for instance the argument $\phi_1^{1:6} = 6\lambda - \lambda' - 5\Omega - \omega$, a force amplitude linear with respect to eccentricity that can be modelled by a series expansion of first order in eccentricity. To illustrate further the far-reaching consequences of the new disturbing function, we determined the inclination dependence of pure eccentricity resonances by measuring their widths with a new two-harmonics pendulum model. We found that contrary to the prediction of the classical disturbing function, those resonances depend strongly on inclination not only in terms of width but also in the nature of libration in resonance. For instance, we showed that below the threshold eccentricity of 0.13 known to exclude asymmetric librations for the pure eccentricity 1:3 outer resonance with Neptune, asymmetric librations do exist
around inclinations of 45° and 90°. That eccentricity threshold was derived numerically for the planar problem and thought to be valid in general as pure eccentricity resonances were believed to be in-dependent of inclination. We now know that the 0.13 threshold is valid only for planar motion. Our preliminary numerical estimations (presented in a forthcoming article) using chaos indicators agree with the width and libration structures shown in Figs 2 and 3. There are two natural extensions of this work: first, the derivation of a disturbing function with an eccentric and inclined perturber. Our original motivation being the dynamical study of small Solar system bodies with arbitrary orbital inclinations, our approximation of a circular planet and our choice of its orbital plane as an inclination reference are quite sufficient. Indeed, we have shown in Namouni & Morais (2015) that typical Solar system planetary eccentricities do not have a significant influence on mean motion resonance capture. Furthermore, as the inclinations of the Solar system planets are small with respect the invariable Laplace surface (Tremaine, Touma & Namouni 2009), they will not influence significantly the dynamics of asteroids with large inclinations. The second extension consists of applying the present disturbing function to derive analytically resonance capture probabilities at arbitrary inclination.

## ACKNOWLEDGEMENTS

We thank Zoran Knežević for his thorough and helpful review of the manuscript. MHMM was supported by grant 2015/17962-5 of São Paulo Research Foundation (FAPESP).

## REFERENCES

Abramowitz M., Stegun I. A., 1972, Handbook of Mathematical Functions. Dover Press, New York

Andoyer H., 1902, Bull. Astron. Ser. I, 20, 321

Beaugé C., 1994, Celest. Mech. Dyn. Astron., 60, 225

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

Bruno A. D., 1994, The Restricted 3-Body Problem: Plane Periodic Orbits. De Gruyter, Berlin

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

Knežević Z., Milani A., 2000, Celest. Mech. Dyn. Astron., 78, 17

Knežević Z., Froeschlé C., Lemaître A., Milani A., Morbidelli A., 1995, A&A, 293, 605

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

Malhotra R., 1996, AJ, 111, 504

Message P. J., 1958, AJ, 63, 443

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

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

Morais M. H. M., Namouni F., 2016, Celest. Mech. Dyn. Astron., 125, 91

Morais M. H. M., Namouni F., 2017, MNRAS, 472, L1

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., 2017a, MNRAS, 467, 2673

Namouni F., Morais M. H. M., 2017b, MNRAS, 471, 2097 (Paper I)

Namouni F., Morais M. H. M., 2017c, JCAM, preprint (arXiv:1704.00550)

Poincaré H., 1902, Bull. Astron. Ser. I, 19, 289

Tremaine S., Touma J., Namouni F., 2009, AJ, 137, 3706

Wiegert P., Connors M., Veillet C., 2017, Nature, 543, 687

Williams J. G., 1969, PhD thesis, Univ. California

Winter O. C., Murray C. D., 1997, A&A, 328, 399

## APPENDIX A: TWO-DIMENSIONAL LAPLACE COEFFICIENTS WITH ARBITRARY REFERENCE INCLINATION

In the following, we examine the properties of the new two-dimensional Laplace coefficients that generalize the classical one-dimensional Laplace coefficients of the classical disturbing function of nearly coplanar orbits (Brouwer & Clemence 1961), and the two-dimensional Laplace coefficients of the disturbing function of nearly polar orbits (Paper I). The new coefficients depend on the arbitrary reference inclination $I_r$ of the disturbing function and are defined as

$$b_s^{jk}(\alpha, I_r) = \frac{1}{\pi^2} \int_0^{2\pi} \int_0^{2\pi} \cos(ju + kv) \, du \, dv \\ \times [1 + \alpha^2 - 2\alpha(\cos u \cos v - \sin u \sin v \cos I_r)]^{-s}. \quad (A1)$$

### A1 Intrinsic properties

We first recall the important property:

$$b_s^{jk}(\alpha, I_r) = 0 \text{ if } j \pm k \text{ is odd.} \quad (A2)$$

The Laplace coefficients are defined for both $\alpha < 1$ and $\alpha > 1$ and satisfy the symmetry relation:

$$b_s^{jk}(p, q, \alpha, I_r) = \alpha^{-2s} b_s^{jk}(p, q, \alpha^{-1}, I_r). \quad (A3)$$

This property is useful when seeking a series expansion with respect to $\alpha > 1$ (see Section A3).

The coefficients $b_s^{jk}(\alpha, I_r)$ are not completely symmetric with respect to the angle integer coefficients $j$ and $k$ (as it is the case for polar motion i.e. $I_r = 90^\circ$). Instead, we have the following identities that depend on the reference inclination $I_r$:

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

$$\begin{aligned} b_s^{(-j)k}(\alpha, I_r) &= b_s^{j(-k)}(\alpha, I_r), \\ &= b_s^{jk}(\alpha, I_r + 180^\circ). \end{aligned} \quad (A5)$$

In the last identity, the sum $I_r + 180^\circ$ is defined modulo $180^\circ$.<fnref n="8" /> The previous properties naturally extend to the functions $A_{i,j,k,l} = \alpha^l D^l b_{i+1/2}^{jk}$ (equation 12) as

$$A_{i,j,k,l}(\alpha, I_r) = 0 \text{ if } j \pm k \text{ is odd,} \quad (A6)$$

$$\begin{aligned} A_{i,-j,-k,l}(\alpha, I_r) &= A_{i,j,k,l}(\alpha, I_r), \\ &= A_{i,k,j,l}(\alpha, I_r), \end{aligned} \quad (A7)$$

$$\begin{aligned} A_{i,-j,k,l}(\alpha, I_r) &= A_{i,j,-k,l}(\alpha, I_r), \\ &= A_{i,j,k,l}(\alpha, I_r + 180^\circ). \end{aligned} \quad (A8)$$

The derivatives of the Laplace coefficients also depend on the reference inclination as follows:

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

<footnote n="8">⁸ One could equally write $I_r - 180^\circ$ instead of $I_r + 180^\circ$ in equation (A5).</footnote>
$$
\begin{equation}
\begin{aligned}
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)}) \right. \\
& \left. \times (1 + \cos I_r) \right. \\
& \left. + (D^{n-1} b_{s+1}^{(j-1)(k+1)} + D^{n-1} b_{s+1}^{(j+1)(k-1)}) \right. \\
& \left. \times (1 - \cos I_r) \right] - 2\alpha s D^{n-1} b_{s+1}^{jk} \\
& - 2(n-1)s D^{n-2} b_{s+1}^{jk}.
\end{aligned}
\tag{A10}
\end{equation}
$$

where $D \equiv d/d\alpha$.

## A2 Relationship to the classical one-dimensional Laplace coefficients

It is also useful to relate the new Laplace coefficients to those of the nearly coplanar disturbing function (Brouwer & Clemence 1961). Simple algebra shows that $b_s^{jk}(\alpha, I_r = 0^\circ) = 0$ unless $k = j$ and that $b_s^{jk}(\alpha, I_r = 180^\circ) = 0$ unless $k = -j$. Furthermore,

$$b_s^{jj}(\alpha, 0) = b_s^{j(-j)}(\alpha, 180^\circ) \quad (A11)$$

$$ \begin{aligned} & = \frac{2}{\pi} \int_{0}^{2\pi} \frac{\cos(ju) \, du}{(1 + \alpha^2 - 2\alpha \cos u)^s}, \\ & = 2b_s^j(\alpha) \end{aligned} \quad (A12) $$

where $b_s^j(\alpha)$ is the classical one-dimensional Laplace coefficient. One may also check that the derivative relationships (A9) and (A10) reduce to the known relations of the one-dimensional Laplace coefficients when $I_r = 0^\circ$. The coplanar identity also extends to the functions $A_{i,j,k,l}$ as

$$A_{i,j,k,l}(\alpha, I_r = 0^\circ) = 0 \quad \text{unless } k = j, \quad (A13)$$

$$A_{i,j,k,l}(\alpha, I_r = 180^\circ) = 0 \quad \text{unless } k = -j. \quad (A14)$$

## A3 Series expansion with respect to $\alpha$

The numerical values of a given $b_s^{jk}$ coefficient and its derivatives may be obtained in two ways. One is to integrate numerically the expression (A1) and use the identities (A9) and (A10) for its derivatives. Whereas the current computing power of a personal computer would be sufficient for small $j$ and $k$, the integral is highly oscillatory for larger integers making its evaluation time consuming. The other way is to seek a series expansion of the new two-dimensional Laplace coefficients with respect to the semimajor axis ratio $\alpha$ in the same way the classical one-dimensional Laplace coefficients are evaluated as a hypergeometric series (for details on the subject, we refer the reader to Brouwer & Clemence 1961, p. 495). Series expansions with respect to the semimajor axis ratio suffer from the well-known feature of slow convergence when $\alpha \ll 1$ is not satisfied, thus requiring large expansion orders. None the less, the series method remains superior in terms of computation time to direct numerical integration of the full integral (A1) provided a sufficiently large expansion order is used.

In order to determine the series expansion of the Laplace coefficient with respect to $\alpha$ to order $N_\alpha$, we may use the definition of Gegenbauer polynomials (also known as ultraspherical polynomials) as the coefficients of the series expansion of $(1 + \alpha^2 - 2\alpha x)^{-s}$ with respect to $\alpha$ (Abramowitz & Stegun 1972). This yields

$$ b_s^{jk} = \frac{1}{\pi^2} \sum_{v=0}^{N_\alpha} \alpha^v \int_0^{2\pi} \int_0^{2\pi} C_v^s (\cos u \cos v - \sin u \sin v \cos I) \cos(ju + kv) \, du \, dv, \quad (A15) $$

where $C_v^s(x)$ are Gegenbauer polynomials. They may be expressed from the hypergeometric function ${}_2F_1(-\nu, \nu + 2s, s + 1/2, z) = \nu! C_v^s(1 - 2z)/(2s)_\nu$ as the polynomial

$$C_v^s(x) = \sum_{n=0}^{\text{Floor}[\nu/2]} \frac{(-1)^n (s)_{(\nu-n)} (2x)^{\nu-2n}}{n!(\nu-2n)!} \quad (A16)$$

where $(s)_\nu = s(s+1)(s+2)...(s+\nu-1)$. In the expression (A15), the semimajor axis ratio $\alpha$ must be smaller than unity. The Laplace coefficients for $\alpha > 1$ are found using the relationship (A3). The integers $j$ and $k$ may also be restricted to $0 \le k \le j$ as the Laplace coefficient of any other values may be brought to this case using the relationships (A4) and (A5). The argument of the Gegenbauer polynomial in the series (A15) may be written as

$$x = \cos u \cos v - \sin u \sin v \cos I, \quad (A17)$$

$$= [\eta_1 \cos(u-v) + \eta_2 \cos(u+v)]/2, \quad (A18)$$

where $\eta_1 = 1 - \cos I$ and $\eta_2 = 1 + \cos I$. The Fourier coefficient of the Gegenbauer polynomials requires integrating the powers $x^m$ where $m = \nu - 2n$. Expressing the cosines in $x$ in terms of the usual imaginary exponentials yields

$$ x^m = 2^{-2m} \sum_{l=0}^{m} \binom{m}{l} \eta_1^{l} \eta_2^{m-l} \sum_{p=0}^{l} \sum_{q=0}^{m-l} \binom{l}{p} \binom{m-l}{q} e^{i[2(q+p)-m]u} e^{i[2(q-p)-m+2l]v}. \quad (A19) $$

Multiplying $x^m$ by $\cos(ju + kv) = [e^{i(ju + kv)} + e^{-i(ju + kv)}]/2$, yields the sum of the two terms:

$$ x_{jk+}^m = 2^{-2m-1} \sum_{l=0}^{m} \binom{m}{l} \eta_1^l \eta_2^{m-l} \sum_{p=0}^{l} \sum_{q=0}^{m-l} \binom{l}{p} \binom{m-l}{q} e^{i[2(q+p)-m+j]u} e^{i[2(q-p)-m+2l+k]v}, \quad (A20) $$

$$ x_{jk-}^m = 2^{-2m-1} \sum_{l=0}^{m} \binom{m}{l} \eta_1^l \eta_2^{m-l} \sum_{p=0}^{l} \sum_{q=0}^{m-l} \binom{l}{p} \binom{m-l}{q} e^{i[2(q+p)-m-j]u} e^{i[2(q-p)-m+2l-k]v}. \quad (A21) $$

Integrating $x_{jk+}^m$ and $x_{jk-}^m$ with respect to $u$ and $v$ between 0 and $2\pi$ requires that the powers of the exponentials be zero because otherwise the integrals vanish as they are $2\pi$-periodic. This gives the following conditions on $p$ and $q$:

$$x_{jk+}^m : 2(q+p) = m-j, \quad 2(q-p) = m-k-2l,$$

$$x_{jk-}^m : 2(q+p) = m+j, \quad 2(q-p) = m+k-2l.$$

The first relation implies that $m$ must be of the same parity as $j$ (and therefore $k$) and that $m \ge j$. The integers $p$ and $q$ are determined uniquely as

$$x_{jk+}^m : p = \frac{2l - j + k}{4}, \quad q = \frac{2(m - l) - (j + k)}{4}, \quad (A22)$$

$$x_{jk-}^m : p = \frac{2l + j - k}{4}, \quad q = \frac{2(m - l) + (j + k)}{4}. \quad (A23)$$

The first and third equations require that $l$ have the same parity as $(j-k)/2$. The second and fourth require that $m-l$ have the same parity as $(j+k)/2$ (which is that of $(j-k)/2$ if $k$ and $j$ are even and its opposite if the integers are odd). This last condition is redundant as when combined with the first, it shows that $m$ has the same parity as $j$ and $k$. For the integers $p$ and $q$ to exist, they must satisfy $0 \le p \le l$ and $0 \le q \le m-l$. This translates into the following relations
for $x_{jk+}^m$:

$$p \ge 0 \text{ implies } l \ge \frac{j-k}{2}, \qquad (A24)$$

$$p \le l \text{ implies } l \ge \frac{k-j}{2}, \qquad (A25)$$

trivial condition as $j \ge k$,

$$q \ge 0 \text{ implies } l \le m - \frac{j+k}{2}, \qquad (A26)$$

$$q \le m-l \text{ implies } l \le m + \frac{j+k}{2}, \qquad (A27)$$

trivial condition as $l \le m$.

Conditions (A24) and (A26) determine the range of the $l$-sum as $(j-k)/2 \le l \le m - (j+k)/2$ with $l$ of the same parity as $(j-k)/2$. For the second term $x_{jk-}^m$, the requirements yield exactly the same summation range for the $l$-sum. Calling $y_{jk\pm}^m = \pi^{-2} \int_0^{2\pi} \int_0^{2\pi} x_{jk\pm}^m du dv$, simple algebra shows that $y_{jk+}^m = y_{jk-}^m$ and

$$y_{jk\pm}^m = \frac{m!}{2^{2m-1}} \sum_{\substack{l \le k \le l' \le m - \frac{l+k}{2} \\ l \text{ same parity as } l \ge \frac{l-k}{2}}} \frac{\eta_1^{l} \eta_2^{m-l}}{[(2l - j + k)/4]![(2l - j - k)/4]!} \times \frac{1}{[(2m - 2l - k - j)/4]![(2m - 2l + k + j)/4]!}$$

The summation over $l$ can be simplified by removing the parity condition through the transformation $l \to 2l + (j-k)/2$ that gives a more compact form of $y_{jk\pm}^m$ as

$$y_{jk\pm}^m = \frac{m! \eta_1^{\frac{l-k}{2}} \eta_2^{\frac{k-l}{2}}}{2^{2m-1}} \sum_{l=0} \frac{\eta_1^{2l} \eta_2^{m-2l}}{l![(2l+j-k)/2]!} \times \frac{1}{[(m-2l-j)/2]![(m-2l+k)/2]!} \qquad (A28)$$

The condition that $m$ must have the same parity as $j$ and satisfy $m \ge j$ restricts the Gegenbauer polynomials (A16) to $0 \le n \le \text{Floor}[(\nu - j)/2]$ and the full two-dimensional Laplace coefficient series (A15) to $\nu \ge j$ and $\nu$ of the same parity as $j$. By substituting twice the sum (A28) into the expression of the Gegenbauer polynomials (A16), and applying the index change $\nu \to 2\nu + j$, the series expansion of the two-dimensional Laplace coefficient with arbitrary inclination with respect to the semimajor axis ratio $\alpha$ to

order $N_\alpha$ (A15) may be written as

$$b_s^{jk} = \frac{(1 - \cos I_r)^{\frac{j-k}{2}} (1 + \cos I_r)^{\frac{j+k}{2}} \alpha^j}{2^{j-2}} \times \sum_{\nu=0}^{\text{Floor}[(N_\alpha - j)/2]} \alpha^{2\nu} \sum_{n=0}^{\nu} \frac{(-1)^n (s)_{(2\nu+j-n)}}{n! 2^{2(\nu-n)}} \times \sum_{l=0}^{\nu-n} \frac{(1 - \cos I_r)^{2l} (1 + \cos I_r)^{2(\nu-n-l)}}{l!(l + (j-k)/2)!} \times \frac{1}{(\nu - n - l)!(\nu - n - l + (j+k)/2)!} \qquad (A29)$$

We remind the reader that expression (A29) assumes $\alpha < 1$ and $0 \le k \le j$. For different parameters, the symmetries of the Laplace coefficients (A3), (A4) and (A5) may be used to get back to the last two conditions.

For $I_r = 0^\circ$ and $j=k$ (as required by A11), all terms in the $l$-sum vanish except $l=0$ thus simplifying the general expression (A29) to

$$b_s^{jj}(I_r = 0^\circ) = 4\alpha^j \sum_{\nu=0}^{\text{Floor}[(N_\alpha - j)/2]} \alpha^{2\nu} \sum_{n=0}^{\nu} \frac{(-1)^n (s)_{(2\nu-n+j)}}{n!(\nu-n)!(\nu-n+j)!}, \qquad (A30)$$

which equals twice the classical series of the one-dimensional Laplace coefficient (A12) (Brouwer & Clemence 1961).

The series expansion of the Laplace coefficient (A29) is used in Section 5 to determine resonance widths that utilize $b_s^{jj}$ for $0^\circ \le I_r \le 180^\circ$. Comparison with the direct numerical integration of the expression (A1) shows that for the 2:1 and hence 1:2 resonance locations, an expansion order of $N_\alpha = 20$ is required to reach a $10^{-6}$ relative error whereas for the 3:1 and hence the 1:3 resonance locations, the relative error is $10^{-9}$ with the same expansion order.

## APPENDIX B: FOURTH-ORDER SERIES EXPANSION WITH ARBITRARY INCLINATION

In this appendix, we produce the explicit force coefficients of the disturbing function with arbitrary reference inclination. Tables B1–B5 contain the force coefficients of the cosine terms $\phi_k^{p:q} = \phi - k\omega$, where $\phi = q\lambda - p\lambda' - (q-p)\Omega$. The indirect part of the disturbing function is given in Table B6 and the secular potential in Table B7. The force coefficients use the quantities $A_{i,p,q,l} = \alpha^l D^l b_{i+1/2}^{jk}$ that appear in the literal expansion (equation 12) and whose properties are further discussed in Appendix A.
<table><caption>Table B1. Force coefficients $c_{mn}^0(p, q, \alpha, I_r)$ of the term $e^m s^n \cos \phi$.</caption><tbody><tr><td>c<sup>0</sup><sub>00</sub></td><td><sup>1</sup>&frasl;<sub>2</sub>A<sub>0,p,q,0</sub></td></tr><tr><td>c<sup>0</sup><sub>01</sub></td><td>-<sup>&alpha;</sup>&frasl;<sub>8</sub>(A<sub>1,p-1,q-1,0</sub> - A<sub>1,p-1,q+1,0</sub> - A<sub>1,p+1,q-1,0</sub> + A<sub>1,p+1,q+1,0</sub>)sin I<sub>r</sub></td></tr><tr><td>c<sup>0</sup><sub>20</sub></td><td><sup>1</sup>&frasl;<sub>8</sub>(-4q<sup>2</sup>A<sub>0,p,q,0</sub> + 2A<sub>0,p,q,1</sub> + A<sub>0,p,q,2</sub>)</td></tr><tr><td>c<sup>0</sup><sub>02</sub></td><td><sup>3&alpha;<sup>2</sup></sup>&frasl;<sub>64</sub>(A<sub>2,p-2,q-2,0</sub> - 2A<sub>2,p-2,q,0</sub> + A<sub>2,p-2,q+2,0</sub> - 2A<sub>2,p,q-2,0</sub> + 4A<sub>2,p,q,0</sub> - 2A<sub>2,p,q+2,0</sub> + A<sub>2,p+2,q-2,0</sub> - 2A<sub>2,p+2,q,0</sub> + A<sub>2,p+2,q+2,0</sub>)sin<sup>2</sup>I<sub>r</sub> - <sup>&alpha;</sup>&frasl;<sub>16</sub>(A<sub>1,p-1,q-1,0</sub> - A<sub>1,p-1,q+1,0</sub> - A<sub>1,p+1,q-1,0</sub> + A<sub>1,p+1,q+1,0</sub>)cos I<sub>r</sub></td></tr><tr><td>c<sup>0</sup><sub>21</sub></td><td><sup>&alpha;</sup>&frasl;<sub>32</sub>([4q<sup>2</sup> - 2](A<sub>1,p-1,q-1,0</sub> - A<sub>1,p-1,q+1,0</sub> - A<sub>1,p+1,q-1,0</sub> + A<sub>1,p+1,q+1,0</sub> + A<sub>1,p-1,q+1,1</sub> + A<sub>1,p+1,q-1,1</sub> - A<sub>1,p+1,q+1,1</sub>) - A<sub>1,p-1,q-1,2</sub> + A<sub>1,p-1,q+1,2</sub> + A<sub>1,p+1,q-1,2</sub> - A<sub>1,p+1,q+1,2</sub>]sin I<sub>r</sub></td></tr><tr><td>c<sup>0</sup><sub>03</sub></td><td>-<sup>5&alpha;<sup>3</sup></sup>&frasl;<sub>256</sub>[3(A<sub>3,p-3,q+1,0</sub> - A<sub>3,p-3,q-1,0</sub> - A<sub>3,p-1,q-3,0</sub> + A<sub>3,p-1,q+3,0</sub> + A<sub>3,p-1,q-3,0</sub> + A<sub>3,p+3,q-1,0</sub> - A<sub>3,p+3,q+1,0</sub>) + 9(A<sub>3,p-1,q-1,0</sub> - A<sub>3,p-1,q+1,0</sub> - A<sub>3,p+1,q-1,0</sub> + A<sub>3,p+1,q+1,0</sub>) - A<sub>3,p+3,q-3,0</sub> - A<sub>3,p-3,q+3,0</sub> + A<sub>3,p-3,q-3,0</sub> + A<sub>3,p+3,q+3,0</sub>]sin<sup>3</sup> I<sub>r</sub> + <sup>3&alpha;<sup>2</sup></sup>&frasl;<sub>64</sub>[A<sub>2,p-2,q-2,0</sub> + A<sub>2,p+2,q-2,0</sub> + A<sub>2,p+2,q+2,0</sub> + A<sub>2,p-2,q+2,0</sub> + 2(2A<sub>2,p,q,0</sub> - A<sub>2,p-2,q,0</sub> - A<sub>2,p,q+2,0</sub> - A<sub>2,p+2,q,0</sub>)]cos I<sub>r</sub> sin I<sub>r</sub></td></tr><tr><td>c<sup>0</sup><sub>40</sub></td><td><sup>1</sup>&frasl;<sub>128</sub>[q<sup>2</sup>(16q<sup>2</sup> - 9)A<sub>0,p,q,0</sub> - 8q<sup>2</sup>(A<sub>0,p,q,1</sub> + A<sub>0,p,q,2</sub>) + 4A<sub>0,p,q,3</sub> + A<sub>0,p,q,4</sub>]</td></tr><tr><td>c<sup>0</sup><sub>04</sub></td><td><sup>3&alpha;<sup>2</sup></sup>&frasl;<sub>32768</sub>[64(A<sub>2,p-2,q-2,0</sub> - 2A<sub>2,p-2,q,0</sub> + A<sub>2,p-2,q+2,0</sub> - 2A<sub>2,p,q-2,0</sub> + 4A<sub>2,p,q,0</sub> - 2A<sub>2,p,q+2,0</sub> + A<sub>2,p+2,q-2,0</sub> - 2A<sub>2,p+2,q,0</sub> + A<sub>2,p+2,q+2,0</sub>) + 35&alpha;<sup>2</sup>(A<sub>4,p-4,q-4,0</sub> + A<sub>4,p+4,q-4,0</sub> + A<sub>4,p-4,q+4,0</sub> + A<sub>4,p+4,q+4,0</sub> - 4(A<sub>4,p+2,q-4,0</sub> + A<sub>4,p+4,q+2,0</sub> + A<sub>4,p+2,q+4,0</sub> + A<sub>4,p+4,q-2,0</sub> + A<sub>4,p-4,q+2,0</sub> + A<sub>4,p-4,q-4,0</sub>) + 6(A<sub>4,p-4,q,0</sub> + A<sub>4,p,q-4,0</sub> + A<sub>4,p,q+4,0</sub> + A<sub>4,p+4,q,0</sub>) + 16(A<sub>4,p-2,q-2,0</sub> + A<sub>4,p-2,q+2,0</sub> + A<sub>4,p+2,q-2,0</sub> + A<sub>4,p+2,q+2,0</sub>) - 24(A<sub>4,p,q-2,0</sub> + A<sub>4,p-2,q,0</sub> + A<sub>4,p,q+2,0</sub> + A<sub>4,p+2,q,0</sub>) + 36A<sub>4,p,q,0</sub>] + <sup>&alpha;</sup>&frasl;<sub>2048</sub>[32(A<sub>1,p+1,q-1,0</sub> - A<sub>1,p-1,q-1,0</sub> + A<sub>1,p-1,q+1,0</sub> - A<sub>1,p+1,q+1,0</sub>) + 15&alpha;<sup>2</sup>(A<sub>3,p+3,q-3,0</sub> - A<sub>3,p-3,q-3,0</sub> - A<sub>3,p+3,q+3,0</sub> + A<sub>3,p-3,q+3,0</sub> + 3(A<sub>3,p-3,q-1,0</sub> - A<sub>3,p-3,q+1,0</sub> + A<sub>3,p-1,q-3,0</sub> - A<sub>3,p-1,q+3,0</sub> - A<sub>3,p+1,q-3,0</sub> + A<sub>3,p+1,q+3,0</sub>) + 9(A<sub>3,p+1,q-1,0</sub> - A<sub>3,p+1,q+1,0</sub> - A<sub>3,p-1,q-1,0</sub> + A<sub>3,p-1,q+1,0</sub>)]cos I<sub>r</sub></td></tr></tbody></table>

<table><caption>Table B1 – continued</caption><tbody><tr><td>c<sup>0</sup><sub>22</sub></td><td>+<sup>&alpha;<sup>2</sup></sup>&frasl;<sub>8192</sub>[48(A<sub>2,p-2,q-2,0</sub> - 2A<sub>2,p-2,q,0</sub> + A<sub>2,p-2,q+2,0</sub> - 2A<sub>2,p,q-2,0</sub> + 4A<sub>2,p,q,0</sub> - 2A<sub>2,p,q+2,0</sub> + A<sub>2,p+2,q-2,0</sub> - 2A<sub>2,p+2,q,0</sub> + A<sub>2,p+2,q+2,0</sub>) + 35&alpha;<sup>2</sup>(4(A<sub>4,p-4,q-2,0</sub> + A<sub>4,p-2,q+4,0</sub> + A<sub>4,p-4,q+4,0</sub> + A<sub>4,p-4,q-4,0</sub> + A<sub>4,p+2,q-4,0</sub> + A<sub>4,p+2,q+4,0</sub> + A<sub>4,p+4,q+2,0</sub>) - A<sub>4,p-4,q-4,0</sub> - A<sub>4,p+4,q-4,0</sub> - 6(A<sub>4,p,q-4,0</sub> + A<sub>4,p-4,q,0</sub> + A<sub>4,p,q+4,0</sub> + A<sub>4,p+4,q,0</sub>) - 16(A<sub>4,p+2,q-2,0</sub> + A<sub>4,p+2,q+2,0</sub> + A<sub>4,p-2,q-2,0</sub> + A<sub>4,p-2,q+2,0</sub>) + 24(A<sub>4,p,q-2,0</sub> + A<sub>4,p,q+2,0</sub> + A<sub>4,p+2,q,0</sub> + A<sub>4,p-2,q,0</sub>) - 36A<sub>4,p,q,0</sub>)cos 2I<sub>r</sub> + <sup>15&alpha;<sup>3</sup></sup>&frasl;<sub>2048</sub>[A<sub>3,p-3,q-3,0</sub> - A<sub>3,p+3,q-3,0</sub> - A<sub>3,p-3,q+3,0</sub> + A<sub>3,p+3,q+3,0</sub> + 3(A<sub>3,p-3,q+1,0</sub> - A<sub>3,p-3,q-1,0</sub> - A<sub>3,p-1,q-3,0</sub> + A<sub>3,p-1,q+3,0</sub> + A<sub>3,p+1,q-3,0</sub> - A<sub>3,p+1,q+3,0</sub> + A<sub>3,p+1,q-1,0</sub> - A<sub>3,p-1,q-3,0</sub> + A<sub>3,p-1,q+3,0</sub>) + 9(A<sub>3,p+1,q+1,0</sub> - A<sub>3,p+1,q-1,0</sub> + A<sub>3,p-1,q-1,0</sub> - A<sub>3,p-1,q+1,0</sub>)cos 3I<sub>r</sub> + <sup>35&alpha;<sup>4</sup></sup>&frasl;<sub>32768</sub>(A<sub>4,p-4,q-4,0</sub> + A<sub>4,p-4,q+4,0</sub> + A<sub>4,p+4,q+4,0</sub> + A<sub>4,p+4,q-4,0</sub> - 4(A<sub>4,p-4,q-2,0</sub> + A<sub>4,p-4,q+2,0</sub> + A<sub>4,p-4,q-2,0</sub> + A<sub>4,p+4,q+2,0</sub> + A<sub>4,p-4,q-4,0</sub> + A<sub>4,p+2,q-4,0</sub>) + 6(A<sub>4,p,q-4,0</sub> + A<sub>4,p-4,q,0</sub> + A<sub>4,p,q+4,0</sub> + A<sub>4,p+4,q,0</sub>) + 16(A<sub>4,p-2,q-2,0</sub> + A<sub>4,p-2,q+2,0</sub> + A<sub>4,p+2,q+2,0</sub> + A<sub>4,p+2,q-2,0</sub>) - 24(A<sub>4,p+2,q,0</sub> + A<sub>4,p-2,q,0</sub> + A<sub>4,p,q-2,0</sub> + A<sub>4,p,q+2,0</sub>) + 36A<sub>4,p,q,0</sub>)cos 4I<sub>r</sub> - <sup>3&alpha;<sup>2</sup></sup>&frasl;<sub>256</sub>[2(2q<sup>2</sup> - 3)(A<sub>2,p-2,q-2,0</sub> - 2A<sub>2,p-2,q,0</sub> + A<sub>2,p-2,q+2,0</sub> - 2A<sub>2,p,q-2,0</sub> + 4A<sub>2,p,q,0</sub> - 2A<sub>2,p+2,q,0</sub> + A<sub>2,p+2,q-2,0</sub> + A<sub>2,p+2,q+2,0</sub>) - 6A<sub>2,p-2,q-1</sub> - A<sub>2,p-2,q-2,2</sub> + 12A<sub>2,p-2,q,1</sub> + 2A<sub>2,p-2,q,2</sub> - 6A<sub>2,p-2,q+2,1</sub> - A<sub>2,p-2,q+2,2</sub> + 12A<sub>2,p,q-2,1</sub> + 2A<sub>2,p,q-2,2</sub> - 24A<sub>2,p,q,1</sub> - 4A<sub>2,p,q,2</sub> + 12A<sub>2,p,q+2,1</sub> + 2A<sub>2,p,q+2,2</sub> - 6A<sub>2,p+2,q-2,1</sub> - A<sub>2,p+2,q-2,2</sub> + 12A<sub>2,p+2,q,1</sub> + 2A<sub>2,p+2,q,2</sub> - 6A<sub>2,p+2,q+2,1</sub> - A<sub>2,p+2,q+2,2</sub>] sin<sup>2</sup> I<sub>r</sub> + <sup>&alpha;</sup>&frasl;<sub>64</sub>[2(2q<sup>2</sup> - 1)(A<sub>1,p-1,q-1,0</sub> - A<sub>1,p-1,q+1,0</sub> - A<sub>1,p+1,q-1,0</sub> + A<sub>1,p+1,q+1,0</sub>) + 4(A<sub>1,p-1,q+1,1</sub> - A<sub>1,p-1,q-1,1</sub> + A<sub>1,p+1,q-1,1</sub> - A<sub>1,p+1,q+1,1</sub>) + A<sub>1,p+1,q-1,2</sub> + A<sub>1,p-1,q+1,2</sub> - A<sub>1,p-1,q-1,2</sub> + A<sub>1,p+1,q+1,2</sub>]cos I<sub>r</sub></td></tr></tbody></table>
<table><caption>Table B2. Force coefficients $c_{mn}^1(p, q, \alpha, I_r)$ of the term $e^m s^n \cos(\phi - \omega)$.</caption>
  <tbody>
    <tr>
      <td>c<sup>1</sup><sub>10</sub></td>
      <td>-1/4[2(1 - q)A<sub>0,p,q-1,0</sub> + A<sub>0,p,q-1,1</sub>]</td>
    </tr>
    <tr>
      <td>c<sup>1</sup><sub>11</sub></td>
      <td>-α/16[(2q - 3)(A<sub>1,p-1,q-2,0</sub> - A<sub>1,p-1,q,0</sub> + A<sub>1,p+1,q,0</sub> - A<sub>1,p+1,q-2,0</sub>) - A<sub>1,p-1,q-2,1</sub> + A<sub>1,p-1,q,1</sub> + A<sub>1,p+1,q-2,1</sub> - A<sub>1,p+1,q,1</sub>]sin I<sub>r</sub></td>
    </tr>
    <tr>
      <td>c<sup>1</sup><sub>12</sub></td>
      <td>(3α<sup>2</sup>/128)[2(2 - q)(A<sub>2,p-2,q-3,0</sub> - 2A<sub>2,p-2,q-1,0</sub> + A<sub>2,p-2,q+1,0</sub> - 2A<sub>2,p,q-3,0</sub> + 4A<sub>2,p,q-1,0</sub> + A<sub>2,p+2,q-3,0</sub> - 2A<sub>2,p,q+1,0</sub> - 2A<sub>2,p+2,q-1,0</sub> + A<sub>2,p+2,q+1,0</sub>) + A<sub>2,p-2,q-3,1</sub> - 2A<sub>2,p-2,q-1,1</sub> + A<sub>2,p-2,q+1,1</sub> - 2A<sub>2,p,q-3,1</sub> + 4A<sub>2,p,q-1,1</sub> - 2A<sub>2,p,q+1,1</sub> + A<sub>2,p+2,q-3,1</sub> - 2A<sub>2,p+2,q-1,1</sub> + A<sub>2,p+2,q+1,1</sub>] sin<sup>2</sup> I<sub>r</sub> + α/32 [(3 - 2q)(A<sub>1,p-1,q-2,0</sub> - A<sub>1,p-1,q,0</sub> - A<sub>1,p+1,q-2,0</sub> + A<sub>1,p+1,q,0</sub>) + A<sub>1,p-1,q-2,1</sub> - A<sub>1,p-1,q,1</sub> - A<sub>1,p+1,q-2,1</sub> + A<sub>1,p+1,q,1</sub>cos I<sub>r</sub>]</td>
    </tr>
    <tr>
      <td>c<sup>1</sup><sub>30</sub></td>
      <td>1/32[2(7q - 4q<sup>2</sup> - 3)A<sub>0,p,q-1,0</sub> + q(4q - 1)A<sub>0,p,q-1,1</sub> + (2q - 4)A<sub>0,p,q-1,2</sub> - A<sub>0,p,q-1,3</sub>]</td>
    </tr>
    <tr>
      <td>c<sup>1</sup><sub>13</sub></td>
      <td>(5α<sup>3</sup>/2048)[(2q - 5)(A<sub>3,p-3,q-4,0</sub> - 3A<sub>3,p-3,q-2,0</sub> + 3A<sub>3,p-3,q,0</sub> - A<sub>3,p-3,q+2,0</sub> - 3A<sub>3,p-1,q-4,0</sub> + 9A<sub>3,p-1,q-2,0</sub> - 9A<sub>3,p-1,q,0</sub> + 3A<sub>3,p-1,q+2,0</sub> + 3A<sub>3,p+1,q-4,0</sub> - 9A<sub>3,p+1,q-2,0</sub> + 9A<sub>3,p+1,q,0</sub> - 3A<sub>3,p+1,q+2,0</sub> - A<sub>3,p+3,q-4,0</sub> - 3A<sub>3,p+3,q,0</sub> + 3A<sub>3,p+3,q-2,0</sub> + A<sub>3,p+3,q+2,0</sub>) + 3(A<sub>3,p-3,q-2,1</sub> - A<sub>3,p-3,q,1</sub> + A<sub>3,p-1,q-4,1</sub> - A<sub>3,p-1,q+2,1</sub> - A<sub>3,p+1,q-4,1</sub> + A<sub>3,p+1,q+2,1</sub> - A<sub>3,p+3,q-2,1</sub> + A<sub>3,p+3,q,1</sub>) + 9(A<sub>3,p+1,q-2,0</sub> - A<sub>3,p-1,q-2,1</sub> + A<sub>3,p-1,q,1</sub> - A<sub>3,p+1,q,1</sub>) + A<sub>3,p+3,q-4,1</sub> + A<sub>3,p-3,q+2,1</sub> - A<sub>3,p+3,q+2,1</sub> - A<sub>3,p-3,q-4,1</sub>] × (sin 3I<sub>r</sub> - 3sin I<sub>r</sub>) + 3α<sup>2</sup>/256 [2(q - 2)(A<sub>2,p-2,q-3,0</sub> - 2A<sub>2,p-2,q-1,0</sub> + A<sub>2,p-2,q+1,0</sub> - 2A<sub>2,p,q-3,0</sub> + 4A<sub>2,p,q-1,0</sub> - 2A<sub>2,p,q+1,0</sub> + A<sub>2,p+2,q-3,0</sub> - 2A<sub>2,p+2,q-1,0</sub> + A<sub>2,p+2,q+1,0</sub>) - A<sub>2,p-2,q-3,1</sub> - A<sub>2,p+2,q+1,1</sub> - A<sub>2,p-2,q+1,1</sub> - A<sub>2,p+2,q-3,1</sub> + 2(A<sub>2,p-2,q-1,1</sub> + A<sub>2,p,q-3,1</sub> + A<sub>2,p,q+1,1</sub> + A<sub>2,p+2,q-1,1</sub>) - 4A<sub>2,p,q-1,1</sub>sin 2I<sub>r</sub>]</td>
    </tr>
    <tr>
      <td>c<sup>1</sup><sub>31</sub></td>
      <td>α/128 [q(7 - 18q + 8q<sup>2</sup>)(A<sub>1,p-1,q-2,0</sub> - A<sub>1,p-1,q,0</sub> - A<sub>1,p+1,q-2,0</sub> + A<sub>1,p+1,q,0</sub>) + (8 - 3q - 4q<sup>2</sup>)(A<sub>1,p-1,q-2,1</sub> - A<sub>1,p-1,q,1</sub> - A<sub>1,p+1,q-2,1</sub> + A<sub>1,p+1,q,1</sub>) + (7 - 2q)(A<sub>1,p-1,q-2,2</sub> - A<sub>1,p-1,q,2</sub> - A<sub>1,p+1,q-2,2</sub> + A<sub>1,p+1,q,2</sub>) + A<sub>1,p-1,q-2,3</sub> - A<sub>1,p-1,q,3</sub> - A<sub>1,p+1,q-2,3</sub> + A<sub>1,p+1,q,3</sub>sin I<sub>r</sub>]</td>
    </tr>
  </tbody>
</table>

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

<table><caption>Table B4. Force coefficients $c_{mn}^3(p, q, \alpha, I_r)$ of the term $e^m s^n \cos(\phi - 3\omega)$.</caption>
  <tbody>
    <tr>
      <td>c<sup>3</sup><sub>30</sub></td>
      <td>1/96 [(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>-α/384 [(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>sin I<sub>r</sub>]</td>
    </tr>
  </tbody>
</table>
<table><caption>**Table B5.** Force coefficients $c_{mn}^4(p, q, \alpha, I_r)$ 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>

<table><caption>**Table B6.** 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>8</sub>[2 s sin I<sub>r</sub> + (s<sup>2</sup> - 2) cos I<sub>r</sub> - 2]e</td></tr><tr><td>λ - λ'</td><td><sup>α</sup>&frasl;<sub>128</sub>[e<sup>4</sup> + 32e<sup>2</sup> - 64 + [e<sup>4</sup> - 16e<sup>2</sup>(s<sup>2</sup> - 2) + 8(s<sup>4</sup> + 4s<sup>2</sup> - 8)]cos I<sub>r</sub> - 32(e<sup>2</sup> - 2)s sin I<sub>r</sub>]</td></tr><tr><td>λ - λ' - 2ω + 2Ω</td><td><sup>α</sup>&frasl;<sub>96</sub>e<sup>2</sup>[(6 + 2e<sup>2</sup> - 3s<sup>2</sup>) cos I<sub>r</sub> - 2(3 + e<sup>2</sup> + 3s sin I<sub>r</sub>)]</td></tr><tr><td>2λ - λ' - ω</td><td><sup>α</sup>&frasl;<sub>16</sub>e[(3e<sup>2</sup> + 2s<sup>2</sup> - 4) cos I<sub>r</sub> - (3e<sup>2</sup> - 4)(s sin I<sub>r</sub> - 1)]</td></tr><tr><td>2λ - λ' - 3ω + 2Ω</td><td><sup>α</sup>&frasl;<sub>48</sub>e<sup>3</sup>(cos I<sub>r</sub> - s sin I<sub>r</sub> - 1)</td></tr><tr><td>3λ - λ' - 2ω</td><td><sup>3α</sup>&frasl;<sub>32</sub>e<sup>2</sup>[(+2e<sup>2</sup> + s<sup>2</sup> - 2) cos I<sub>r</sub> + 2(e<sup>2</sup> + s sin I<sub>r</sub> - 1)]</td></tr><tr><td>3λ - λ' - 4ω + 2Ω</td><td>-<sup>3α</sup>&frasl;<sub>128</sub>e<sup>4</sup> sin<sup>2</sup>(I<sub>r</sub>/2)</td></tr><tr><td>4λ - λ' - 3ω</td><td>-<sup>α</sup>&frasl;<sub>6</sub>e<sup>3</sup>(1 + cos I<sub>r</sub> - s sin I<sub>r</sub>)</td></tr><tr><td>5λ - λ' - 4ω</td><td>-<sup>125α</sup>&frasl;<sub>384</sub>e<sup>4</sup> cos<sup>2</sup>(I<sub>r</sub>/2)</td></tr><tr><td>λ' + ω - 2Ω</td><td><sup>3α</sup>&frasl;<sub>8</sub>e[2 + (s<sup>2</sup> - 2) cos I<sub>r</sub> + 2 s sin I<sub>r</sub>]</td></tr><tr><td>λ + λ' - 2Ω</td><td>-<sup>α</sup>&frasl;<sub>128</sub>[64 - 32e<sup>2</sup> - e<sup>4</sup> + [e<sup>4</sup> - 16e<sup>2</sup>(s<sup>2</sup> - 2) + 8(4s<sup>2</sup> + s<sup>4</sup> - 8)]cos I<sub>r</sub> - 32(e<sup>2</sup> - 2)s sin I<sub>r</sub>]</td></tr><tr><td>λ + λ' - 2ω</td><td>-<sup>α</sup>&frasl;<sub>96</sub>e<sup>2</sup>[(6 + 2e<sup>2</sup> - 3s<sup>2</sup>) cos I<sub>r</sub> + 2(3 + e<sup>2</sup> - 3s sin I<sub>r</sub>)]</td></tr><tr><td>2λ + λ' - 3ω</td><td>-<sup>α</sup>&frasl;<sub>48</sub>e<sup>3</sup>(1 + cos I<sub>r</sub> - s sin I<sub>r</sub>)</td></tr><tr><td>2λ + λ' - ω - 2Ω</td><td><sup>α</sup>&frasl;<sub>16</sub>e[(4 - 3e<sup>2</sup> - 2s<sup>2</sup>) cos I<sub>r</sub> + (3e<sup>2</sup> - 4)(1 + s sin I<sub>r</sub>)]</td></tr><tr><td>3λ + λ' - 4ω</td><td>-<sup>3α</sup>&frasl;<sub>128</sub>e<sup>4</sup>α cos<sup>2</sup>(I<sub>r</sub>/2)</td></tr><tr><td>3λ + λ' - 2ω - 2Ω</td><td>-<sup>3α</sup>&frasl;<sub>32</sub>e<sup>2</sup>[2 - 2e<sup>2</sup> + (2e<sup>2</sup> + s<sup>2</sup> - 2) cos I<sub>r</sub> + 2 s sin I<sub>r</sub>]</td></tr><tr><td>4λ + λ' - 3ω - 2Ω</td><td><sup>α</sup>&frasl;<sub>6</sub>e<sup>3</sup>(cos I<sub>r</sub> - s sin I<sub>r</sub> - 1)</td></tr><tr><td>5λ + λ' - 4ω - 2Ω</td><td>-<sup>125α</sup>&frasl;<sub>384</sub>e<sup>4</sup>α sin<sup>2</sup>(I<sub>r</sub>/2)</td></tr></tbody></table>

<table><caption>**Table B7.** Force coefficients $s_{mn}^k(\alpha, I_r)$ of the secular term $e^m s^n \cos(k\omega)$.</caption><tbody><tr><td>k = 0</td><td><sup>α</sup>&frasl;<sub>8</sub>(A<sub>1,1,-1,0</sub> - A<sub>1,1,1,0</sub>) sin I<sub>r</sub></td></tr><tr><td>s<sup>0</sup><sub>01</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>20</sub></td><td><sup>α</sup>&frasl;<sub>64</sub>[4(A<sub>1,1,-1,0</sub> - A<sub>1,1,1,0</sub>) cos I<sub>r</sub> + 3α(2A<sub>2,0,0,0</sub> + A<sub>2,2,-2,0</sub> - 4A<sub>2,2,0,0</sub> + A<sub>2,2,2,0</sub>)sin<sup>2</sup> I<sub>r</sub>]</td></tr><tr><td>s<sup>0</sup><sub>02</sub></td><td><sup>α<sup>2</sup></sup>&frasl;<sub>256</sub> sin I<sub>r</sub>[12(2A<sub>2,0,0,0</sub> + A<sub>2,2,-2,0</sub> - 4A<sub>2,2,0,0</sub> + A<sub>2,2,2,0</sub>)cos I<sub>r</sub> + 5α(9A<sub>3,1,-1,0</sub> - 9A<sub>3,1,1,0</sub> + A<sub>3,3,-3,0</sub> - 6A<sub>3,3,-1,0</sub> + 6A<sub>3,3,1,0</sub> - A<sub>3,3,3,0</sub>)sin<sup>2</sup> I<sub>r</sub>]</td></tr><tr><td>s<sup>0</sup><sub>21</sub></td><td><sup>α</sup>&frasl;<sub>32</sub>(2A<sub>1,1,-1,0</sub> + 4A<sub>1,1,-1,1</sub> + A<sub>1,1,-1,2</sub> - 2A<sub>1,1,1,0</sub> - 4A<sub>1,1,1</sub> - A<sub>1,1,1,2</sub>)sin I<sub>r</sub></td></tr><tr><td>s<sup>0</sup><sub>22</sub></td><td><sup>α</sup>&frasl;<sub>512</sub>[3α(12A<sub>2,0,0,0</sub> + 12A<sub>2,0,0,1</sub> + 2A<sub>2,0,0,2</sub> + 6A<sub>2,2,-2,0</sub> + 6A<sub>2,2,-2,1</sub> + A<sub>2,2,-2,2</sub> - 24A<sub>2,2,0,0</sub> - 24A<sub>2,2,0,1</sub> - 4A<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>) - 3α(12A<sub>2,0,0,0</sub> + 12A<sub>2,0,0,1</sub> + 2A<sub>2,0,0,2</sub> + 6A<sub>2,2,-2,0</sub> + 6A<sub>2,2,-2,1</sub> + A<sub>2,2,-2,2</sub> - 24A<sub>2,2,0,0</sub> - 24A<sub>2,2,0,1</sub> - 4A<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>)cos(2I<sub>r</sub>) + 8(2A<sub>1,1,-1,0</sub> + 4A<sub>1,1,-1,1</sub> + A<sub>1,1,-1,2</sub> - 2A<sub>1,1,1,0</sub> - 4A<sub>1,1,1,1</sub> - A<sub>1,1,1,2</sub>)sin I<sub>r</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>s<sup>0</sup><sub>04</sub></td><td><sup>3α<sup>2</sup></sup>&frasl;<sub>32768</sub>[64(2A<sub>2,0,0,0</sub> + A<sub>2,2,-2,0</sub> - 4A<sub>2,2,0,0</sub> + A<sub>2,2,2,0</sub>) + 7α<sup>2</sup>[90A<sub>4,0,0,0</sub> + 40(2A<sub>4,2,-2,0</sub> + 2A<sub>4,2,2,0</sub> - A<sub>4,4,-2,0</sub> - A<sub>4,4,2,0</sub>) - 60(4A<sub>4,2,0,0</sub> + A<sub>4,4,0,0</sub>) + 5(A<sub>4,4,4,0</sub> + A<sub>4,4,-4,0</sub>)] + <sup>α</sup>&frasl;<sub>2048</sub>[32(A<sub>1,1,-1,0</sub> - A<sub>1,1,1,0</sub>) + 15α<sup>2</sup>(9A<sub>3,1,-1,0</sub> - 9A<sub>3,1,1,0</sub> - 6A<sub>3,3,-1,0</sub> + 6A<sub>3,3,1,0</sub> - A<sub>3,3,3,0</sub>) + A<sub>3,3,-3,0</sub>)]cos I<sub>r</sub> - <sup>α<sup>2</sup></sup>&frasl;<sub>8192</sub>[48(4A<sub>2,2,0,0</sub> - 2A<sub>2,0,0,0</sub> - A<sub>2,2,-2,0</sub> - A<sub>2,2,2,0</sub>) + 35α<sup>2</sup>[A<sub>4,4,-4,0</sub> + A<sub>4,4,4,0</sub> - 8A<sub>4,4,-2,0</sub> - 8A<sub>4,4,2,0</sub> + 12A<sub>4,4,0,0</sub> + 16A<sub>4,2,-2,0</sub> + 16A<sub>4,2,2,0</sub> + 18A<sub>4,0,0,0</sub> - 48A<sub>4,2,0,0</sub>]]cos 2I<sub>r</sub> - <sup>15α<sup>3</sup></sup>&frasl;<sub>2048</sub>(9A<sub>3,1,-1,0</sub> - 9A<sub>3,1,1,0</sub> + A<sub>3,3,-3,0</sub> - 6A<sub>3,3,-1,0</sub> + 6A<sub>3,3,1,0</sub> - A<sub>3,3,3,0</sub>)cos 3I<sub>r</sub> + <sup>35α<sup>4</sup></sup>&frasl;<sub>32768</sub>(A<sub>4,4,4,0</sub> + A<sub>4,4,-4,0</sub> - 8A<sub>4,4,-2,0</sub> - 8A<sub>4,4,2,0</sub> + 12A<sub>4,4,0,0</sub> + 16A<sub>4,2,-2,0</sub> + 16A<sub>4,2,2,0</sub> + 18A<sub>4,0,0,0</sub> - 48A<sub>4,2,0,0</sub>)cos 4I<sub>r</sub>]</td></tr></tbody></table>
<table><caption>**Table B7 – continued**</caption><tbody><tr><td>k = 2</td><td>s<sub>20</sub><sup>2</sup></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>s<sub>21</sub><sup>2</sup></td><td>-<sup>&alpha;</sup>&frasl;<sub>64</sub>[12(A<sub>1,1,-1,0</sub> - A<sub>1,1,1,0</sub> - A<sub>1,3,-1,0</sub> + A<sub>1,3,1,0</sub>) + 8(A<sub>1,1,-1,1</sub> - A<sub>1,1,1</sub> - A<sub>1,3,-1,1</sub> + A<sub>1,3,1,1</sub>) + A<sub>1,1,-1,2</sub> - A<sub>1,1,1,2</sub> - A<sub>1,3,-1,2</sub> + A<sub>1,3,1,2</sub>sin I<sub>r</sub>]</td></tr><tr><td>s<sub>22</sub><sup>2</sup></td><td>-<sup>&alpha;</sup>&frasl;<sub>512</sub>[3&alpha;(10[4A<sub>2,2,2,0</sub> + 2A<sub>2,2,2,1</sub> + 4A<sub>2,0,0,0</sub> + 2A<sub>2,0,0,1</sub> + 4A<sub>2,2,-2,0</sub> + 2A<sub>2,2,-2,1</sub> + 4A<sub>2,4,0,0</sub> + 2A<sub>2,4,0,1</sub> - 2A<sub>2,4,-2,0</sub> - A<sub>2,4,-2,1</sub> - 2A<sub>2,4,2,0</sub> - A<sub>2,4,2,1</sub> - 12A<sub>2,2,0,0</sub> - 6A<sub>2,2,0,1</sub>] - 6A<sub>2,2,0,2</sub> + 2[A<sub>2,0,0,2</sub> + A<sub>2,2,-2,2</sub> + A<sub>2,2,2,2</sub> + A<sub>2,4,0,2</sub>] - A<sub>2,4,2,2</sub> - A<sub>2,4,-2,2</sub>sin<sup>2</sup>I<sub>r</sub> + 4(12[A<sub>1,1,-1,0</sub> - A<sub>1,1,0</sub> - A<sub>1,3,-1,0</sub> + A<sub>1,3,1,0</sub>] + 8[A<sub>1,1,-1,1</sub> - A<sub>1,1,1,1</sub> + A<sub>1,3,1,1</sub> - A<sub>1,3,-1,1</sub>] + A<sub>1,1,-1,2</sub> - A<sub>1,1,1,2</sub> - A<sub>1,3,-1,2</sub> + A<sub>1,3,1,2</sub>cos I<sub>r</sub>]</td></tr><tr><td>k = 4</td><td>s<sub>40</sub><sup>2</sup></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<sub>40</sub><sup>4</sup></td><td></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>

This paper has been typeset from a TeX/LA TeX file prepared by the author.