Vol.0 (2002x) No.0, 000–000

11institutetext: School of Astronomy and Space Science, Nanjing University,
Nanjing 210023, China<EMAIL_ADDRESS>
22institutetext: Key Laboratory of Modern Astronomy and Astrophysics in
Ministry of Education, Nanjing University, Nanjing 210023, China

# A new expansion of planetary disturbing function and applications to
interior, co-orbital and exterior resonances with planets

Hanlun Lei

(Received 2021 month day; accepted 2021 month day)

###### Abstract

In this study, a new expansion of planetary disturbing function is developed
for describing the resonant dynamics of minor bodies with arbitrary
inclinations and semimajor axis ratios. In practice, the disturbing function
is expanded around circular orbits in the first step and then, in the second
step, the resulting mutual interaction between circular orbits is expanded
around a reference point. As usual, the resulting expansion is presented in
the Fourier series form, where the force amplitudes are dependent on the
semimajor axis, eccentricity and inclination and the harmonic arguments are
linear combinations of the mean longitude, longitude of pericenter and
longitude of ascending node of each mass. The resulting new expansion is valid
for arbitrary inclinations and semimajor axis ratios. In the case of mean
motion resonant configuration, the disturbing function can be easily averaged
to produce the analytical expansion of resonant disturbing function. Based on
the analytical expansion, the Hamiltonian model of mean motion resonances is
formulated, and the resulting analytical developments are applied to Jupiter’s
inner and co-orbital resonances and Neptune’s exterior resonances. Analytical
expansion is validated by comparing the analytical results with the associated
numerical outcomes.

###### keywords:

Celestial Mechanics.

##  1 Introduction

Expanding the planetary disturbing function as a power series in orbital
elements is an essential and fundamental issue in the development of
perturbation theories in celestial mechanics. To understand the dynamics of
resonances and secular evolutions, it is a key step to study the properties of
the planetary disturbing function experienced by the object of interest
(Murray & Dermott, 1999). In history, there are various types of expansions
for planetary disturbing function: (a) Laplace-type expansions where the force
amplitudes are expressed in terms of Laplace coefficients (Brouwer & Clemence,
1961; Ellis & Murray, 2000; Murray & Dermott, 1999; Morais & Namouni, 2013;
Namouni & Morais, 2017b, a), (b) Legendre-type expansions where Legendre
polynomials are used in the derivation (Kaula, 1962), (c) asymmetric and local
expansions where the disturbing function is expanded around a specific
resonant center (Ferraz-Mello, 1987; Ferraz-Mello & Sato, 1989; Yokoyama,
1994; Roig et al., 1998), (d) Taylor-type expansions where the disturbing
function is expanded around an artificially chosen reference point in Taylor
series (Páez & Locatelli, 2015; Morais, 1999; Petrovskaya, 1970, 1972; Beaugé,
1996), and (e) semianalytical expansions where the force amplitudes are
determined by numerical Fourier analysis (Klioner, 2000). Among these
expansions, the Laplace- and Legendre-type expansions are classical and they
are widely used in previous studies.

The Laplace-type expansion of planetary disturbing function starts from Peirce
(1849), who derived a sixth-order expansion. As an extension, Newcomb (1895)
performed a seventh-order expansion. Brown & Shook (1933) provided a clear and
detailed derivation for the expansion up to the second-order in terms of the
classical orbital elements. Brouwer & Clemence (1961) formulated a third-order
expansion for the planetary disturbing function, which has been widely adopted
as a standard and low-order expansion. In more recent times, higher-order
explicit expansions of the planetary disturbing function can be found in Ellis
& Murray (2000) and Murray & Dermott (1999), who expanded the disturbing
function as a power series in eccentricities and inclinations of the objects
involved. Particularly, Ellis & Murray (2000) presented a new algorithm for
the determination of the terms associated with a certain argument in the
expansion of planetary disturbing function. Recently, Morais & Namouni (2013)
and Namouni & Morais (2017b, a) made further extensions of the classical
expansions to retrograde orbits, polar orbits and orbits with arbitrary
inclinations. Due to the choice of the reference orbit as a coplanar circular
orbit, the classical expansions are valid in low-inclination and low-
eccentricity configurations. However, the Laplace coefficients and their
derivatives with respect to the semimajor axis ratio arising in the Laplace-
type expansions are divergent when the semimajor axis ratio is approaching
unity. This results in the fact that the Laplace-type expansions are divergent
for co-orbital configurations.

Concerning the Legendre-type expansion, Kaula (1962) expanded the disturbing
function as a power series of the semimajor axis ratio ($\alpha$) between two
objects in triple systems. Please refer to Murray & Dermott (1999) for
explicit expansions of such an expansion. Obviously, the Legendre-type
expansion is applicable for objects with arbitrary inclinations but with
semimajor axis ratio much smaller than unity. It has been widely applied to
varieties of triple systems with hierarchical configurations (Kozai, 1962;
Beaugé et al., 2006; Naoz et al., 2011, 2013; Naoz, 2016; Lei et al., 2018;
Lei, 2019a). However, when the configuration of system is not highly
hierarchical, the disturbing function needs to be truncated at a high enough
order in the semimajor axis ratio in order to achieve a certain precision and,
in this case, the convergence becomes very slow. In particular, when the
semimajor axis ratio is close to unity, the Legendre-type expansion fails to
converge.

According to the aforementioned discussions, both the Laplace- and Legendre-
type expansions of planetary disturbing function are valid for the resonant
dynamics in these configurations where the semimajor axis ratio is not close
to unity. Concerning this issue, we may ask: is there an expansion valid for
minor bodies with arbitrary semimajor axis ratios? To this end, a new
expansion is developed in this study. Based on such a ‘general’ expansion, it
becomes possible to analytically explore dynamical structures of minor bodies
located inside the interior, co-orbital and exterior mean motion resonances at
arbitrary inclinations.

The structure of the remaining part is organized as follows. In Section 2, the
disturbing function for the minor bodies in planetary triple systems is
briefly introduced, and the new expansion of planetary disturbing function is
developed in Section 3. In Section 4, the resonant disturbing function are
presented under the assumption that the perturber is moving in a circular
orbit and then the analytical results obtained by the new expansion is
compared with numerical integration results. The Hamiltonian dynamics are
discussed in Section 5 and the practical applications of the analytical
developments to the interior, co-orbital and exterior resonances are reported
in Section 6. At last, conclusions are summarized in Section 7.

##  2 Planetary disturbing function

The planetary system considered in this study is composed of a central body
with mass $m_{0}$ (i.e. a central star), a planet with mass $m_{p}$ and an
asteroid with mass $m$. Usually, the mass of the asteroid $m$ is much smaller
than $m_{0}$ and $m_{p}$, so that it is reasonable to approximate the asteroid
as a test particle.

To describe the orbits, we introduce an inertial coordinate system originated
at the central body, where the $x$–$y$ plane is aligned with the invariant
plane of the system, the $x$-axis points towards an arbitrary direction in the
invariant plane, and the $z$-axis goes along the angular momentum vector of
the planet moving around the central body. Under this reference frame, the
state of the asteroid (or planet) is characterized by orbital elements: the
semimajor axis $a$ ($a_{p}$), eccentricity $e$ ($e_{p}$), inclination $I$
($I_{p}$), longitude of the ascending node $\Omega$ ($\Omega_{p}$), argument
of pericenter $\omega$ ($\omega_{p}$) and the mean anomaly $M$ ($M_{p}$) or
the true anomaly $f$ ($f_{p}$). Unless otherwise stated, in the entire study
we will adopt the notations without any subscript to stand for the elements of
the asteroid and use the ones with subscript $p$ to represent the elements of
the planet.

The dynamical model taken in this investigation can be treated as a perturbed
Keplerian problem, in which the motion of the asteroid moving around the
central body is perturbed by the gravitational attraction of the planet. The
disturbing function, governing the evolution of the asteroid, can be written
as (Murray & Dermott, 1999)

${\cal R}={\cal
G}{m_{p}}\left({\frac{1}{\Delta}-\frac{r}{{r_{p}^{2}}}\cos\psi}\right),$ (1)

where ${\cal G}$ is the universal gravitational constant, $r$ and $r_{p}$ are,
respectively, the radial distances of the asteroid and planet relative to the
central body. The distance between the asteroid and planet, $\Delta$, is
determined by

$\Delta={\left({{r^{2}}+r_{p}^{2}-2r{r_{p}}\cos\psi}\right)^{1/2}}$

where $\psi$ is the separation angle between the asteroid and planet relative
to the central body, given by

$\cos{\psi}={\sin^{2}}\frac{I}{2}\cos\left({f+{\theta_{p}}+\omega-\Omega}\right)+{\cos^{2}}\frac{I}{2}\cos\left({f-{\theta_{p}}+\omega+\Omega}\right),$
(2)

with ${\theta_{p}}={f_{p}}+{\omega_{p}}+{\Omega_{p}}$ as the true longitude of
the planet.

Usually, the disturbing function can be separated into the direct part
standing for the direct gravitational attraction from the planet and the
indirect part corresponding to the perturbation of the planet to the origin of
the selected coordinate system, given by (Ellis & Murray, 2000)

${\cal R}={\cal G}{m_{p}}\left({{{\cal R}_{D}}+{{\cal R}_{I}}}\right),$ (3)

where ${\cal R}_{D}$ and ${\cal R}_{I}$ are

${{\cal R}_{D}}=\frac{1}{\Delta},\quad{{\cal
R}_{I}}=-\frac{r}{{r_{p}^{2}}}\cos\psi.$ (4)

##  3 A new expansion of planetary disturbing function

Before expanding the planetary disturbing function, it is necessary to
introduce two small parameters with the same order of magnitude as the
eccentricities of the asteroid and planet, given by (Murray & Dermott, 1999;
Ellis & Murray, 2000)

$\varepsilon=\frac{r}{a}-1\sim{\cal
O}(e),\quad{\varepsilon_{p}}=\frac{{{r_{p}}}}{{{a_{p}}}}-1\sim{\cal
O}(e_{p}).$ (5)

With these two small parameters, the direct part of disturbing function can be
expanded around $\varepsilon=\varepsilon_{p}=0$ (i.e. around circular orbits)
as a formal series in the following form:

$\displaystyle{{\cal R}_{D}}$
$\displaystyle=\sum\limits_{n=0}^{N}{\sum\limits_{m=0}^{n}{\frac{{{\varepsilon^{m}}\varepsilon_{p}^{n-m}}}{{n!}}{a^{m}}a_{p}^{n-m}\frac{{{{\rm
d}^{n}}}}{{{\rm d}{a^{m}}{\rm
d}a_{p}^{n-m}}}\left({\frac{1}{{{\Delta_{0}}}}}\right)}}$ (6)
$\displaystyle=\sum\limits_{n=0}^{N}{\sum\limits_{m=0}^{n}{\frac{{{{\left({r-a}\right)}^{m}}{{\left({{r_{p}}-{a_{p}}}\right)}^{n-m}}}}{{n!}}\frac{{{{\rm
d}^{n}}}}{{{\rm d}{a^{m}}{\rm
d}a_{p}^{n-m}}}\left({\frac{1}{{{\Delta_{0}}}}}\right)}},$

where $N$ specifies the truncated order in terms of
$\max\left(\varepsilon,{\varepsilon_{p}}\right)$ or
$\max\left(e,{e_{p}}\right)$, and the term $\frac{1}{{{\Delta_{0}}}}$ is given
by

$\frac{1}{{{\Delta_{0}}}}={\left[{{a^{2}}+a_{p}^{2}-2a{a_{p}}\cos\psi}\right]^{-\frac{1}{2}}}.$
(7)

Obviously, the term $\frac{1}{{{\Delta_{0}}}}$ stands for the mutual
interaction between two inclined circular orbits.

It is to be noted that, in the classical Laplace-type expansions, the first
step is to expand the disturbing function around a prograde or retrograde
coplanar orbit (Murray & Dermott, 1999; Ellis & Murray, 2000; Morais &
Namouni, 2013) or around an inclined reference orbit (Namouni & Morais, 2017b,
a) as a power series of inclination and then, in the second step, the
disturbing function is expanded around circular orbits as power series of
eccentricities. Thus, the step of series expansion represented by equation (6)
is in accordance with the second step of the traditional Laplace-type
expansions. Through this step of expansion, the disturbing function between
two elliptic orbits transforms into a summation of disturbing functions
between circular orbits.

Observing equation (6), we can find that the difficulty of expanding the
direct part of disturbing function ${{\cal R}_{D}}$ lies in the expansion of
the term $\frac{1}{{{\Delta_{0}}}}$. Let us remind that, in the classical
Laplace-type expansions, the term $\frac{1}{{{\Delta_{0}}}}$ given by equation
(7) is expanded by using two-dimensional Laplace coefficients as follows
(Namouni & Morais, 2017a):

$\displaystyle\frac{1}{{{\Delta_{0}}}}$
$\displaystyle=\frac{1}{{{a_{p}}}}\frac{1}{{\sqrt{1+{\alpha^{2}}-2\alpha\cos\psi}}}$
$\displaystyle=\frac{1}{{{a_{p}}}}\sum\limits_{\scriptstyle-\infty<j,k<\infty\hfill\atop\scriptstyle\bmod(j+k,2)=0\hfill}{\frac{1}{4}b_{1/2}^{jk}(\alpha,I)\cos\left[{j\left({\Omega-{\lambda_{p}}}\right)+k\left({f+\omega}\right)}\right]},$

where $\alpha=a/a_{p}$ is the semimajor axis ratio between asteroid and planet
and $b_{1/2}^{jk}(\alpha,I)$ are the two-dimensional Laplace coefficients. In
practical simulations, the Laplace coefficients $b_{1/2}^{jk}(\alpha,I)$ are
approximated by series expansions in $\alpha$ up to order $N_{\alpha}$. It is
known that the Laplace coefficients are divergent when $\alpha$ is approaching
unity ($\alpha=1$ corresponds to co-orbital configuration). Therefore, the
resulting expansions of disturbing function based on Laplace coefficients are
divergent in co-orbital or nearly co-orbital configurations.

To describe all types of configurations including the interior ($\alpha<1$),
exterior ($\alpha>1$) and co-orbital ($\alpha=1$) resonances, it is required
to develop a ‘general’ expansion for the planetary disturbing function which
can be used for configurations with arbitrary semimajor axis ratios. To this
end, we introduce a new small parameter to make it be possible to perform
Taylor expansion for the planetary disturbing function whose convergence
depends on the magnitude of the introduced small parameter solely. In the
entire process, we avoid to use Laplace coefficients, making the new expansion
of disturbing function be easy to realize in computer code.

For convenience, we organize equation (7) in the following form:

$\frac{1}{{{\Delta_{0}}}}=\frac{1}{{a+{a_{p}}}}{\left[{1-\frac{{2a{a_{p}}}}{{{{\left({a+{a_{p}}}\right)}^{2}}}}\left({1+\cos\psi}\right)}\right]^{-\frac{1}{2}}}=\frac{1}{{a+{a_{p}}}}{\left[{1-x}\right]^{-\frac{1}{2}}},$
(8)

where the variable $x$ is given by

$x=\frac{{2a{a_{p}}}}{{{{\left({a+{a_{p}}}\right)}^{2}}}}\left({1+\cos\psi}\right).$
(9)

Evidently, the domain of definition of $x$ is $x\in\left[{0,1}\right]$. In
particular, the variable $x$ is equal to zero when the separation angle is
$\psi=\pi$, and the variable $x$ is equal to unity when $\psi=0$ and $a=a_{p}$
(the condition with $\psi=0$ and $a=a_{p}$ corresponds to collision points for
two circular inclined orbits). Obviously, equation (8) is an increasing
function of $x$, so that the direct part of disturbing function, ${{\cal
R}_{D}}$, is also an increasing function of $x$, meaning that there is a
positive correlation between ${{\cal R}_{D}}$ and $x$. Thus, for a given pair
of $(a,a_{p})$, when the two objects involved are located on the opposite
sides of the central star ($\psi=\pi$), the disturbing function ${{\cal
R}_{D}}$ takes the minimum and, when the two objects are located on the same
side ($\psi=0$), the disturbing function ${{\cal R}_{D}}$ takes the maximum.

Considering the domain of definition of $x$, we can perform Taylor expansion
for equation (8) around the origin $x=0$ as follows:

$\frac{1}{{{\Delta_{0}}}}=\frac{1}{{a+{a_{p}}}}\sum\limits_{k=0}^{{k_{\max}}}{\frac{{\left({2k-1}\right)!!}}{{\left({2k}\right)!!}}{x^{k}}},$

where $k_{\max}$ specifies the truncated order of the expansion in terms of
$x$. Except for the collision point at which $x=1$, such a Taylor expansion is
uniformly convergent. However, when the variable $x$ is relatively large or
close to unity, the convergence of the Taylor expansion becomes very slow.
Thus, this is not an ideal expansion.

Observing the expression given by equation (9), we find that the variable $x$
oscillates around a fixed value for two circular orbits specified by a given
pair of semimajor axes $(a,a_{p})$. Naturally, this fixed value, denoted by
$x_{c}$, can be taken as a reference point for $x$ and then equation (8) can
be expanded in the vicinity of $x=x_{c}$ as a Taylor series in the following
form:

$\frac{1}{{{\Delta_{0}}}}=\frac{1}{{a+{a_{p}}}}\sum\limits_{k=0}^{{k_{\max}}}{\frac{{\left({2k-1}\right)!!}}{{\left({2k}\right)!!}}\frac{{{{\left({x-{x_{c}}}\right)}^{k}}}}{{{{\left({1-{x_{c}}}\right)}^{1/2+k}}}}},$
(10)

where $k_{\max}$ specifies the truncated order of the expansion in terms of
the deviation of $x$ relative to its reference point, namely $\delta
x=x-x_{c}$.

Obviously, the accuracy of the Taylor expansion is determined by the magnitude
of $\delta x$ and the truncated order $k_{\max}$. A better choice of the
reference point $x_{c}$ results in a smaller magnitude of $\delta x$, implying
that a lower truncated order $k_{\max}$ is required to achieve a given
accuracy. Thus, the current difficulty lies in the determination of the
reference point. Observing the expression of equation (9), we can find a good
choice for the reference point, given by

${x_{c}}=\frac{{2{a^{0}}{a_{p}^{0}}}}{{{{\left({a^{0}+{a_{p}^{0}}}\right)}^{2}}}}=\frac{{2{\alpha_{0}}}}{{{{\left({1+{\alpha_{0}}}\right)}^{2}}}},$
(11)

where the initial semimajor axis ratio, $\alpha_{0}$, is defined by
$\alpha_{0}=a^{0}/a_{p}^{0}$, where $a^{0}$ and $a_{p}^{0}$ are, respectively,
the initial values of the semimajor axes of the asteroid and planet. It is not
difficult to observe that the deviation of $x$ relative to its reference point
$x_{c}$ is

$\delta
x=x-{x_{c}}=\frac{{2\alpha}}{{{{\left({1+\alpha}\right)}^{2}}}}\left({1+\cos\psi}\right)-\frac{{2{\alpha_{0}}}}{{{{\left({1+{\alpha_{0}}}\right)}^{2}}}}\approx\frac{{2\alpha}}{{{{\left({1+\alpha}\right)}^{2}}}}\cos\psi$

which satisfies $\delta x\in[-0.5,0.5]$ regardless of the values of semimajor
axis ratios and inclinations, ensuring that the Taylor expansion given by
equation (10) is uniformly convergent. Remember that this is a key feature of
our expansion method, making the final expansion be valid for arbitrary
inclinations and semimajor axis ratios.

In particular, when the asteroid is located inside the $p_{0}$:$q_{0}$ mean
motion resonance with respect to the planet, the reference point for $x$ can
be approximated as

${x_{c}}=\frac{{2{{\left({\frac{{{q_{0}}}}{{{p_{0}}}}}\right)}^{2/3}}}}{{{{\left[{1+{{\left({\frac{{{q_{0}}}}{{{p_{0}}}}}\right)}^{2/3}}}\right]}^{2}}}}.$
(12)

In order to validate the Taylor expansion performed by equation (10), in
Figure 1 we make an accuracy analysis for the Taylor expansion of the function
$f(x)=\frac{1}{\sqrt{1-x}}$ (equation (8) shows
$\frac{1}{\Delta_{0}}=\frac{1}{a+a_{p}}f(x)$) truncated at orders
$k_{\max}=10,20,30$ in the case of the reference point at $x_{c}=0$ (origin)
and $x_{c}=\frac{{2{\alpha_{0}}}}{{{{\left({1+{\alpha_{0}}}\right)}^{2}}}}$.
The relative deviation is defined by

$\Delta f(x)=\frac{\left|f_{\rm approx}-f_{\rm accurate}\right|}{f_{\rm
accurate}},$

where $f_{\rm approx}$ is the magnitude obtained by Taylor expansion and
$f_{\rm accurate}$ stands for the exact value of $f(x)$. From Figure 1, we can
see that, at the same truncated order, the expansion around
$x_{c}=\frac{{2{\alpha_{0}}}}{{{{\left({1+{\alpha_{0}}}\right)}^{2}}}}$ has
better accuracy than that around $x_{c}=0$. This is in agreement with our
expectation.

Figure 1: Accuracy analysis for the Taylor expansion of the function
$f(x)=\frac{1}{\sqrt{1-x}}$ in the case of the reference point at $x_{c}=0$
(_left panel_) and
$x_{c}=\frac{{2{\alpha_{0}}}}{{{{\left({1+{\alpha_{0}}}\right)}^{2}}}}$
(_right panel_). $\Delta f(x)$ stands for the relative deviation of $f(x)$
(please see the text for its definition). The red points stand for the
location of reference points. We take the 3:1 MMR as an example
($\alpha_{0}=0.48$) to do this analysis.

Substituting equation (9) into equation (10) and performing Newton’s binomial
expansion to the result can lead to the following expression:

$\frac{1}{{{\Delta_{0}}}}=\sum\limits_{k=0}^{{k_{\max}}}{\sum\limits_{q=0}^{k}{\sum\limits_{l=0}^{q}{\frac{{{2^{q}}\left({2k-1}\right)!!}}{{\left({2k}\right)!!}}{{\left({-1}\right)}^{k-q}}{k\choose
q}{q\choose
l}\frac{{x_{c}^{k-q}}}{{{{\left({1-{x_{c}}}\right)}^{1/2+k}}}}\frac{{{a^{q}}a_{p}^{q}}}{{{{\left({a+{a_{p}}}\right)}^{2q+1}}}}{{\cos}^{l}}\psi}}},$
(13)

where $i\choose j$ is the binomial coefficient, defined by ${i\choose
j}=\frac{{i!}}{{j!(i-j)!}}$. For the sake of brevity, we denote the function
associated with the semimajor axes $a$ and $a_{p}$ by

${f_{q}}(a,{a_{p}})=\frac{{{a^{q}}a_{p}^{q}}}{{{{\left({a+{a_{p}}}\right)}^{2q+1}}}}$
(14)

and define the differential operator with respect to the semimajor axes of the
asteroid and planet as

${{\rm D}_{{k_{1}},{k_{2}}}}=\frac{{{{\rm d}^{{k_{1}}+{k_{2}}}}}}{{{\rm
d}{a^{{k_{1}}}}{\rm d}a_{p}^{{k_{2}}}}}.$ (15)

Through some mathematical derivation, the high-order partial derivatives of
${f_{q}}(a,{a_{p}})$ with respect to $a$ and $a_{p}$ can be calculated by

${{\rm
D}_{{k_{1}},{k_{2}}}}{f_{q}}=\sum\limits_{{l_{2}}=0}^{{k_{2}}}{\sum\limits_{{l_{1}}=0}^{{k_{1}}}{\frac{{{{\left({-1}\right)}^{{l_{2}}+{l_{1}}}}\left({2q+{l_{2}}+{l_{1}}}\right)!}}{{{{\left({a+{a_{p}}}\right)}^{2q+1+{l_{2}}+{l_{1}}}}\left({2q}\right)!}}{{k_{2}}\choose{l_{2}}}{{k_{1}}\choose{l_{1}}}\left[{\frac{{{{\rm
d}^{{k_{1}}-{l_{1}}}}}}{{{\rm
d}{a^{{k_{1}}-{l_{1}}}}}}{a^{q}}}\right]\left[{\frac{{{{\rm
d}^{{k_{2}}-{l_{2}}}}}}{{{\rm d}a_{p}^{{k_{2}}-{l_{2}}}}}a_{p}^{q}}\right]}}.$
(16)

Replacing equation (13) in equation (6) and combining the differential
operator defined by equation (15), we can obtain the direct part of disturbing
function as follows:

$\displaystyle{{\cal R}_{D}}$
$\displaystyle=\sum\limits_{n=0}^{N}{\sum\limits_{m=0}^{n}{\sum\limits_{k=0}^{{k_{\max}}}{\sum\limits_{q=0}^{k}{\sum\limits_{l=0}^{q}{\frac{{{2^{q}}\left({2k-1}\right)!!}}{{\left({2k}\right)!!}}{{\left({-1}\right)}^{k-q}}{k\choose
q}{q\choose
l}\frac{{x_{c}^{k-q}}}{{{{\left({1-{x_{c}}}\right)}^{1/2+k}}}}}}}}}$ (17)
$\displaystyle\times\frac{{{{\left({r-a}\right)}^{m}}{{\left({{r_{p}}-{a_{p}}}\right)}^{n-m}}}}{{n!}}\left[{{{\rm
D}_{m,n-m}}{f_{q}}\left({a,{a_{p}}}\right)}\right]{\cos^{l}}\psi,$

where the term ${\cos^{l}}\psi$ can be expanded to be

$\displaystyle{\cos^{l}}\psi$
$\displaystyle=\sum\limits_{t=0}^{l}{\sum\limits_{{t_{1}}=0}^{t}{\sum\limits_{{t_{2}}=0}^{l-t}{\frac{1}{{{2^{l}}}}{l\choose
t}{t\choose{{t_{1}}}}{{l-t}\choose{{t_{2}}}}{{\sin}^{2\left({l-t}\right)}}({I\mathord{\left/{\vphantom{I2}}\right.\kern-1.2pt}2}){\cos}^{2t}({I\mathord{\left/{\vphantom{I2}}\right.\kern-1.2pt}2})}}}$
(18)
$\displaystyle\times\cos\left[{\left({l-2{t_{1}}-2{t_{2}}}\right)\left({f+\omega}\right)+\left({l-2t+2{t_{1}}-2{t_{2}}}\right)\left({{\theta_{p}}-\Omega}\right)}\right]$

and the terms ${\left({r-a}\right)^{m}}$ and
${\left({{r_{p}}-{a_{p}}}\right)^{n-m}}$ can be expanded as follows:

$\displaystyle{\left({r-a}\right)^{m}}=\sum\limits_{{t_{3}}=0}^{m}{{{\left({-1}\right)}^{m-{t_{3}}}}{a^{m}}{m\choose
t_{3}}{{\left({\frac{r}{a}}\right)}^{{t_{3}}}}},$ (19)
$\displaystyle{\left({{r_{p}}-{a_{p}}}\right)^{n-m}}=\sum\limits_{{t_{4}}=0}^{n-m}{{{\left({-1}\right)}^{n-m-{t_{4}}}}{n-m\choose
t_{4}}a_{p}^{n-m}{{\left({\frac{{{r_{p}}}}{{{a_{p}}}}}\right)}^{{t_{4}}}}}.$

By substituting Eqs. (18) and (19) into equation (17), the direct part of
disturbing function becomes

$\displaystyle{{\cal R}_{D}}$
$\displaystyle=\sum\limits_{n=0}^{N}{\sum\limits_{m=0}^{n}{\sum\limits_{k=0}^{{k_{\max}}}{\sum\limits_{q=0}^{k}{\sum\limits_{l=0}^{q}{\sum\limits_{t=0}^{l}{\sum\limits_{{t_{1}}=0}^{t}{\sum\limits_{{t_{2}}=0}^{l-t}{\sum\limits_{{t_{3}}=0}^{m}{\sum\limits_{{t_{4}}=0}^{n-m}{\kappa\left[{{{\rm
D}_{m,n-m}}{f_{q}}}\right]{a^{m}}a_{p}^{n-m}\frac{{{{\sin}^{2l}}\left({{I\mathord{\left/{\vphantom{I2}}\right.\kern-1.2pt}2}}\right)}}{{{{\tan}^{2t}}\left({{I\mathord{\left/{\vphantom{I2}}\right.\kern-1.2pt}2}}\right)}}}}}}}}}}}}$
(20)
$\displaystyle\times{\left({\frac{r}{a}}\right)^{{t_{3}}}}{\left({\frac{{{r_{p}}}}{{{a_{p}}}}}\right)^{{t_{4}}}}\cos\left[{\left({l-2{t_{1}}-2{t_{2}}}\right)\left({f+\omega}\right)+\left({l-2t+2{t_{1}}-2{t_{2}}}\right)\left({{\theta_{p}}-\Omega}\right)}\right]$

where $\kappa$ is given by

$\kappa=\frac{{{2^{q-l}}x_{c}^{k-q}\left({2k-1}\right)!!}}{{{{\left({1-{x_{c}}}\right)}^{1/2+k}}n!\left({2k}\right)!!}}{\left({-1}\right)^{n+k-q-{t_{3}}-{t_{4}}}}{k\choose
q}{q\choose l}{t\choose{{t_{1}}}}{{l-t}\choose{{t_{2}}}}{l\choose
t}{m\choose{{t_{3}}}}{{n-m}\choose{{t_{4}}}}.$

The trigonometric functions related to the true anomaly arising in equation
(20) can be transformed into the forms of the mean anomaly by the following
relation (Kaula, 1961, 1962):

$\displaystyle{\left({\frac{r}{a}}\right)^{{t_{3}}}}{\left({\frac{{{r_{p}}}}{{{a_{p}}}}}\right)^{{t_{4}}}}\cos\left[{\left({l-2{t_{1}}-2{t_{2}}}\right)\left({f+\omega}\right)+\left({l-2t+2{t_{1}}-2{t_{2}}}\right)\left({{\theta_{p}}-\Omega}\right)}\right]$
(21)
$\displaystyle=\sum\limits_{{s_{1}}=-\infty}^{\infty}{\sum\limits_{{s_{2}}=-\infty}^{\infty}{X_{{s_{1}}}^{{t_{3}},\left({l-2{t_{1}}-2{t_{2}}}\right)}\left(e\right)X_{{s_{2}}}^{{t_{4}},\left({l-2t+2{t_{1}}-2{t_{2}}}\right)}\left({{e_{p}}}\right)}}$
$\displaystyle\times\cos\left[{{s_{1}}M+{s_{2}}{M_{p}}+\left({l-2{t_{1}}-2{t_{2}}}\right)\omega+\left({l-2t+2{t_{1}}-2{t_{2}}}\right)\left({{\varpi_{p}}-\Omega}\right)}\right]$

The mathematical expression given by equation (21) has been used in previous
studies (Murray & Dermott, 1999; Ellis & Murray, 2000; Beaugé & Michtchenko,
2003; Beaugé et al., 2006; Lei et al., 2019; Lei, 2019a, b). The Hansen
coefficient $X_{c}^{a,b}(e)$ arising in equation (21) is a function of the
eccentricity $e$, given by (Hughes, 1981)

$X_{c}^{a,b}(e)={e^{\left|{c-b}\right|}}\sum\limits_{s=0}^{\infty}{Y_{s+t,s+u}^{a,b}{e^{2s}}},$

where $t=\max(0,c-b)$, $u=\max(0,b-c)$, and $Y_{s+t,s+u}^{a,b}$ is the Newcomb
operator calculated in a recurrence manner (Hughes, 1981; Murray & Dermott,
1999; Ellis & Murray, 2000). As for the Hansen coefficient $X_{c}^{a,b}(e)$,
its lowest power in eccentricity is ${\left|{c-b}\right|}$. It should be
mentioned that the expansion given by equation (21) is convergent under the
condition that the eccentricities ($e$ and $e_{p}$) are smaller than the
critical value $e_{c}=0.6627$ (Wintner, 1947).

Substituting equation (21) into equation (20), we can organize the direct part
of the disturbing function as a Fourier series in terms of the classical
angles $M$, $M_{p}$, $\Omega$, $\omega$ and ${{\varpi_{p}}}$ in the following
form:

$\displaystyle{{\cal R}_{D}}$
$\displaystyle=\sum\limits_{n=0}^{N}{\sum\limits_{m=0}^{n}{\sum\limits_{k=0}^{{k_{\max}}}{\sum\limits_{q=0}^{k}{\sum\limits_{l=0}^{q}{\sum\limits_{t=0}^{l}{\sum\limits_{{t_{1}}=0}^{t}{\sum\limits_{{t_{2}}=0}^{l-t}{\sum\limits_{{t_{3}}=0}^{m}{\sum\limits_{{t_{4}}=0}^{n-m}{\sum\limits_{{s_{1}}=-\infty}^{\infty}{\sum\limits_{{s_{2}}=-\infty}^{\infty}{\kappa\times\left[{{{\rm
D}_{m,n-m}}{f_{q}}}\right]}}}}}}}}}}}}$ (22)
$\displaystyle\times{a^{m}}a_{p}^{n-m}X_{{s_{1}}}^{{t_{3}},\left({l-2{t_{1}}-2{t_{2}}}\right)}\left(e\right)X_{{s_{2}}}^{{t_{4}},\left({l-2t+2{t_{1}}-2{t_{2}}}\right)}\left({{e_{p}}}\right)\frac{{{{\sin}^{2l}}\left({{I\mathord{\left/{\vphantom{I2}}\right.\kern-1.2pt}2}}\right)}}{{{{\tan}^{2t}}\left({{I\mathord{\left/{\vphantom{I2}}\right.\kern-1.2pt}2}}\right)}}$
$\displaystyle\times\cos\left[{{s_{1}}M+{s_{2}}{M_{p}}+\left({l-2{t_{1}}-2{t_{2}}}\right)\omega+\left({l-2t+2{t_{1}}-2{t_{2}}}\right)\left({{\varpi_{p}}-\Omega}\right)}\right].$

For the indirect part, its Fourier series form can be easily derived to be

$\displaystyle{{\cal R}_{I}}$
$\displaystyle=-\frac{a}{{a_{p}^{2}}}\sum\limits_{{s_{1}}=-\infty}^{\infty}{\sum\limits_{{s_{2}}=-\infty}^{\infty}{X_{{s_{1}}}^{1,1}\left(e\right)X_{{s_{2}}}^{-2,1}\left({{e_{p}}}\right)}}$
(23)
$\displaystyle\times\left[\begin{array}[]{l}{\sin^{2}}\frac{I}{2}\cos\left({{s_{1}}M+{s_{2}}{M_{p}}+\omega+{\varpi_{p}}-\Omega}\right)\\\
+{\cos^{2}}\frac{I}{2}\cos\left({{s_{1}}M-{s_{2}}{M_{p}}+\omega-{\varpi_{p}}+\Omega}\right)\end{array}\right].$

At last, we arrive at the final literal expansion of the planetary disturbing
function in the following form:

$\displaystyle{\cal R}$ $\displaystyle={\cal
G}{m_{p}}\sum\limits_{n=0}^{N}{\sum\limits_{m=0}^{n}{\sum\limits_{k=0}^{{k_{\max}}}{\sum\limits_{q=0}^{k}{\sum\limits_{l=0}^{q}{\sum\limits_{t=0}^{l}{\sum\limits_{{t_{1}}=0}^{t}{\sum\limits_{{t_{2}}=0}^{l-t}{\sum\limits_{{t_{3}}=0}^{m}{\sum\limits_{{t_{4}}=0}^{n-m}{\sum\limits_{{s_{1}}=-\infty}^{\infty}{\sum\limits_{{s_{2}}=-\infty}^{\infty}\kappa}}}}}}}}}}}$
(24) $\displaystyle\times\left[{{{\rm
D}_{m,n-m}}{f_{q}}}\right]{a^{m}}a_{p}^{n-m}X_{{s_{1}}}^{{t_{3}},\left({l-2{t_{1}}-2{t_{2}}}\right)}\left(e\right)X_{{s_{2}}}^{{t_{4}},\left({l-2t+2{t_{1}}-2{t_{2}}}\right)}\left({{e_{p}}}\right)\frac{{{{\sin}^{2l}}\left({{I\mathord{\left/{\vphantom{I2}}\right.\kern-1.2pt}2}}\right)}}{{{{\tan}^{2t}}\left({{I\mathord{\left/{\vphantom{I2}}\right.\kern-1.2pt}2}}\right)}}$
$\displaystyle\times\cos\left[\begin{array}[]{l}{s_{1}}\lambda+{s_{2}}{\lambda_{p}}+\left({l-2{t_{1}}-2{t_{2}}-{s_{1}}}\right)\varpi\\\
+\left({l-2t+2{t_{1}}-2{t_{2}}-{s_{2}}}\right){\varpi_{p}}+\left({2t+4{t_{2}}-2l}\right)\Omega\end{array}\right]$
$\displaystyle-\frac{{{\cal
G}{m_{p}}}}{{{a_{p}}}}\left({\frac{a}{{{a_{p}}}}}\right)\sum\limits_{{s_{1}}=-\infty}^{\infty}{\sum\limits_{{s_{2}}=-\infty}^{\infty}{X_{{s_{1}}}^{1,1}\left(e\right)X_{{s_{2}}}^{-2,1}\left({{e_{p}}}\right)}}$
$\displaystyle\times\left\\{\begin{array}[]{l}{\sin^{2}}\frac{I}{2}\cos\left[{{s_{1}}\lambda+{s_{2}}{\lambda_{p}}+\left({1-{s_{1}}}\right)\varpi+\left({1-{s_{2}}}\right){\varpi_{p}}-2\Omega}\right]\\\
+{\cos^{2}}\frac{I}{2}\cos\left[{{s_{1}}\lambda-{s_{2}}{\lambda_{p}}+\left({1-{s_{1}}}\right)\varpi-\left({1-{s_{2}}}\right){\varpi_{p}}}\right]\end{array}\right\\}$

where $\varpi=\omega+\Omega$ and $\varpi_{p}=\omega_{p}+\Omega_{p}$ are the
longitudes of pericenter, and $\lambda=M+\varpi$ and
${\lambda_{p}}={M_{p}}+{\varpi_{p}}$ are the mean longitudes of the asteroid
and planet, respectively. It should be mentioned that the expansion given by
equation (24) is convergent when $e<0.6627$ and $e_{p}<0.6627$ are satisfied
due to the convergence of the transformation between mean and true anomalies
given by equation (21).

In the final expansion, the harmonic arguments are linear combinations of the
mean longitudes, longitudes of pericenter and longitudes of ascending node of
the asteroid and planet, denoted by

$\varphi={k_{1}}\lambda+{k_{2}}\varpi+{k_{3}}\Omega+{k_{4}}{\lambda_{p}}+{k_{5}}{\varpi_{p}}+{k_{6}}{\Omega_{p}},$

where the coefficients $k_{i}(i=1,2,...6)$ are all integers. As a result, the
expansion of disturbing function can be written in an elegant form:

${\cal R}=\sum\limits_{{\bm{k}}}{{{\cal
C}_{{\bm{k}}}}\left({a,e,I,{a_{p}},{e_{p}}}\right)\cos{\varphi}},$ (25)

where ${\bm{k}}=\left\\{{k_{1}},{k_{2}},{k_{3}},{k_{4}},{k_{5}}\right\\}$ and
the explicit expressions of the force amplitudes ${{{\cal
C}_{{\bm{k}}}}\left({a,e,I,{a_{p}},{e_{p}}}\right)}$ can be directly derived
from equation (24). Observing the literal expansion given by equation (24), we
can summarize the following properties:

  * •

The expansion of the planetary disturbing function contains only cosine terms
of argument $\varphi$.

  * •

The coefficient $k_{6}$ is zero, meaning that the perturber’s longitude of
ascending node vanishes from the disturbing function. This is due to the
choice of the reference plane as the orbit of the perturber (i.e. the
inclination of the perturber’s orbit is zero).

  * •

The d’Alembert relation holds for the coefficients in the cosine arguments,
namely $\sum\limits_{i}{{k_{i}}=0}$.

  * •

The coefficient of the longitude of ascending node $\Omega$ is always an even
number.

  * •

In the expression of the force amplitude, the power of $\sin(I/2)$ (or
$\cos(I/2)$) is always an even number, indicating that the force amplitude is
an even function with respect to the mutual inclination.

  * •

The relationship between the coefficient of $\Omega$ and the power of
$\sin(\frac{I}{2})$ holds $\left|2t+4t_{2}-2l\right|\leq 2(l-t)$ (see the
expression of the direct part), meaning that the power of $\sin(\frac{I}{2})$
arising in the force amplitude is greater than or equal to the absolute value
of the coefficient of $\Omega$ in the cosine argument (such relationship also
holds for the indirect part).

Regarding the expression of Hansen coefficient, we can determine its lowest
order terms in eccentricity as $X_{c}^{a,b}\left(e\right)={\cal
O}\left({{e^{\left|{b-c}\right|}}}\right)$. Applying this property to the
coefficients appearing in equation (24) yields

$\displaystyle
X_{{s_{1}}}^{{t_{3}},\left({l-2{t_{1}}-2{t_{2}}}\right)}\left(e\right)={\cal
O}\left({{e^{\left|{l-2{t_{1}}-2{t_{2}}-{s_{1}}}\right|}}}\right),$
$\displaystyle
X_{{s_{2}}}^{{t_{4}},\left({l-2t+2{t_{1}}-2{t_{2}}}\right)}\left({{e_{p}}}\right)={\cal
O}\left({{e_{p}^{\left|{l-2t+2{t_{1}}-2{t_{2}}-{s_{2}}}\right|}}}\right),$
$\displaystyle X_{{s_{1}}}^{1,1}\left(e\right)={\cal
O}\left({{e^{\left|{1-{s_{1}}}\right|}}}\right),\quad
X_{{s_{2}}}^{-2,1}\left({{e_{p}}}\right)={\cal
O}\left({e_{p}^{\left|{1-{s_{2}}}\right|}}\right).$

Thus, it is not difficult to observe that, in the expansion given by equation
(24), the lowest powers of $e$ and $e_{p}$ arising in the force amplitude are
the absolute values of the coefficients of $\varpi$ and $\varpi_{p}$,
respectively (i.e. $\left|k_{2}\right|$ and $\left|k_{5}\right|$). According
to this property, it is possible to see that (a) the amplitude of a certain
harmonic argument has positive correlation with the eccentricity $e$ ($e_{p}$)
if the number $k_{2}$ ($k_{5}$) is not equal to zero, and (b) when the
eccentricity $e$ ($e_{p}$) becomes zero, the terms associated with the
arguments $\varphi$ with non-zero $k_{2}$ ($k_{5}$) would vanish from the
disturbing function.

In particular, when the inclinations and eccentricities are small (in this
situation the classical Laplace-type expansion works very well), it is
possible to reserve the lowest-order terms of $e$ ($e_{p}$) and $\sin(I/2)$ in
equation (24) to approximate the disturbing function:

${\cal R}\approx\sum\limits_{\bm{k}}{{{\tilde{\cal
C}}_{\bm{k}}}\left({a,{a_{p}}}\right){e^{\left|{{k_{2}}}\right|}}e_{p}^{\left|{{k_{5}}}\right|}{{\sin}^{\left|{{k_{3}}}\right|}}(I/2)\cos\left({{k_{1}}\lambda+{k_{2}}\varpi+{k_{3}}\Omega+{k_{4}}{\lambda_{p}}+{k_{5}}{\varpi_{p}}}\right)},$
(26)

which has a similar formal expression to the ones given in previous studies
(Murray & Dermott, 1999; Morbidelli, 2002).

Especially, when the planet is assumed to move around the central body in a
circular orbit (i.e. $e_{p}=0$), the resulting dynamical model, describing the
motion of the asteroid, reduces to the circular restricted three-body problem
(CRTBP) (Szebehely, 1967). In this case, the Hansen coefficients
$X_{{s_{2}}}^{{t_{4}},\left({l-2t+2{t_{1}}-2{t_{2}}}\right)}\left({{e_{p}}}\right)$
in the expansion given by equation (24) are different from zero only when
$s_{2}$ is equal to ${l-2t+2{t_{1}}-2{t_{2}}}$, and the coefficients
${X_{{s_{2}}}^{-2,1}\left({{e_{p}}}\right)}$ are different from zero only when
$s_{2}$ is equal to unity. Consequently, the disturbing function for
describing the motion of asteroids in the CRTBP can be simplified to be

$\displaystyle{\cal R}$ $\displaystyle=\frac{{{\cal
G}{m_{p}}}}{{{a_{p}}}}\sum\limits_{n=0}^{N}{\sum\limits_{k=0}^{{k_{\max}}}{\sum\limits_{q=0}^{k}{\sum\limits_{l=0}^{q}{\sum\limits_{m=0}^{n}{\sum\limits_{t=0}^{l}{\sum\limits_{{t_{1}}=0}^{t}{\sum\limits_{{t_{2}}=0}^{l-t}{\sum\limits_{s=-\infty}^{\infty}{{\kappa_{0}}X_{s}^{m,\left({l-2{t_{1}}-2{t_{2}}}\right)}\left(e\right)\frac{{{\alpha^{n}}}}{{n!}}}}}}}}}}}$
(27) $\displaystyle\times\left[{{{\rm
D}_{n}}{f_{q}}(\alpha)}\right]\frac{{{{\sin}^{2l}}\left({\frac{I}{2}}\right)}}{{{{\tan}^{2t}}\left({\frac{I}{2}}\right)}}\cos\left[\begin{array}[]{l}s\lambda+\left({l-2t+2{t_{1}}-2{t_{2}}}\right){\lambda_{p}}\\\
+\left({l-2{t_{1}}-2{t_{2}}-s}\right)\varpi+\left({2t-2l+4{t_{2}}}\right)\Omega\end{array}\right]$
$\displaystyle-\frac{{{\cal
G}{m_{p}}}}{{{a_{p}}}}\left({\frac{a}{{{a_{p}}}}}\right)\sum\limits_{s=-\infty}^{\infty}{X_{s}^{1,1}\left(e\right)\left\\{\begin{array}[]{l}{\sin^{2}}({I\mathord{\left/{\vphantom{I2}}\right.\kern-1.2pt}2})\cos\left[{s\lambda+{\lambda_{p}}+\left({1-s}\right)\varpi-2\Omega}\right]\\\
+{\cos^{2}}({I\mathord{\left/{\vphantom{I2}}\right.\kern-1.2pt}2})\cos\left[{s\lambda-{\lambda_{p}}+\left({1-s}\right)\varpi}\right]\end{array}\right\\}}$

where $\alpha$ is the semimajor axis ratio between the asteroid and planet
(i.e. $\alpha=a/a_{p}$) and $\kappa_{0}$ is

${\kappa_{0}}=\frac{{{{\left({-1}\right)}^{n-m+k-q}}{2^{q-l}}\left({2k-1}\right)!!}}{{\left({2k}\right)!!}}{k\choose
q}{q\choose l}{n\choose m}{t\choose{{t_{1}}}}{{l-t}\choose{{t_{2}}}}{l\choose
t}\frac{{x_{c}^{k-q}}}{{{{\left({1-{x_{c}}}\right)}^{1/2+k}}}}.$

In practical simulations, the upper limit of $\left|s\right|$ arising in the
summations is also taken as 20. In equation (27), the functions associated
with $\alpha$ are denoted by

${f_{q}}(\alpha)=\frac{{{\alpha^{q}}}}{{{{\left({1+\alpha}\right)}^{2q+1}}}},\quad{{\rm
D}_{n}}{f_{q}}(\alpha)=\frac{{{{\rm d}^{n}}}}{{{\rm
d}{\alpha^{n}}}}{f_{q}}(\alpha),$

where the high-order derivatives are calculated by

${{\rm
D}_{n}}{f_{q}}(\alpha)=\sum\limits_{m=0}^{n}{{{\left({-1}\right)}^{m}}\frac{{\left({2q+m}\right)!}}{{\left({2q}\right)!}}{n\choose
m}{{\left({1+\alpha}\right)}^{-2q-1-m}}\left[{\frac{{{{\rm d}^{n-m}}}}{{{\rm
d}{\alpha^{n-m}}}}{\alpha^{q}}}\right]}.$

Remember that, in equation (24) or (27), there are two numbers characterizing
the truncated orders of series expansion: $N$ and $k_{\max}$. The first number
$N$ specifies the maximum order in terms of the small parameter:
$\max\left(\varepsilon,{\varepsilon_{p}}\right)$. As discussed before,
$\varepsilon$ and $\varepsilon_{p}$ have the same order of magnitude as the
eccentricities of the asteroid and planet, thus the number $N$ corresponds to
the order in terms of $\max\left(e,e_{p}\right)$. The second number $k_{\max}$
determines the maximum order of the Taylor expansion in the deviation of $x$
relative to its reference point, namely $\delta x=x-x_{c}$.

Note that the convergence of the new expansion developed in this study is not
restricted by the mutual inclination and the semimajor axis ratio between the
asteroid and planet because the small parameter arising in the series
expansion always satisfies $\delta x\in[-0.5,0.5]$ regardless of mutual
inclination and semimajor axis ratio. Thus, it becomes possible to use such a
‘general’ expansion of disturbing function to formulate dynamical models in
describing the interior ($\alpha<1$), co-orbital ($\alpha\approx 1$) and
exterior ($\alpha>1$) mean motion resonances of minor bodies with arbitrary
inclinations in planetary systems.

##  4 Resonant disturbing function

In this section, we assume that the planet moves in a circular orbit around
the central body (the model corresponds to the well-known CRTBP). In this
case, the disturbing function given by equation (27) can be utilized to
formulate resonant models. Regarding an asteroid located inside the
$p_{0}$:$q_{0}$ mean motion resonance with respect to the planet, there is
more than one critical argument (with different $k_{0}$), given by

$\sigma_{{k_{0}}}^{{p_{0}}:{q_{0}}}={q_{0}}\lambda-{p_{0}}{\lambda_{p}}+\left({{p_{0}}-{q_{0}}}\right)\Omega+{k_{0}}\omega,$
(28)

where $k_{0}$ has the same parity of $p_{0}+q_{0}$. In previous works, several
special resonances associated with $p_{0}$:$q_{0}$ resonances have been
discussed, such as the pure eccentricity (or inclination) resonances and
retrograde resonances (Murray & Dermott, 1999; Morais & Namouni, 2013; Namouni
& Morais, 2017a). When $k_{0}={{p_{0}}-{q_{0}}}$, the angle defined by
equation (28) becomes the critical argument of pure eccentricity resonance,
denoted by
$\sigma_{1}={q_{0}}\lambda-{p_{0}}{\lambda_{p}}+\left({{p_{0}}-{q_{0}}}\right)\varpi$.
When $k_{0}=0$ (it requires that $p_{0}+q_{0}$ is even), the angle defined by
equation (28) represents the critical argument of pure inclination resonance,
denoted by
$\sigma_{2}={q_{0}}\lambda-{p_{0}}{\lambda_{p}}+\left({{p_{0}}-{q_{0}}}\right)\Omega$
(it is noted that, when $p_{0}-q_{0}=1$, the pure inclination argument is
defined by
$\sigma_{2}=2({q_{0}}\lambda-{p_{0}}{\lambda_{p}}+\left({{p_{0}}-{q_{0}}}\right)\Omega)$).
When $k_{0}=-p_{0}-q_{0}$, the angle defined by equation (28) corresponds to
the critical argument of the retrograde (pure eccentricity) resonance, denoted
by
$\sigma_{3}={q_{0}}\lambda-{p_{0}}{\lambda_{p}}+\left({{p_{0}}-{q_{0}}}\right)\Omega-\left({{p_{0}}+{q_{0}}}\right)\omega$.
It is known that, for the minor bodies located on coplanar or nearly coplanar
orbits, the pure eccentricity resonance with $\sigma_{1}$ as the critical
argument has the dominant strength and, for the bodies moving in retrograde
coplanar orbits, the retrograde resonance with $\sigma_{3}$ as the critical
argument has the dominant strength. However, when the minor bodies moves on an
inclined orbit, the force amplitudes of the arguments
$\sigma_{{k_{0}}}^{{p_{0}}:{q_{0}}}$ with different $k_{0}$ arising in the
disturbing function may have comparable amplitudes, indicating that it is not
suitable to use a single critical argument in describing the resonant
behavior. Please refer to Figure 16 in Lei (2019b) for detailed discussions on
single critical argument and characteristic angle for a certain mean motion
resonance.

Among all the critical arguments of the $p_{0}$:$q_{0}$ resonance, there is a
common part, given by
$\sigma={q_{0}}\lambda-{p_{0}}{\lambda_{p}}+\left({{p_{0}}-{q_{0}}}\right)\Omega$,
which corresponds to the special case with $k_{0}=0$ in equation (28). In Lei
(2019b), the common angle $\sigma$ is defined as the characteristic angle of
the $p_{0}$:$q_{0}$ resonance (notice that the expression of $\sigma$ is
determined by the numbers $p_{0}$ and $q_{0}$).

According to the mathematical expression of $\sigma$, the mean longitude of
the minor body can be expressed by means of $\sigma$, $\lambda_{p}$ and
$\Omega$ as

$\lambda=\frac{1}{{{q_{0}}}}\left[{\sigma+{p_{0}}{\lambda_{p}}-\left({{p_{0}}-{q_{0}}}\right)\Omega}\right].$

As a result, the disturbing function given by equation (27) becomes a function
of the angles $\sigma$, $\omega$, $\Omega$ and $\lambda_{p}$. For a minor body
located inside a mean motion resonance, the associated resonant angle
($\sigma$) becomes a semi-secular angle variable compared to $\lambda_{p}$. In
addition, it is known that the angles $\omega$ and $\Omega$ are slow variables
in the timescale of mean motion resonance (Gallardo, 2019, 2020). Thus, among
all these angles ($\sigma$, $\omega$, $\Omega$, $\lambda_{p}$), only the mean
longitude of the planet ($\lambda_{p}$) is a fast variable. Usually, the terms
in the disturbing function involving the fast angle $\lambda_{p}$ produce
short-term influences upon the motion of the minor body. During the long-term
evolution of the minor body, it is often to remove those short-term effects by
means of averaging technique. Such an averaging process is called “secular
approximation” (Naoz et al., 2013). Alternatively, the short-period effects
can be averaged out by means of von Zeiple’s method (Brouwer, 1959) or
Hori–Deprit method (Hori, 1966; Deprit, 1969). As a result, filtering out the
short-term effects (i.e. averaging the disturbing function over $\lambda_{p}$)
can be achieved by

${{\cal
R}^{*}}\left(a,e,I,{\sigma,\omega}\right)=\frac{1}{{2q_{0}\pi}}\int\limits_{0}^{2q_{0}\pi}{{\cal
R}\left(a,e,I,{\sigma,\omega,\Omega,{\lambda_{p}}}\right){\rm
d}{\lambda_{p}}},$ (29)

which leads to the resonant disturbing function. After averaging the
disturbing function over $q_{0}$ orbital periods of the planet, both the angle
variables $\lambda_{p}$ and $\Omega$ disappear from the averaged expression,
so that the number of degree of freedom is reduced by two. As a result, the
resonant model determined by ${{\cal
R}^{*}}\left(a,e,I,{\sigma,\omega}\right)$ is of two degrees of freedom.
However, there is only one integral of motion (Hamiltonian of system), showing
that the resulting two degree-of-freedom dynamical model determined by ${{\cal
R}^{*}}\left(a,e,I,{\sigma,\omega}\right)$ is not integrable.

Replacing equation (27) in equation (29), we can easily obtain the resonant
disturbing functions associated with the interior ($p_{0}>q_{0}$ and
$\alpha<1$), co-orbital ($p_{0}=q_{0}$ and $\alpha\approx 1$) and exterior
($p_{0}<q_{0}$ and $\alpha>1$) mean motion resonances for minor bodies with
arbitrary inclinations. According to the number of $p_{0}$, there are two
cases for the resonant disturbing function: $p_{0}=1$ and $p_{0}\neq 1$. When
$p_{0}=1$, the resonant disturbing function for describing the 1:$q_{0}$
resonances becomes

$\displaystyle{{\cal R}^{{*}}}$ $\displaystyle=\frac{{{\cal
G}{m_{p}}}}{{{a_{p}}}}\sum\limits_{n\geq 0}^{N}{\sum\limits_{k\geq
0}^{{k_{\max}}}{\sum\limits_{q\geq 0}^{k}{\sum\limits_{l\geq
0}^{q}{\sum\limits_{m\geq 0}^{n}{\sum\limits_{t\geq
0}^{l}{\sum\limits_{{t_{1}}\geq 0}^{t}{\sum\limits_{{t_{2}}\geq
0}^{l-t}{{\kappa_{0}}\times
X_{{q_{0}}\left({2t-l-2{t_{1}}+2{t_{2}}}\right)}^{m,\left({l-2{t_{1}}-2{t_{2}}}\right)}\left(e\right)}}}}}}}}\frac{{{\alpha^{n}}}}{{n!}}$
(30) $\displaystyle\times\left[{{{\rm
D}_{n}}{f_{q}}(\alpha)}\right]\frac{{{{\sin}^{2l}}\left({\frac{I}{2}}\right)}}{{{{\tan}^{2t}}\left({\frac{I}{2}}\right)}}\cos\left\\{\begin{array}[]{l}\left({2t-l-2{t_{1}}+2{t_{2}}}\right)\sigma\\\
+\left[{\left({l-2{t_{1}}-2{t_{2}}}\right)-{q_{0}}\left({2t-l-2{t_{1}}+2{t_{2}}}\right)}\right]\omega\end{array}\right\\}$
$\displaystyle-\frac{{{\cal
G}{m_{p}}}}{{{a_{p}}}}\left({\frac{a}{{{a_{p}}}}}\right)\left\\{\begin{array}[]{l}X_{{q_{0}}}^{1,1}\left(e\right){\cos^{2}}({I\mathord{\left/{\vphantom{I2}}\right.\kern-1.2pt}2})\cos\left[{\sigma+\left({1-{q_{0}}}\right)\omega}\right]\\\
+X_{-{q_{0}}}^{1,1}\left(e\right){\sin^{2}}({I\mathord{\left/{\vphantom{I2}}\right.\kern-1.2pt}2})\cos\left[{\sigma-\left({1+{q_{0}}}\right)\omega}\right]\end{array}\right\\}$

where $\sigma={q_{0}}\lambda-{\lambda_{p}}+\left({1-{q_{0}}}\right)\Omega$
and, when $p_{0}\neq 1$ (in this case, the indirect part has no contribution
to the resonant disturbing function), the resonant disturbing function can be
organized as

$\displaystyle{{\cal R}^{*}}$ $\displaystyle=\frac{{{\cal
G}{m_{p}}}}{{{a_{p}}}}\sum\limits_{n\geq 0}^{N}{\sum\limits_{k\geq
0}^{{k_{\max}}}{\sum\limits_{q\geq 0}^{k}{\sum\limits_{l\geq
0}^{q}{\sum\limits_{m\geq 0}^{n}{\sum\limits_{t\geq
0}^{l}{\sum\limits_{{t_{1}}\geq 0}^{t}{\sum\limits_{\scriptstyle{t_{2}}\geq
0\atop\scriptstyle\bmod\left[{\left({2t-l-2{t_{1}}+2{t_{2}}}\right),{p_{0}}}\right]=0}^{l-t}{{\kappa_{0}}}}}}}}}}$
(31) $\displaystyle\times\frac{{{\alpha^{n}}}}{{n!}}\left[{{{\rm
D}_{n}}{f_{q}}(\alpha)}\right]X_{\frac{{{q_{0}}}}{{{p_{0}}}}\left({2t-l-2{t_{1}}+2{t_{2}}}\right)}^{m,\left({l-2{t_{1}}-2{t_{2}}}\right)}\left(e\right)\frac{{{{\sin}^{2l}}\left({\frac{I}{2}}\right)}}{{{{\tan}^{2t}}\left({\frac{I}{2}}\right)}}$
$\displaystyle\times\cos\left\\{{\frac{1}{{{p_{0}}}}\left({2t-l-2{t_{1}}+2{t_{2}}}\right)\sigma+\left[{\left({l-2{t_{1}}-2{t_{2}}}\right)-\frac{{{q_{0}}}}{{{p_{0}}}}\left({2t-l-2{t_{1}}+2{t_{2}}}\right)}\right]\omega}\right\\}$

where $\bmod(A,B)=0$ means that $A$ is divisible by $B$. It is noted that the
expansion of resonant disturbing function given by equation (30) or (31) can
be easily realized in computer language. By the way, direct numerical
integration of equation (29) can produce numerical results of the resonant
disturbing function, as adopted by Gallardo (2006, 2019, 2020) in their
studies. In the direct numerical integration, the original form of the
disturbing function given by equation (1) can be used. In the current work, we
will take the numerical integration results as references in order to validate
our analytical expansion.

Obviously, the resonant disturbing function ${\cal R}^{*}$ given by equation
(30) or (31) is a function of $a$, $e$, $I$, $\sigma$ and $\omega$. For the
sake of brevity, the resonant disturbing function can be denoted by an elegant
form,

${{\cal
R}^{*}\left(a,e,I,{\sigma,\omega}\right)}=\sum\limits_{k=0}^{\infty}{\sum\limits_{{k_{1}}}{{\cal
C}_{k,{k_{1}}}^{\cal R}(a,e,I)\cos\left({k\sigma+{k_{1}}\omega}\right)}},$
(32)

where ${k_{1}}\in\mathbb{Z}$ and it has the same parity of $k(p_{0}+q_{0})$.
The coefficients ${{\cal C}_{{k},{k_{1}}}^{\cal R}}(a,e,I)$ are derived from
the resonant disturbing function given by equation (30) or (31).

Figure 2: Relative errors of the new expansion truncated at different orders
in eccentricity for approximating the resonant disturbing function. Relative
errors are measured by taking the numerical results produced by direct
numerical integration as the accurate value of resonant disturbing function.
The new expansion is truncated at orders $N=2,3,4$ in eccentricity and, in
each case, the expansion order in $\delta x$ is fixed at $k_{\max}=30$. The
left panel is for Jupiter’s inner 3:1 resonance and the right one is for
Jupiter’s inner 2:1 resonance. In simulations, the eccentricity is fixed at
$e=0.3$, the inclination at $I=60^{\circ}$, the argument of pericenter at
$\omega=90^{\circ}$ and the semimajor axis is taken as the value at the
resonant center (simulations with other parameters can be performed in a
similar manner).

To validate the new expansion shown by Eqs. (30) and (31), it needs to discuss
the precision of the analytical developments truncated at different orders. To
this end, the new expansion is truncated at orders $N=2,3,4$ in eccentricity
and $k_{\max}=30$ in $\delta x$ and they are used to approximate the resonant
disturbing function associated with Jupiter’s interior 3:1 and 2:1 resonances.
In simulations, the truncated order of $\delta x$ is assumed at $k_{\max}=30$,
the eccentricity is fixed at $e=0.3$, the inclination is taken as
$I=60^{\circ}$ and the argument of pericenter is taken as $\omega=90^{\circ}$
(these parameters are chosen artificially). The precision of $f(={\cal
R}^{*})$ is measured by the relative error, defined by

${\Delta f}=\frac{{\left|{{f_{\rm approximate}}-{f_{\rm
accurate}}}\right|}}{{\left|{{f_{\rm accurate}}}\right|}}$

where $f_{\rm approximate}$ and $f_{\rm accurate}$ represents the approximate
and accurate values of $f$, respectively. $f_{\rm approximate}$ is computed by
analytical expansion and $f_{\rm accurate}$ is produced by direct numerical
integration. The relative errors of ${\cal R}^{*}$ are reported in Figure 2,
which shows that, in general, the series expansion truncated at a higher order
$N$ in eccentricity has a lower relative error (or higher precision) than the
lower-order expansions. This is expected by us.

Figure 3: Comparisons of the numerical results obtained by direct numerical
integration and the analytical results obtained by the new expansion for
computing the resonant disturbing function associated with the Jupiter’s inner
3:1 and 2:1 resonances (_upper panels_) and exterior 1:2 and 1:3 resonances
(_bottom panels_). The numerical results are taken as references for measuring
the accuracy of our analytical developments. The new expansion is truncated at
the fourth order in eccentricity (i.e. $N=4$) and at order $k_{\max}=30$ in
$\delta x$. In simulations, the eccentricity is assumed at $e=0.3$, the
inclination is fixed at $I=60^{\circ}$, the argument of pericenter is set as
$\omega=90^{\circ}$ and the semimajor axis is taken as the value of the
resonant center.

Then, we compare the analytical results obtained by the new expansion
truncated at order $N=4$ in eccentricity and $k_{\max}=30$ in $\delta x$ with
the numerical results obtained through direct numerical integration.
Concerning the resonant disturbing function associated with Jupiter’s inner
3:1 and 2:1 resonances as well as Jupiter’s exterior 1:2 and 1:3 resonances,
in Figure 3 we make a direct comparison between the numerical results with the
analytical results obtained by the new expansion. Please refer to the caption
of Figure 3 for the detailed setting of other parameters. From Figure 3, we
can observe: (a) our new expansion could catch the peak and valley positions
of resonant disturbing function accurately in comparison to the numerical
results (notice that the valley and peak positions correspond to the stable
and unstable equilibrium points of the resonant model and the difference of
resonant disturbing function evaluated at the peak and valley positions,
$\Delta{\cal R}^{*}$, showing the resonant strength, is positively correlated
to the resonant width, as discussed in Gallardo (2020)); (b) from a
quantitative viewpoint, the analytical results obtained by the new expansion
are in coincident with the numerical results in the entire interval of
$\sigma$.

Figure 4: Resonant disturbing function for co-orbital motion in the
Sun–Jupiter system under different settings of $\omega$ ($\omega$ is taken as
$0^{\circ}$, $45^{\circ}$, $90^{\circ}$ and $135^{\circ}$ from the top left
panel to the bottom right panel), computed by means of analytical and
numerical approaches. The numerical results are shown in blue dotted lines and
the analytical results obtained by the new expansion truncated at order $N=4$
in eccentricity and $k_{\max}=30$ in $\delta x$ are shown in red lines. In
simulations, the eccentricity is assumed at $e=0.3$, the inclination at
$I=60^{\circ}$ and the semimajor axis is taken as the value at the resonant
center. The regions in shadow represent the minimum distance between the
orbits of the asteroid and planet is smaller than $3R_{H}$ where $R_{H}$ is
the Hill radius of the planet. It is observed that outside the shadow region
the analytical results agree well with the numerical results.

At last, we apply the new expansion truncated at orders $N=4$ and
$k_{\max}=30$ to the co-orbital resonances with Jupiter ($\alpha=1$), as shown
in Figure 4. The numerical results are also provided for the convenience of
comparison. In Figure 4, the domain of $\sigma$ is marked in shadow if the
minimum distance between the asteroid and planet is smaller than three times
of the Jupiter’s Hill radius $R_{H}$. From Figure 4, it is observed that (a)
in the shadow area the analytical developments underestimate the resonant
disturbing function, this is because the analytical results are always smaller
than the numerical results in the shadow region, and (b) outside the shadow
area the analytical results produced by the new expansion are in good
agreement with the numerical results. It is known that, when the asteroid is
located in the shadow area, the perturbation theory usually fails to work due
to the strong perturbation coming from the planet. Thus, when we are using
perturbation theory to study the long-term dynamical behaviors we need to
ensure that the asteroids considered are located outside the shadow area.

In summary, the comparisons made in Figures 3 and 4 strongly support the
conclusion that the new expansion is valid for minor bodies located inside the
interior, co-orbital and exterior resonances (for better accuracy it is
required that the minimum distance between the asteroid and planet is greater
than $3R_{H}$). In the coming section, we will use our analytical developments
to formulate Hamiltonian model for mean motion resonances.

##  5 Hamiltonian model of mean motion resonances

For convenience of studying the Hamiltonian dynamics, we adopt a set of
modified Delaunay variables, given by (Morbidelli, 2002)

$\displaystyle\Lambda=\sqrt{\mu
a},\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\lambda=M+\varpi,$ (33)
$\displaystyle P=\sqrt{\mu a}\left({1-\sqrt{1-{e^{2}}}}\right),\quad\quad\quad
p=-\varpi,$ $\displaystyle Q=\sqrt{\mu a\left({1-{e^{2}}}\right)}\left({1-\cos
i}\right),\quad q=-\Omega,$
$\displaystyle{\Lambda_{p}},\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad{\lambda_{p}}=M_{p}+{\varpi}_{p},$

where $\mu=\sqrt{{\cal G}m_{0}}$ is the gravitational parameter of the central
mass, $\Lambda_{p}$ is the conjugate momentum of the mean longitude of planet
$\lambda_{p}$. With such a set of canonical variables, the Hamiltonian of
planetary system can be written as (Morbidelli, 2002)

${\cal H}=-\frac{{{\mu^{2}}}}{{2{\Lambda^{2}}}}+{n_{p}}{\Lambda_{p}}-{\cal
R}\left({\Lambda,P,Q,\lambda,p,q,{\lambda_{p}}}\right),$ (34)

where $n_{p}$ is the mean motion of the planet and the disturbing function
${\cal R}\left({\Lambda,P,Q,\lambda,p,q,{\lambda_{p}}}\right)$ is equal to
${\cal R}\left(a,e,i,\sigma,\omega\right)$, given by equation (27). In order
to describe the dynamics associated with the $p_{0}$:$q_{0}$ resonance, we
introduce a linear canonical transformation, given by

$\displaystyle{\Sigma_{{1}}}$
$\displaystyle=\frac{1}{q_{0}}\Lambda,\quad{\sigma_{{1}}}=q_{0}\lambda-{p_{0}}{\lambda_{p}}-\left({{p_{0}}-{q_{0}}}\right)q=\sigma,$
(35) $\displaystyle{\Sigma_{2}}$
$\displaystyle=-P,\quad{\sigma_{2}}=q-p=\omega,$ $\displaystyle{\Sigma_{3}}$
$\displaystyle=-P-Q-\frac{{{p_{0}}-{q_{0}}}}{q_{0}}\Lambda,\quad{\sigma_{3}}=-q=\Omega,$
$\displaystyle{\Sigma_{4}}$
$\displaystyle={\Lambda_{p}}+\frac{{{p_{0}}}}{q_{0}}\Lambda,\quad{\sigma_{4}}={\lambda_{p}},$

which can be realized by the following generating function:

${\cal
S}=\left[{{q_{0}}\lambda-\left({{p_{0}}-{q_{0}}}\right)q-{p_{0}}{\lambda_{p}}}\right]{\Sigma_{1}}+\left({q-p}\right){\Sigma_{2}}-q{\Sigma_{3}}+{\lambda_{p}}{\Sigma_{4}}.$

As a result, the Hamiltonian can be written as follows:

${\cal
H}=-\frac{{{\mu^{2}}}}{{2{{\left({q_{0}{\Sigma_{1}}}\right)}^{2}}}}+{n_{p}}\left({{\Sigma_{4}}-{p_{0}}{\Sigma_{1}}}\right)-{\cal
R}\left({{\Sigma_{1}},{\Sigma_{2}},{\Sigma_{3}},{\Sigma_{4}},{\sigma_{1}},{\sigma_{2}},{\sigma_{3}},{\sigma_{4}}}\right).$
(36)

In the new set of canonical variables, $\sigma_{1}=\sigma$ corresponds to the
characteristic angle of the $p_{0}$:$q_{0}$ resonance, and $\sigma_{2}=\omega$
corresponds to the argument of pericenter. It should be mentioned that the
canonical transformation given by equation (35) is valid for all the interior
($p_{0}>q_{0}$), co-orbital ($p_{0}=q_{0}$) and exterior ($p_{0}<q_{0}$)
resonances.

In particular, when the minor body is located inside a mean motion resonance,
the disturbing function is replaced by the resonant disturbing function given
by equation (32). Consequently, the averaged resonant Hamiltonian becomes

${{\cal
H}^{*}}=-\frac{{{\mu^{2}}}}{{2{{\left({q_{0}{\Sigma_{1}}}\right)}^{2}}}}-{p_{0}}{n_{p}}{\Sigma_{1}}-{{\cal
R}^{*}}\left({{\Sigma_{1}},{\Sigma_{2}},{\Sigma_{3}},{\Sigma_{4}},{\sigma_{1}},{\sigma_{2}}}\right),$
(37)

where the constant term $n_{p}\Sigma_{4}$ is eliminated from the Hamiltonian
and the canonical equations of motion can be written as

$\displaystyle{\dot{\sigma_{1}}}=\frac{{\partial{{\cal
H}^{*}}}}{{\partial{\Sigma_{1}}}},\quad{{\dot{\Sigma}}_{1}}=-\frac{{\partial{{\cal
H}^{*}}}}{{\partial{\sigma_{1}}}},$ (38)
$\displaystyle{\dot{\sigma_{2}}}=\frac{{\partial{{\cal
H}^{*}}}}{{\partial{\Sigma_{2}}}},\quad{{\dot{\Sigma}}_{2}}=-\frac{{\partial{{\cal
H}^{*}}}}{{\partial{\sigma_{2}}}},$

which correspond to a dynamical system with two degrees of freedom
($\sigma_{1}$ and $\sigma_{2}$ are the associated angular coordinates). In the
resonant model specified by equation (37), $\sigma_{3}$ and $\sigma_{4}$ are
cycle coordinates, so that their conjugate momenta $\Sigma_{3}$ and
$\Sigma_{4}$ become integrals of motion, given by

${\Sigma_{3}}=-P-Q-\frac{{{p_{0}}-{q_{0}}}}{q_{0}}\Lambda=\sqrt{\mu
a}\left[{\sqrt{1-{e^{2}}}\cos I-\frac{{{p_{0}}}}{q_{0}}}\right]={\rm const}$
(39)

and

${\Sigma_{4}}={\Lambda_{p}}+\frac{{{p_{0}}}}{q_{0}}\Lambda={\Lambda_{p}}+\frac{{{p_{0}}}}{q_{0}}\sqrt{\mu
a}={\rm const}.$ (40)

The integral of motion given by equation (39) shows that there is coupled
oscillations among semimajor axis ratio, eccentricity and inclination of the
test particle, meaning that exchange between the Keplerian energy
($K=-\frac{\mu}{2a}$) and the angular momentum along the $z$ direction
($H=\sqrt{\mu a(1-e^{2})}\cos{I}$) occurs in the long-term evolution.

In terms of the classical elements, the averaged resonant Hamiltonian given by
equation (37) can be further expressed as

$\displaystyle{\cal H^{*}}=$
$\displaystyle-\frac{\mu}{{2a}}-{n_{p}}\frac{{{p_{0}}}}{q_{0}}\sqrt{\mu
a}-{{\cal R}^{*}}\left({a,e,i,\sigma,\omega}\right)$ (41) $\displaystyle=$
$\displaystyle-\frac{\mu}{{2a}}-{n_{p}}\frac{{{p_{0}}}}{q_{0}}\sqrt{\mu
a}-\sum\limits_{k=0}^{\infty}{\sum\limits_{{k_{1}}}{{\cal C}_{k,{k_{1}}}^{\cal
R}(a,e,I)\cos\left({k\sigma+{k_{1}}\omega}\right)}}$

where the resonant disturbing function ${{\cal
R}^{*}}\left({a,e,i,\sigma,\omega}\right)$ is given by equation (32). It
should be noted that the form given by equation (41) is analogous to the one
given by Gallardo (2020) (see equation 5 in his work). The difference is that
the resonant disturbing function in the present work is produced by analytical
developments while, in Gallardo (2020), the resonant disturbing function is
obtained by means of direct numerical integration.

If we are only interested in the dynamics during the timescale of mean motion
resonances, the argument of pericenter $\omega$ can be assumed as a constant
(Gallardo, 2006, 2019, 2020; Lei, 2019b), so that the dynamical model reduces
to a system with a single degree of freedom, specified by the resonant
Hamiltonian

${{\cal H}^{*}}=-\frac{\mu}{{2a}}-{n_{p}}\frac{{{p_{0}}}}{{{q_{0}}}}\sqrt{\mu
a}-\sum\limits_{k=0}^{\infty}{\left[{{{\cal
C}_{k}}\left({a,e,I,\omega}\right)\cos k\sigma+{{\cal
S}_{k}}\left({a,e,I,\omega}\right)\sin k\sigma}\right]},$ (42)

where the coefficients are given by

$\displaystyle{{\cal
C}_{k}}\left({a,e,I,\omega}\right)=\sum\limits_{{k_{1}}}{{\cal
C}_{k,{k_{1}}}^{\cal R}\left({a,e,I}\right)\cos\left({{k_{1}}\omega}\right)},$
$\displaystyle{{\cal
S}_{k}}\left({a,e,I,\omega}\right)=-\sum\limits_{{k_{1}}}{{\cal
C}_{k,{k_{1}}}^{\cal R}\left({a,e,I}\right)\sin\left({{k_{1}}\omega}\right)}.$

In the single degree-of-freedom resonant model, the equations of motion are

${\dot{\sigma}_{1}}=\frac{{\partial{{\cal
H}^{*}}}}{{\partial{\Sigma_{1}}}},\quad{\dot{\Sigma}_{1}}=-\frac{{\partial{{\cal
H}^{*}}}}{{\partial{\sigma_{1}}}},$ (43)

which can be written in an equivalent form:

$\displaystyle\frac{{{\rm d}\sigma}}{{{\rm d}t}}=$ $\displaystyle
2{q_{0}}\sqrt{{a\mathord{\left/{\vphantom{a\mu}}\right.\kern-1.2pt}\mu}}\frac{{\partial{{\cal
H}^{*}}}}{{\partial a}}=\frac{{{q_{0}}\sqrt{\mu
a}}}{{{a^{2}}}}-{n_{p}}{p_{0}}-2{q_{0}}\sqrt{{a\mathord{\left/{\vphantom{a\mu}}\right.\kern-1.2pt}\mu}}\frac{{\partial{{\cal
R}^{*}}}}{{\partial a}},$ (44) $\displaystyle\frac{{{\rm d}a}}{{{\rm d}t}}=$
$\displaystyle-2{q_{0}}\sqrt{{a\mathord{\left/{\vphantom{a\mu}}\right.\kern-1.2pt}\mu}}\frac{{\partial{{\cal
H}^{*}}}}{{\partial\sigma}}=2{q_{0}}\sqrt{{a\mathord{\left/{\vphantom{a\mu}}\right.\kern-1.2pt}\mu}}\frac{{\partial{{\cal
R}^{*}}}}{{\partial\sigma}}.$

The second equation of equation (44) indicates that the dependence of the
resonant disturbing function ${\cal R}^{*}$ (or Hamiltonian ${\cal H}^{*}$) on
$\sigma$ determines the dynamical behavior of semimajor axis, as pointed out
in Gallardo (2020). Under the resonant model, the equilibrium points can be
obtained by solving the frozen equation:

$\frac{{{\rm d}\sigma}}{{{\rm d}t}}=\frac{{{\rm d}a}}{{{\rm d}t}}=0.$ (45)

The first equation $\frac{{{\rm d}\sigma}}{{{\rm d}t}}\approx
q_{0}{\dot{\lambda}}-p_{0}{\dot{\lambda}_{p}}=0$ shows that the equilibrium
points are approximately located at

$a={a_{0}}={\left({\frac{{q_{0}^{2}}}{{p_{0}^{2}}}\frac{\mu}{{n_{p}^{2}}}}\right)^{\frac{1}{3}}}.$

The stability condition of equilibrium points requires that the second
derivative of ${\cal H}^{*}$ with respect to $\sigma$ evaluated at the
considered point is smaller than zero, expressed by

$\frac{{{\partial^{2}}{{\cal
H}^{*}}}}{{\partial{\sigma^{2}}}}=\sum\limits_{k=1}^{\infty}{{k^{2}}\left[{{{\cal
C}_{k}}\left({a,e,I,\omega}\right)\cos k\sigma+{{\cal
S}_{k}}\left({a,e,I,\omega}\right)\sin k\sigma}\right]}\;{{<}}\;{{0}},$ (46)

meaning that the stable equilibrium points correspond to the local maxima of
the resonant Hamiltonian (or local minima of the resonant disturbing
function). Usually, the stable equilibrium points, denoted by
$(a_{0},\sigma_{s})$, correspond to the resonant centers and the unstable
ones, denoted by $(a_{0},\sigma_{u})$, correspond to saddle points. The level
curves of the resonant Hamiltonian passing through the saddles points provide
the dynamical separatrices, which divide the phase space into libration and
circulation regions. According to the discussion in Gallardo (2020), the
libration region is specified by the so-called resonance’s half width $\Delta
a=a_{\rm sep}-a_{0}$ ($a_{\rm sep}$ is defined by the separatrix), calculated
by

$\displaystyle\Delta a$ $\displaystyle=a_{\rm sep}-a_{0}$ (47)
$\displaystyle=\frac{{2\sqrt{6}}}{{3n}}{\left[{{{\cal
H}^{*}}\left({{a_{0}},{e_{0}},{I_{0}},{\sigma_{s}},{\omega_{0}}}\right)-{{\cal
H}^{*}}\left({{a_{0}},{e_{0}},{I_{0}},{\sigma_{u}},{\omega_{0}}}\right)}\right]^{\frac{1}{2}}}$
$\displaystyle=\frac{{2\sqrt{6}}}{{3n}}{\left[{{{\cal
R}^{*}}\left({{a_{0}},{e_{0}},{I_{0}},{\sigma_{u}},{\omega_{0}}}\right)-{{\cal
R}^{*}}\left({{a_{0}},{e_{0}},{I_{0}},{\sigma_{s}},{\omega_{0}}}\right)}\right]^{\frac{1}{2}}}$
$\displaystyle=\frac{{2\sqrt{6}}}{{3n}}{\left\\{{\sum\limits_{k=0}^{\infty}{\left[{{{\cal
C}_{k}}\left({\cos k{\sigma_{u}}-\cos k{\sigma_{s}}}\right)+{{\cal
S}_{k}}\left({\sin k{\sigma_{u}}-\sin
k{\sigma_{s}}}\right)}\right]}}\right\\}^{\frac{1}{2}}}$

which is in agreement with the expression derived from Hamiltonian approach by
Gallardo (2020) (see equation 18 in his work) and in agreement with the
expression of resonance’s half width derived from the multi-harmonics pendulum
model by Lei (2019b) (please refer to equation 26 in his work). Remind that
the resonant disturbing function in Gallardo (2020) is produced by means of
direct numerical integration and, in the present work, the resonant disturbing
function (or resonant Hamiltonian) can be obtained by means of both the
analytical expansions and numerical integration. In Section 6, we will use
equation (47) to produce the analytical and numerical resonant widths in terms
of $\Delta a$ and compare them in order to validate our analytical
developments.

##  6 Applications

In this section, we will apply the analytical developments performed in
Sections 4 and 5 to Jupiter’s inner and co-orbital resonances and to Neptune’s
exterior mean motion resonances. In particular, we will make direct
comparisons between analytical and numerical results.

###  6.1 Jupiter’s inner and co-orbital resonances

The resonant Hamiltonian given by equation (42) shows that the resonant model
is specified by a pair of pseudo conjugate variables $(a,\sigma)$ (notice that
there is a one to one correspondence between the semimajor axis and the
conjugate variable of $\sigma$). Thus, it is possible to understand the global
dynamics of mean motion resonances by plotting the level curves of the
resonant Hamiltonian in the space spanned by $a$ and $\sigma$, which
correspond to the so-called (pseudo) phase-space structures.

Figure 5: Level curves of the resonant Hamiltonian (phase-space structures)
shown in the space spanned by the characteristic angle $\sigma$ and semimajor
axis $a$ for Jupiter’s 3:1 and 2:1 resonances. In both plots, the eccentricity
is assumed at $e=0.3$, the inclination at $I=60^{\circ}$ and the argument of
pericenter at $\omega=0^{\circ}$. The green lines stand for the orbits
numerically propagated under the Sun–Jupiter–asteroid dynamical model. The
lines shown in red correspond to the dynamical separatrices passing through
saddle points and the blue dots stand for the locations of resonant centers.
The resonant half width, measuring the distance between the separatrix and the
resonant center evaluated at $\sigma=\sigma_{s}$, is denoted by $\Delta
a=a_{\rm sep}-a_{0}$. The mathematical expression of $\Delta a$ is given by
equation (47).

In Figure 5, we report the phase-space structures associated with Jupiter’s
inner 3:1 and 2:1 resonances. For Jupiter’s inner 3:1 resonance, the
characteristic argument is defined by $\sigma=\lambda-3\lambda_{p}+2\Omega$
and, for Jupiter’s inner 2:1 resonance, the characteristic argument is defined
by $\sigma=\lambda-2\lambda_{p}+\Omega$. In simulations, the eccentricity is
fixed at $e=0.3$, the inclination at $60^{\circ}$ and the argument of
pericenter at $\omega=0^{\circ}$ (notice that the phase structures with other
settings of elements could be produced in a similar manner). In the phase-
space structures, the resonant centers (corresponding to stable equilibrium
points in the resonant model) are marked in blue dots, and the level curves
passing through the saddle points stand for the dynamical separatrices shown
in red lines, which divide the total space into two types of regions:
circulation and libration. For the 3:1 resonance, the resonant center is
located at $\sigma=\lambda-3\lambda_{p}+2\Omega=\pi$ and, for the 2:1
resonance, the resonant center is located at
$\sigma=\lambda-2\lambda_{p}+\Omega=0$. In addition, the resonant half width
measures the distance between $a_{\rm sep}$ and $a_{0}$, which is explicitly
marked by $\Delta a$. The mathematical expression of $\Delta a$ is provided by
equation (47).

To validate the phase portraits obtained from our new expansion of planetary
disturbing function, orbits with the same initial conditions are numerically
propagated under the Sun–Jupiter–asteroid system and they are shown in green
lines, as shown in Figure 5. The periodic oscillations are caused by the
short-period effects arising in the $N$-body model. Evidently, the numerical
trajectories could follow closely along the level curves of resonant
Hamiltonian, showing that our expansion of planetary disturbing function is
valid to predict resonant dynamics of asteroids.

Figure 6: Resonant widths for Jupiter’s inner 3:1 resonance as functions of
the mutual inclination $I$. In the left panel, analytical results of resonant
width are reported for the cases of
$\omega=0^{\circ},30^{\circ},60^{\circ},90^{\circ}$. The resonant width in
terms of the variation of semimajor axis is marked. In the right panel, a
comparison is made between the analytical and numerical results of resonant
width for the case of $\omega=90^{\circ}$. In both plots, the eccentricity is
fixed at $e=0.3$.

As stated in Section 5, the resonant disturbing function arising in equation
(47) can be produced by means of direct numerical integration (the associated
resonant width and center are called numerical results in the following
discussions) and by means of the new expansion represented by equation (30) or
(30) (the associated resonant width and center are called analytical results).

Figure 7: Resonant widths for Jupiter’s inner 2:1 resonance as functions of
the mutual inclination $I$. In the left panel, analytical results of resonant
width are reported for the cases of
$\omega=0^{\circ},30^{\circ},60^{\circ},90^{\circ}$. The resonant width in
terms of the variation of semimajor axis is marked. In the right panel, a
comparison is made between the analytical and numerical results of resonant
width for the case of $\omega=90^{\circ}$. In both plots, the eccentricity is
fixed at $e=0.3$.

Regarding Jupiter’s inner 3:1 and 2:1 resonances with $e=0.3$, Figures 6 and 7
present the analytical resonant widths as functions of the mutual inclination
$I$ for the cases of $\omega=0^{\circ},30^{\circ},60^{\circ},90^{\circ}$ in
the left panels and provide a direct comparison between the analytical and
numerical results of the resonant width for the case of $\omega=90^{\circ}$ in
the right panels. The resonant full width in terms of the variation of
semimajor axis (it is equal to $2\Delta a$) is explicitly marked.

From Figures 6 and 7, it is observed that (a) the influence of $\omega$ upon
the resonant width is significant, (b) when the mutual inclination is zero
(prograde) or $\pi$ (retrograde), the resonant widths with different $\omega$
are coincident (this is because, at inclination of $I=0$ or $I=\pi$, only the
prograde or retrograde pure eccentricity resonance dominates the dynamics, as
discussed in Section 4), and (c) the comparisons made in the right panels show
that the analytical results agree well with the numerical ones. It is worth
mentioning that Gallardo (2020) has studied the influence of $\omega$ upon
mean motion resonances by introducing the index of fragility: the index is
higher, the associated resonance is more fragile (i.e. the resonant structure
is more easy to be broken by changing $\omega$).

When $\omega=90^{\circ}$, we can observe from the right panels of Figures 6
and 7 that (a) the dynamical separatrices associated with the 3:1 resonance
have a gourd shape with three necks and it is found that the resonant center
switches between $0$ and $\pi$ passing through each neck (it is noted that
around the first neck there is a small interval of inclination
$I\in[37^{\circ},42^{\circ}]$ at which there are two pairs of dynamical
separatrices corresponding to two resonant centers located at $\sigma=0$ and
$\sigma=\pi$), (b) the dynamical separatrices associated with the 2:1
resonance have a complicated shape and it is found that, in the interval of
$I\in[54^{\circ},137^{\circ}]$, there are two pair of dynamical separatrices
bounding two different resonant centers located at $\sigma=0$ and $\sigma=\pi$
and, in the remaining interval, there is only one resonance center at either
$\sigma=0$ or $\sigma=\pi$, and (c) for the 2:1 resonance, the resonant widths
at $I=0$ and $I=\pi$ are slightly underestimated by the analytical
developments.

Regarding Jupiter’s inner 3:1 and 2:1 resonances with $\omega=0^{\circ}$ and
$\omega=90^{\circ}$, Lei (2019b) has developed a multi-harmonics pendulum
model for describing the resonant dynamics and reported the curves of
dynamical separatrices (see the bottom-left panels of Figures 5, 11, 18 and 20
in his work), which are in perfect agreement with the analytical (and
numerical) results shown in Figures 6 and 7 (it is noted that, in Lei (2019b),
the Laplace-type expansions of disturbing function are used for formulating
the resonant model).

Figure 8: Analytical and numerical results of the location of libration center
(_left panel_) and the associated resonant width (_right panel_) for Jupiter’s
1:1 resonance. In simulations, the eccentricity is fixed at $e=0.3$ and the
argument of pericenter is at $\omega=0^{\circ}$. The points inside shadow
regions stand for the parameters at which the minimum distance between the
asteroid and planet is smaller than three times of the planet’s Hill radius
(i.e. $R_{\min}<3R_{H}$). The resonant width in terms of the variation of
semimajor axis is marked. It is observed from the right panel that inside the
shaded area the resonant width is significantly underestimated by the
analytical developments due to the underestimation of resonant disturbing
function at the saddle points, as shown in Figure 4.

For Jupiter’s 1:1 resonance, the characteristic argument is defined by
$\sigma=\lambda-\lambda_{p}$ (this argument is also called the synodic angle
between the asteroid and planet in previous studies). In practical
applications of our analytical developments to this case, we fix the argument
of the pericenter at $\omega=0^{\circ}$ and the eccentricity at $e=0.3$
(simulations with other parameters can be performed in a similar manner). It
is known that, for the 1:1 resonance, there are asymmetric libration centers,
one of them is located around $\sigma=60^{\circ}$ (usually called the $L_{4}$
point) and the other one is located around $\sigma=-60^{\circ}$ (usually
called the $L_{5}$ point). These asymmetric libration centers are symmetric
with respect to $\sigma=\pi$, and they have the same dynamical behavior. Thus,
in the practical discussions, we only focus on the one around the usual
$L_{4}$ point. In Figure 8, the location of the resonant center as a function
of mutual inclination is shown in the left panel and the associated resonant
width is presented in the right panel. In both plots, the analytical results
are shown in red and the numerical results are shown in black. In particular,
the shadow region covers those parameters at which the minimum distance
between the asteroid and planet is smaller than three times of Jupiter’s Hill
radius (i.e. $R_{\min}<3R_{H}$). As discussed in Section 4, the resonant
disturbing function is underestimated by the analytical developments when the
asteroid is located inside the shadow region, resulting in the fact that the
resonant width is also underestimated by the analytical developments, as shown
in the right panel of Figure 8. However, when the asteroid in located outside
the shadow region, the analytical results are in good agreement with the
numerical ones. From the left panel of Figure 8, the asymmetric libration
center exists if and only if the inclination is smaller than $155^{\circ}$
and, in this region, it is observed from the right panel of Figure 8 that
there are two pairs of dynamical separatrices bounding the asymmetric
libration center. When the inclination is greater than $155^{\circ}$, the
asymmetric libration center disappears and is replaced by the symmetric
libration center at $\sigma=0$ and, in this region, there is one pair of
dynamical separatrices bounding the libration center, as shown in the right
panel of Figure 8.

###  6.2 Neptune’s exterior resonances

The analytical developments discussed in Section 5 are applied to the
Sun–Neptune system for studying the resonant dynamics of exterior mean motion
resonances with Neptune, including the 1:2, 1:3, 2:3, 3:5, 4:7, 3:7, 2:5 and
3:8 resonances.

Figure 9: Analytical and numerical results of the location of libration center
(_left panels_) and the associated resonant width (_right panels_) for
Neptune’s exterior 1:2 resonance. In the upper panels, the eccentricity is
fixed at $e=0.1$ and, in the bottom panels, the eccentricity is assumed at
$e=0.2$. The resonant width in terms of the variation of semimajor axis is
marked.

Regarding Neptune’s exterior 1:2 and 1:3 resonances, their characteristic
arguments are $\sigma=2\lambda-\lambda_{p}-\Omega$ and
$\sigma=3\lambda-\lambda_{p}-2\Omega$, respectively. Figure 9 presents the
locations of libration center and the associated resonant width as functions
of the mutual inclination for the 1:2 resonance, and Figure 10 shows the
corresponding results for the 1:3 resonance.

For the case of the 1:2 resonance as shown in Figure 9, two cases with $e=0.1$
and $e=0.2$ are taken into account. It is observed from the left panels that
the asymmetric libration centers exist when the inclination is smaller than
$140^{\circ}$ and, when the inclination is greater than $140^{\circ}$, the
asymmetric libration centers disappear and are replaced by symmetric libration
center located at $\sigma=\pi$. It is noted that in the interval
$I\in[0^{\circ},140^{\circ})$ there are two asymmetric libration centers which
are symmetric with respect to $\pi$. In our study, we only consider one of
them (smaller than $\pi$) because the other one has the same dynamical
behavior. From Figure 9, we can observe: (a) for both cases with $e=0.1$ and
$e=0.2$, the analytical results could be in good agreement with the numerical
results for the location of libration center and the associated resonant
width, (b) with inclination changing from $0^{\circ}$ to $140^{\circ}$, the
$\sigma$ of the libration center first decreases and then increases up to
$\sigma=\pi$ which is the location of the symmetric libration center, (c) in
the interval of $I\in[0^{\circ},140^{\circ})$, there are two pairs of
dynamical separatrices (corresponding to the inner and outer boundaries)
bounding the asymmetric libration center, (d) when the inclination is greater
than $140^{\circ}$, there is only one pair of separatrices bounding the
symmetric libration center, and (e) the resonant widths associated with
$e=0.2$ are greater than the ones associated with $e=0.1$, because there is a
positive correlation between the force amplitude in the expansions and the
eccentricity as discussed in Section 3.

Figure 10: Analytical and numerical results of the location of libration
center (_left panels_) and the associated resonant width (_right panels_) for
Neptune’s exterior 1:3 resonance. In the upper panels, the eccentricity is
fixed at $e=0.2$ and, in the bottom panels, the eccentricity is assumed at
$e=0.3$. The resonant width in terms of the variation of semimajor axis is
marked.

For the case of the 1:3 resonance as shown in Figure 10, two cases with
$e=0.2$ and $e=0.3$ are taken into consideration. It is observed from the left
panels of Figure 10 that, for the case of $e=0.2$, the asymmetric libration
centers disappear when the inclination is greater than $104^{\circ}$ and, for
the case of $e=0.3$, the asymmetric libration centers disappear when the
inclination is greater than $130^{\circ}$. Similar to the case of the 1:2
resonance, only one of the asymmetric libration centers is considered in this
study (the other one has the same dynamical behavior). From Figure 10, we can
observe that (a) the analytical results are in good agreement with the
numerical ones, (b) resonant width in the case of $e=0.3$ is greater than that
in the case of $e=0.2$, and (c) there are two pairs of dynamical separatrices
bounding the asymmetric libration center and one pair of separatrices bounding
the symmetric libration center.

Figure 11: Analytical results of the resonant width associated with Neptune’s
exterior 2:3, 3:5, 4:7, 3:7, 2:5 and 3:8 resonances. The eccentricity is fixed
at $e=0.2$ and the argument of pericenter is fixed at $\omega=0^{\circ}$.

Then, our analytical developments are applied to Neptune’s exterior 2:3, 3:5,
4:7, 3:7, 2:5 and 3:8 resonances with eccentricity at $e=0.2$ and argument of
pericenter at $\omega=0^{\circ}$. The resonant widths as functions of the
mutual inclination are reported in Figure 11. For these considered resonances,
all of them have symmetric libration center located at $\sigma=\pi$. From
Figure 11, it is observed that (a) for the 2:3, 3:5 and 4:7 resonances, the
resonant width is a decreasing function of the mutual inclination, (b) for the
3:7, 2:5 and 3:8 resonances, the resonant width first increases and then
decreases with the inclination, and (c) when the inclination is close to
$\pi$, the resonant widths of all the considered resonances are approaching
zero, indicating that their retrograde resonant strengths are very weak (this
point is different from that of inner resonances, as shown in the previous
subsection).

##  7 Conclusions

In this work, a new expansion of the planetary disturbing function has been
developed. In the first step, the original disturbing function is expanded
around circular orbits and it is expressed as a formal series in the
eccentricities of the objects involved. After the first step, the disturbing
function becomes a summation of the terms associated with the mutual
interaction between circular orbits. In the second step, we define a parameter
$x$ and expand the core function (i.e. the mutual interaction between circular
orbits) around the reference point $x_{c}$ as Taylor series of $\delta
x=x-x_{c}$. The disturbing function is finally organized as a Fourier series
form, where the force amplitudes are related to the semimajor axis,
eccentricity and inclination and the harmonic arguments are linear
combinations of the mean longitude, longitude of pericenter and longitude of
the ascending node of each mass. In the final expansion, there are two numbers
in specifying the orders of expansion: the truncated order in eccentricities
denoted by $N$ and the expansion order in $\delta x$ denoted by $k_{\max}$.

The advantages of the new expansion developed in the present work lies in the
following two aspects: (a) in the process of expansion, Laplace coefficients
are not used, so that the convergence problem arising from the series
expansion of Laplace coefficients in semimajor axis ratio can be avoided, and
(b) the new expansion is convergent regardless of the values of the mutual
inclination and semimajor axis ratio, so it becomes possible to utilize the
new expansion of planetary disturbing function to study the dynamics of minor
bodies located inside the interior, co-orbital and exterior resonances at
arbitrary inclinations.

Based on the new expansion of planetary disturbing function, the resonant
Hamiltonian is formulated through the linear and canonical transformations of
the modified Delaunay variables. It shows that the resonant model corresponds
to a dynamical model with two degrees of freedom ($\sigma$ and $\omega$ are
angular variables of the system). Considering the fact that the angle $\omega$
is much slower than $\sigma$, it is possible to make an assumption that
$\omega$ is a constant during the timescale of mean motion resonances, so that
the dynamical model naturally reduces to a system with a single degree of
freedom (only $\sigma$ is the angle coordinate). Under this assumption, the
level curves of resonant Hamiltonian in the space ($\sigma$, $a$), correspond
to the (pseudo) phase-space structures, which show the global dynamics of mean
motion resonances. In the Hamiltonian model of mean motion resonances, the
expression of resonant half width in terms of variation of semimajor axis is
provided.

The analytical developments are applied to Jupiter’s inner and co-orbital
resonances and Neptune’s exterior resonances. As for Jupiter’s inner 3:1 and
2:1 resonances, the analytical results of resonant width at different $\omega$
are produced and, in particular, a direct comparison is made between the
analytical and numerical results for the case of $\omega=90^{\circ}$, showing
that the analytical and numerical results are in good agreement. Regarding
Jupiter’s co-orbital resonance and Neptune’s exterior 1:2 and 2:3 resonances,
the location of asymmetric and symmetric libration centers and the associated
resonant width are produced by means of analytical and numerical approaches,
and comparisons between the analytical and numerical results shows that our
new expansion is valid. Furthermore, the analytical developments are applied
to Neptune’s exterior 2:3, 3:5, 4:7, 3:7, 2:5 and 3:8 resonances and their
resonant widths as functions of mutual inclination are reported.

###### Acknowledgements.

This work is performed with the financial support of the National Natural
Science Foundation of China (No. 12073011) and the National Key R&D Program of
China (No. 2019YFA0706601).

## References

  * Beaugé (1996) Beaugé, C. 1996, Celest. Mech. Dyn. Astron., 64, 313
  * Beaugé & Michtchenko (2003) Beaugé, C., & Michtchenko, T. 2003, Mon. Not. R. Astron. Soc., 341, 760
  * Beaugé et al. (2006) Beaugé, C., Nesvornỳ, D., & Dones, L. 2006, Astron. J., 131, 2299
  * Brouwer (1959) Brouwer, D. 1959, Astron. J., 64, 378
  * Brouwer & Clemence (1961) Brouwer, D., & Clemence, G. M. 1961, Methods of celestial mechanics (Academic Press)
  * Brown & Shook (1933) Brown, E. W., & Shook, C. A. 1933, Planetary theory (Cambridge Univ. Press, Cambridge)
  * Deprit (1969) Deprit, A. 1969, Celest. Mech., 1, 12
  * Ellis & Murray (2000) Ellis, K. M., & Murray, C. D. 2000, Icarus, 147, 129
  * Ferraz-Mello (1987) Ferraz-Mello, S. 1987, Astron. Astrophys., 183, 397
  * Ferraz-Mello & Sato (1989) Ferraz-Mello, S., & Sato, M. 1989, Astron. Astrophys., 225, 541
  * Gallardo (2006) Gallardo, T. 2006, Icarus, 184, 29
  * Gallardo (2019) Gallardo, T. 2019, Icarus, 317, 121
  * Gallardo (2020) Gallardo, T. 2020, Celest. Mech. Dyn. Astron., 132
  * Hori (1966) Hori, G.-i. 1966, Publ. Astron. Soc. Jpn., 18, 287
  * Hughes (1981) Hughes, S. 1981, Celest. Mech., 25, 101
  * Kaula (1961) Kaula, W. M. 1961, Geophysical J. Int., 5, 104
  * Kaula (1962) Kaula, W. M. 1962, Astron. J., 67, 300
  * Klioner (2000) Klioner, S. A. 2000, Celest. Mech. Dyn. Astron., 77, 215
  * Kozai (1962) Kozai, Y. 1962, Astron. J., 67, 591
  * Lei (2019a) Lei, H. 2019a, Mon. Not. R. Astron. Soc., 490, 4756
  * Lei (2019b) Lei, H. 2019b, Mon. Not. R. Astron. Soc., 487, 2097
  * Lei et al. (2018) Lei, H., Circi, C., & Ortore, E. 2018, Mon. Not. R. Astron. Soc., 481, 4602
  * Lei et al. (2019) Lei, H., Circi, C., & Ortore, E. 2019, Mon. Not. R. Astron. Soc., 485, 2731
  * Morais (1999) Morais, M. 1999, Astron. Astrophys., 350, 318
  * Morais & Namouni (2013) Morais, M., & Namouni, F. 2013, Celest. Mech. Dyn. Astron., 117, 405
  * Morbidelli (2002) Morbidelli, A. 2002, Modern celestial mechanics: aspects of solar system dynamics (Taylor & Francis)
  * Murray & Dermott (1999) Murray, C. D., & Dermott, S. F. 1999, Solar system dynamics (Cambridge university press)
  * Namouni & Morais (2017a) Namouni, F., & Morais, M. H. M. 2017a, Mon. Not. R. Astron. Soc., 474, 157
  * Namouni & Morais (2017b) Namouni, F., & Morais, M. H. M. 2017b, Mon. Not. R. Astron. Soc., 471, 2097
  * Naoz (2016) Naoz, S. 2016, Ann. Rev. Astron. Astrophys., 54, 441
  * Naoz et al. (2011) Naoz, S., Farr, W. M., Lithwick, Y., Rasio, F. A., & Teyssandier, J. 2011, Nature, 473, 187
  * Naoz et al. (2013) Naoz, S., Farr, W. M., Lithwick, Y., Rasio, F. A., & Teyssandier, J. 2013, Mon. Not. R. Astron. Soc., 431, 2155
  * Newcomb (1895) Newcomb, S. 1895, Astron. Papers Am. Ephem., 5, 5
  * Páez & Locatelli (2015) Páez, R. I., & Locatelli, U. 2015, Mon. Not. R. Astron. Soc., 453, 2177
  * Peirce (1849) Peirce, B. 1849, AJ, 1, 1
  * Petrovskaya (1970) Petrovskaya, M. 1970, Celest. Mech., 3, 121
  * Petrovskaya (1972) Petrovskaya, M. 1972, Celest. Mech., 6, 328
  * Roig et al. (1998) Roig, F., Simula, A., Ferraz-Mello, S., & Tsuchida, M. 1998, Astron. Astrophys., 339
  * Szebehely (1967) Szebehely, V. 1967, Theory of Orbits (Academic Press, New York)
  * Wintner (1947) Wintner, A. 1947, The analytical foundations of celestial mechanics (Princeton University Press)
  * Yokoyama (1994) Yokoyama, T. 1994, Celest. Mech. Dyn. Astron., 60, 387