# Dynamical structures of retrograde resonances: analytical and numerical
studies

Hanlun Lei1,2, Jian Li1,2
1 School of Astronomy and Space Science, Nanjing University, Nanjing 210023,
China
2 Key Laboratory of Modern Astronomy and Astrophysics in Ministry of
Education, Nanjing University, Nanjing 210023, China E-mail:
<EMAIL_ADDRESS><EMAIL_ADDRESS>
(Accepted. Received; in original form)

###### Abstract

In this work, retrograde mean motion resonances (MMRs) are investigated by
means of analytical and numerical approaches. Initially, we define a new
resonant angle to describe the retrograde MMRs and then perform a series of
canonical transformations to formulate the resonant model, in which the phase
portrait, resonant centre and resonant width can be analytically determined.
To validate the analytical developments, the non-perturbative analysis is made
by taking advantage of Poincaré surfaces of section. Some modifications are
introduced in the production of Poincaré sections and, in particular, it
becomes possible to make direct comparisons between the analytical and
numerical results. It is found that there exists an excellent correspondence
between the phase portraits and the associated Poincaré sections, and the
analytical results agree well with the numerical results in terms of the
resonant width and the location of resonant centre. Finally, the numerical
approach is utilized to determine the resonant widths and resonant centres
over the full range of eccentricity. In particular, seven known examples of
retrograde asteroids including 2015 BZ509, 2008 SO218, 1999 LE31, 2000 DG8,
2014 AT28, 2016 LS and 2016 JK24 are found inside the libration zones of
retrograde MMRs with Jupiter. The results obtained in this work may be helpful
for understanding the dynamical evolution for asteroids inside retrograde
MMRs.

###### keywords:

celestial mechanics–minor planets, asteroids, general–planets and satellites:
dynamical evolution and stability

††pagerange: Dynamical structures of retrograde resonances: analytical and
numerical studies–A††pubyear: 2021

##  1 Introduction

In recent years, a growing number of minor objects have been discovered on
retrograde orbits. In particular, among the 1048123 asteroids discovered so
far, there are 112 retrograde asteroids in our Solar
system111https://minorplanetcenter.net//iau/MPCORB.html, retrieved 3 February
2020. Gallardo (2019a) confirmed that, all along the Solar system, there is a
stability stripe around inclination of $\sim$$150^{\circ}$, where the
planetary perturbations produce smaller dynamical effects. A remark is that
the known examples of retrograde asteroids are concentrated around such an
inclination $i$$\sim$$150^{\circ}$ (Gallardo, 2019a). From the viewpoint of
stability, the configuration of mean motion resonance (MMR) can provide a
protection mechanism against planetary perturbations. Naturally, MMRs play a
fundamental role in the long-term stability especially for those high-
eccentricity objects in planetary systems. Thus, it becomes of great
significance to study dynamical structures of retrograde MMRs which could help
to understand the dynamical origin and evolution for those asteroids inside
retrograde MMRs.

In our Solar system, an increasing number of asteroids are found inside
retrograde MMRs. For example, Morais & Namouni (2013a) identified a set of
asteroids among Centaurs and Damocloids inside retrograde MMRs with Jupiter
and Saturn, including 2006 BZ8 (in 2/$-5$ resonance with Jupiter), 2008 SO218
(in 1/$-2$ resonance with Jupiter) and 2009 QY6 (in 2/$-3$ resonance with
Saturn). As stated by Morais & Namouni (2013a), these retrograde asteroids are
the first examples of Solar system objects in retrograde resonances. Several
years later, thanks to an improved orbit determination, Wiegert et al. (2017)
confirmed that asteroid 2015 BZ509 is the first asteroid inside a retrograde
co-orbital resonance with Jupiter and the authors further predicted that
retrograde co-orbital asteroids of Jupiter and other giant planets may be more
common than previously expected. Simulations made by Namouni & Morais (2018)
showed that asteroid 2015 BZ509 sits near the peak of co-orbital capture
efficiency. Regarding the transneptunian objects (TNOs), Gladman et al. (2009)
discovered the first retrograde TNO (2008 KV42) with an orbital inclination of
$103^{\circ}$, and then Chen et al. (2016) identified a new retrograde TNO
2011 KT19 with an inclination of $110^{\circ}$. One year later, Morais &
Namouni (2017) identified that 2011 KT19 is currently inside a 7:9 polar
resonance with Neptune. In a recent survey of retrograde asteroids performed
by Li et al. (2019), the authors confirmed that 2011 KT19 is currently in a
retrograde 7:9 resonance with Neptune.

Recently, Connors & Wiegert (2018) reported that asteroid 2007 VW266 (a
retrograde object near Jupiter’s orbit) is inside a retrograde 13:14 resonance
with Jupiter. Li et al. (2018) numerically explored the resonant behaviors of
the clones for the minor bodies among Centaurs and Damocloids and they
identified another four Centaurs (including 2006 RJ2, 2006 BZ8, 2017 SV13 and
2012 YE8), which are potential candidates inside retrograde co-orbital
resonances with Saturn. Furthermore, Li et al. (2019) numerically investigated
the possible retrograde resonant configurations of minor objects in the Solar
system and they identified 38 asteroids to be trapped inside retrograde MMRs
with planets.

To understand the dynamical differences of prograde and retrograde MMRs,
Morais & Giuppone (2012) numerically investigated the stability of prograde
and retrograde planets and they showed that retrograde planets are stable up
to distances closer to the perturber than the prograde counterpart. They
pointed out that the enhanced stability of retrograde planets with respect to
prograde planets is caused by the difference of dynamical structures between
retrograde and prograde resonances: for a certain $p$:$q$ resonance, the
strength (or width) of resonance is proportional to $e^{|p-q|}$ for the
prograde resonance while the strength is proportional to $e^{p+q}$ for the
retrograde resonance. As a continuation, Morais & Namouni (2013b) further
studied the dynamics of retrograde resonances in detail by using literal
expansion of disturbing function for non-coorbital resonances and using
numerical averaging of disturbing function for coorbital resonances, and the
numerical technique based on computing surfaces of section is taken to explore
phase-space structures in the vicinity of the retrograde 2:1, 1:1 and 1:2
resonances. Recently, in the framework of planar circular restricted three-
body problem, Li et al. (2020) studied the dynamics of the exterior retrograde
1:n resonances, and they showed that there is no asymmetric libration for the
retrograde 1:n resonances because of the dominant contribution of the first-
order harmonic arising in the literal expansion of disturbing function.

The dynamical system theory based on families of periodic orbits in the
restricted three-body problem provides a new approach to understand the
dynamics of retrograde MMRs. For example, in order to understand the capture
mechanism of asteroid 2015 BZ509 and similar co-orbital objects, Morais &
Namouni (2019) investigated families of periodic orbits of the retrograde co-
orbital problem and studied their stability and bifurcations. Recently, in the
Sun–Jupiter system, Kotoulas & Voyatzis (2020a) produced families of planar
resonant retrograde periodic orbits in the planar circular and elliptic
restricted three-body problems for the retrograde 2/1, 3/2, 4/3, 3/1, 5/3,
7/5, 4/1, 5/2 and 7/4 resonances. Focusing on the dynamics of resonant TNOs,
Kotoulas & Voyatzis (2020b) identified planar and three-dimensional retrograde
periodic orbits corresponding to the retrograde 1/2, 2/3 and 3/4 resonances
with Neptune.

Regarding retrograde MMRs, it is of significance to know (a) the place where
the resonance occurs (usually it is different from the nominal resonance
location), (b) the width that measures the size of resonance zones, and (c)
the dynamical structures in phase space. To this end, in the current work, we
investigate retrograde MMRs by means of analytical and numerical approaches.

In the analytical study, a new resonant angle is defined, and then a resonant
Hamiltonian model is formulated for retrograde MMRs by performing a series of
canonical transformations. Based on the resonant model, it is possible to
produce phase portraits, where the resonant centres and resonant widths can be
identified. In the numerical study, the non-perturbative approach based on
computing Poincaré sections is adopted. The technique of Poincaré sections has
been widely used in studying MMRs, e.g. Malhotra (1996); Winter & Murray
(1997a, b); Morais & Giuppone (2012); Morais & Namouni (2013b); Wang &
Malhotra (2017); Malhotra et al. (2018) and Malhotra & Zhang (2020). About the
choice of Poincaré sections, Wang & Malhotra (2017) made an important
improvement: recording the states of test particles at every successive
perihelion passage. As stated by Malhotra & Zhang (2020), such an improvement
yields a more direct visualization and physical interpretation about the
resonance zones arising in Poincaré sections. In the present work, we further
introduce two slight modifications in the production of Poincaré sections: (a)
recording the states of test particles when the ‘faster’ angular variable is
equal to zero and (b) taking the motion integral $\Gamma_{2}$ ($\Gamma_{2}$ is
the motion integral of the resonant model) as the conserved quantity. These
two minor changes make it possible to establish a correspondence between the
phase portraits in the resonant model and the Poincaré sections in the non-
averaged model. Most importantly, it is possible for us to make direct
comparisons between analytical and numerical results for a certain resonance.
This comparison is important in terms of validating the analytical model
formulated in this work and also understanding the structures arising in
Poincaré sections. As an extension, the numerical approach based on Poincaré
sections is utilized to identify resonant widths over the full range of
eccentricity, so the results provide global dynamical structures for the
retrograde MMRs of interest. In particular, seven known examples of asteroids
inside retrograde MMRs with Jupiter (2015 BZ509, 2008 SO218, 1999 LE31, 2000
DG8, 2014 AT28, 2016 LS and 2016 JK24) are considered and it is found that all
of them are located inside the libration zones determined by analyzing
Poincaré sections.

The remaining part of this work is organized as follows. In Section 2, the
Hamiltonian function for the planar circular restricted three-body problem is
briefly introduced, and in Section 3, the Hamiltonian model is formulated for
retrograde MMRs. In Section 4, the resonant model is applied to the retrograde
2:1 and 1:2 resonances with Jupiter and analytical results are presented. In
Section 5, the non-perturbative technique based on computing Poincaré sections
is introduced and applied to produce numerical results for retrograde 2:1 and
1:2 resonances with Jupiter. Direct comparisons between analytical and
numerical results are given in Section 6. In Section 7, the numerical approach
is applied to identifying resonance widths over the full range of eccentricity
for the retrograde 2:1, 1:1, 1:2, 1:3 and 1:4 resonances with Jupiter. At
last, the summary and discussion are provided in Section 8.

Figure 1: Schematic diagram for celestial longitude and latitude of the test
particle at the pericentre.

##  2 Hamiltonian function

In this work, we focus on the retrograde motion of asteroids in the planar
circular restricted three-body problem with the Sun and a planet as the
massive and secondary primaries. In this dynamical model, the primaries move
around their barycentre on circular orbits and the asteroid is regarded as a
test particle which has no gravitational influence upon the motion of the
primaries. The central body (i.e., the Sun) holds the mass of $m_{0}$ and the
planet with mass of $m_{p}$ plays the role of a perturber. All the objects are
moving on the same plane and, in particular, the perturber moves around the
Sun in counterclockwise direction and the test particle moves around the Sun
in clockwise direction (i.e., the asteroid is moving on a retrograde orbit
with inclination of $180^{\circ}$).

For the sake of computational accuracy, it is usual to normalize time and
space variables by taking the total mass of the Sun and planet as the unit of
mass, their distance as the unit of length, the orbital period of the planet
divided by $2\pi$ as the unit of time. Under such a system of normalized
units, the universal gravitational constant $\cal G$, the mean motion
frequency of the planet $n_{p}$ and the radius of the planet’s orbit
$r_{p}=a_{p}$ are all unitary. In practical simulations, the Sun–Jupiter
system is taken as the fundamental model, in which the length unit is 5.2
$\rm{au}$ and the time unit is 688.955 $\rm{d}$. The normalized mass of the
Sun is $m_{0}=0.9990461188$ and that of the Jupiter is $m_{p}=1-m_{0}$. The
system parameters adopted in this work are the same as the ones used in Lei &
Li (2020), where the first-order prograde MMRs are investigated.

To describe the orbits of the planet and the test particle, we choose the Sun-
centred reference frame with the orbital plane of the planet as the
fundamental plane. In the planar problem, the orbit of the planet is described
by the semimajor axis $a_{p}$ ($a_{p}=r_{p}=1$ in normalized unit) and the
mean longitude $\lambda_{p}$, and the orbit of the test particle is
characterized by the semimajor axis $a$, eccentricity $e$, the longitude of
perihelion $\varpi$ (defined below) and the mean longitude $\lambda=M+\varpi$.

In the Sun-centred reference frame, the motion of the test particle is
governed by the planetary disturbing function, given by (the normalized
variable $r_{p}=1$ is used in the following derivation)

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

where $\Delta$ is the mutual distance between the planet and the test
particle, given by

$\Delta=\sqrt{1+{r^{2}}-2r\cos\psi}$

with $\psi$ as the relative angle between the position vectors of the planet
and the test particle. In the co-planar retrograde configuration (in this
case, the inclination of the test particle’s orbit relative the fundamental
plane is $i=180^{\circ}$), the relative angle $\psi$ can be expressed as

$\psi=f+{\lambda_{p}}-\Omega+\omega=f+{\lambda_{p}}-\varpi,$

where $f$ is the true anomaly and $\varpi$ is defined by
$\varpi=\Omega-\omega$ in retrograde configurations, which can be found in
Shevchenko (2016). In previous studies, $\varpi$ is usually defined by
$\varpi=\omega-\Omega$ in the retrograde case (Morais & Namouni, 2013a, b),
which has an opposite sign in comparison to our definition.

Regarding the definition of $\varpi(=\Omega-\omega)$ adopted in this work,
some discussions are made here. On one hand, in the Sun-centred reference
frame, the eccentricity vector can be expressed by

${\bm{e}}=e\left(\begin{array}[]{l}\cos\Omega\cos\omega-\sin\Omega\sin\omega\cos
i\\\ \sin\Omega\cos\omega+\cos\Omega\sin\omega\cos i\\\ \sin\omega\sin
i\end{array}\right).$

Particularly, in the co-planar (prograde or retrograde) configurations, the
eccentricity vector becomes
${\bm{e}}=e\left(\cos{\varpi},\sin{\varpi},0\right)$ where
$\varpi=\Omega+\omega$ in the prograde co-planar case and
$\varpi=\Omega-\omega$ in the retrograde co-planar case. On the other hand,
when the test particle is located at the perihelion, the celestial longitude
is usually defined by $\alpha=\Omega+\arctan(\tan{\omega}\cos{i})$ (please see
Fig. 1 for the definition of $\alpha$). Naturally, it becomes
$\alpha=\Omega+\omega$ in the co-planar prograde case (i.e., $i=0$) and
$\alpha=\Omega-\omega$ in the co-planar retrograde case (i.e., $i=\pi$). Thus,
in the retrograde occasion, the value of $\varpi(=\Omega-\omega)$ measures the
relative angle between the test particle’s eccentricity vector and the
$x$-axis of the defined coordinate system, showing that
$\varpi(=\Omega-\omega)$ stands for the longitude of perihelion. Naturally,
the definition of $\lambda=M+\varpi(=M+\Omega-\omega)$ is the mean longitude
in the retrograde configuration. Thus, we can find the first advantage of the
definition $\varpi(=\Omega-\omega)$ adopted in this work: both the angles
$\varpi(=\Omega-\omega)$ and $\lambda=M+\varpi(=M+\Omega-\omega)$ have clear
physical implications. In the following discussions, we will provide another
two advantages for the definition $\varpi(=\Omega-\omega)$ used in this work.

Starting from the co-planar and retrograde assumptions (the assumptions have
been used in the expression of $\psi$), we directly follow the procedures
given by Murray & Dermott (1999) and Ellis & Murray (2000) to derive the
expansion of planetary disturbing function in a formal series of the orbital
elements,

$\displaystyle{\cal R}=$ $\displaystyle{\cal
G}{m_{p}}\sum\limits_{n=0}^{N}{\sum\limits_{j=-\infty}^{\infty}{\sum\limits_{m=0}^{n}{\sum\limits_{s=-\infty}^{\infty}{{{\left({-1}\right)}^{n-m}}{A_{n,j}}\left(\alpha\right){n\choose
m}}}}}$ (2) $\displaystyle\times
X_{s}^{m,j}\left(e\right)\cos\left[{s\left(\lambda-\varpi\right)+j\left(\lambda_{p}-\varpi\right)}\right]$
$\displaystyle-{\cal
G}{m_{p}}a\sum\limits_{s=-\infty}^{\infty}{X_{s}^{1,1}\left(e\right)\cos\left[{s\left(\lambda-\varpi\right)+\left(\lambda_{p}-\varpi\right)}\right]}$

where $N$ is the truncated order in the eccentricity (the influence of $N$
upon the accuracy of resonant disturbing function will be reported in Fig. 2).
In equation (2), the Hansen coefficients $X_{s}^{a,b}(e)$ are functions of the
eccentricity $e$ and they can be calculated in a recursive manner (Hughes,
1981; Murray & Dermott, 1999), and ${{A_{n,j}}\left(\alpha\right)}$ is a
function of the semimajor axis ratio $\alpha=a/a_{p}$ ($\alpha=a/a_{p}=a$ in
normalized units), defined by

${A_{n,j}}\left(\alpha\right)=\frac{{{\alpha^{n}}}}{{n!}}\left[{\frac{{{{\rm
d}^{n}}}}{{{\rm
d}{\alpha^{n}}}}\frac{1}{2}b_{{1\mathord{\left/{\vphantom{12}}\right.\kern-1.2pt}2}}^{\left(j\right)}(\alpha)}\right]$

with $b_{1/2}^{j}(\alpha)$ as the Laplace coefficients (Murray & Dermott,
1999; Ellis & Murray, 2000). Since the coefficients $b_{1/2}^{j}(\alpha)$ are
divergent when the semimajor axis ratio is close to unity, the expansion of
disturbing function based on Laplace coefficients is not able to deal with co-
orbital resonances (Murray & Dermott, 1999; Morais & Namouni, 2013b).

It is noted that the literal expansion of disturbing function in the co-planar
retrograde case given by equation (2) has a similar expression (but not the
same) to the one in the co-planar prograde case (please refer to Lei & Li
(2020) for the explicit expansion of disturbing function in the prograde
configuration).

As for the expansion of disturbing function in the retrograde configuration
(orbiting in opposite directions), Morais & Namouni (2013b) derived it
directly from the standard expansion of the three-dimensional and prograde
disturbing function by performing variable substitution:
$i^{*}=180^{\circ}-i$, $\lambda_{p}^{*}=-\lambda_{p}$, $\omega^{*}=\omega-\pi$
and $\Omega^{*}=-\Omega-\pi$. Under the new system of notations, the longitude
of perihelion $\varpi$ is transformed into $\varpi^{*}=\omega-\Omega$ and the
mean longitude $\lambda$ into $\lambda^{*}=M+\omega-\Omega$. As discussed
before, $\varpi^{*}$ and $\lambda^{*}$ are no longer the longitude of
perihelion and mean longitude. The same method utilized to expand the
disturbing function in the retrograde configuration can be found in Li et al.
(2020). Although the expansion presented in this work (see equation (2) for
the explicit expression) and the ones given in previous works (Morais &
Namouni, 2013b; Li et al., 2020) are derived from different approaches, we
believe they should be equivalent.

Observing the explicit expression of disturbing function given by equation
(2), we can see that the argument arising in the cosine terms of disturbing
function holds

$\theta=s(\lambda-\varpi)+j(\lambda_{p}-\varpi),$

which can naturally satisfy the d’Alembert rule (Murray & Dermott, 1999;
Morbidelli, 2002). This is the second advantage for the definition of
$\varpi(=\Omega-\omega)$ adopted in this work.

For convenience, let’s denote the independent angles arising in equation (2)
by

$\theta_{1}=\lambda-\varpi,\quad\theta_{2}=\lambda_{p}-\varpi,$

where $\theta_{1}$ stands for the mean anomaly of the test particle and
$\theta_{2}$ represents the angular separation between the perturber and the
test particle’s perihelion relative to the Sun. As a result, the expansion of
disturbing function given by equation (2) can be organized in a compact form
as follows:

${\cal R}=\sum\limits_{{k_{1}},{k_{2}}}{{{\cal
D}_{{k_{1}},{k_{2}}}}\left({a,e}\right)\cos\left({{k_{1}}{\theta_{1}}+{k_{2}}{\theta_{2}}}\right)},k_{1}\in\mathbb{N},k_{2}\in\mathbb{Z}$
(3)

where ${\cal D}_{{k_{1}},{k_{2}}}$ are the coefficients related to the
semimajor axis $a$ and eccentricity $e$ and the explicit expressions can be
directly derived from equation (2).

Based on the disturbing function, the Hamiltonian function of system, can be
written as

${\cal H}=-\frac{\mu}{{2a}}+{n_{p}}{\Lambda_{p}}-{\cal
R}\left({a,e,\theta_{1},\theta_{2}}\right)$ (4)

where $\mu$ is the gravitational parameter of the Sun, given by $\mu={\cal
G}m_{0}$, the mean motion frequency of the planet is $n_{p}=1$ and
${\Lambda_{p}}$ is the momentum conjugated to the mean longitude of the planet
$\lambda_{p}$.

To formulate the dynamical model, let us adopt the modified Delaunay variables
as follows:

$\displaystyle\Lambda=\sqrt{\mu a},\quad\lambda=M+\varpi,$ (5) $\displaystyle
P=\sqrt{\mu a}\left(1+\sqrt{1-e^{2}}\right),\quad p=-\varpi,$
$\displaystyle\Lambda_{p},\quad\lambda_{p},$

where the pair $({\Lambda_{p}},{\lambda_{p}})$ is introduced to describe the
mean motion of the perturber.

The canonical variables given by equation (5) is different from the
traditional modified Delaunay variables due to the different definition of
$\varpi$ adopted in this work. Please see Morbidelli (2002) for the classical
version of modified Delaunay’s variables.

Using a linear transformation, it is possible to introduce the following set
of variables,

$\displaystyle{\Theta_{1}}=\Lambda,\quad{\theta_{1}}=\lambda+p=\lambda-\varpi,$
(6)
$\displaystyle{\Theta_{2}}=P-\Lambda,\quad{\theta_{2}}={\lambda_{p}}+p={\lambda_{p}}-\varpi,$
$\displaystyle{\Theta_{3}}={\Lambda_{p}}-P+\Lambda,\quad{\theta_{3}}={\lambda_{p}},$

which is a canonical transformation with the generating function,

${\cal
S}=\left({\lambda+p}\right){\Theta_{1}}+\left({{\lambda_{p}}+p}\right){\Theta_{2}}+{\lambda_{p}}{\Theta_{3}}.$

Thus, the Hamiltonian represented by equation (4) can be expressed as a
function of the new set of canonical variables,

$\displaystyle{\cal H}$
$\displaystyle=-\frac{{{\mu^{\rm{2}}}}}{{2\Theta_{1}^{2}}}+{\Theta_{2}}-{\cal
R}\left({{\Theta_{1}},{\Theta_{2}},{\theta_{1}},{\theta_{2}}}\right)$ (7)
$\displaystyle=-\frac{{{\mu^{\rm{2}}}}}{{2\Theta_{1}^{2}}}+{\Theta_{2}}-\sum\limits_{{k_{1}},{k_{2}}}{{{\cal
D}_{{k_{1}},{k_{2}}}}\left({\Theta_{1},\Theta_{2}}\right)\cos\left({{k_{1}}{\theta_{1}}+{k_{2}}{\theta_{2}}}\right)}$

where the constant terms have been removed from the Hamiltonian and the
equality $n_{p}=1$ is used. Evidently, the dynamical model determined by
equation (7) is of two degrees of freedom with $\theta_{1}$ and $\theta_{2}$
as the angular coordinates. The Hamiltonian canonical relations yield the
equations of motion as follows:

$\displaystyle\frac{{{\rm d}{\theta_{1}}}}{{{\rm d}t}}=$
$\displaystyle\frac{{\partial{\cal
H}}}{{\partial{\Theta_{1}}}},\quad\frac{{{\rm d}{\Theta_{1}}}}{{{\rm
d}t}}=-\frac{{\partial{\cal H}}}{{\partial{\theta_{1}}}},$ (8)
$\displaystyle\frac{{{\rm d}{\theta_{2}}}}{{{\rm d}t}}=$
$\displaystyle\frac{{\partial{\cal
H}}}{{\partial{\Theta_{2}}}},\quad\frac{{{\rm d}{\Theta_{2}}}}{{{\rm
d}t}}=-\frac{{\partial{\cal H}}}{{\partial{\theta_{2}}}}.$

For convenience, let’s denote the state of test particles by

${\bm{X}}=\left\\{\theta_{1},\theta_{2},\Theta_{1},\Theta_{2}\right\\}.$ (9)

Thus, the equation of motion given by equation (8) can be expressed in a
vectorial form,

$\dot{\bm{X}}={\bm{F}}\left({\bm{X}}\right).$ (10)

The dynamical model shown by equations (9–10) are to be used in defining and
producing Poincaré surfaces of section, as discussed in Section 5.

##  3 Hamiltonian model of retrograde MMRs

In this section, we formulate the Hamiltonian model of retrograde MMRs by
defining an appropriate resonant angle and performing canonical
transformations. For a retrograde $k_{p}$:$k$ resonance, we define the
critical argument as

$\displaystyle\sigma$ $\displaystyle=\frac{1}{k_{\max}}\varphi$ (11)
$\displaystyle=\frac{1}{k_{\max}}\left[k\lambda-
k_{p}\lambda_{p}+(k_{p}-k)\varpi\right]$
$\displaystyle=\frac{1}{{{k_{\max}}}}\left({k(\lambda-\varpi)-{k_{p}}(\lambda_{p}-\varpi)}\right)$
$\displaystyle=\frac{1}{{{k_{\max}}}}\left({k{\theta_{1}}-{k_{p}}{\theta_{2}}}\right)$

where $\varphi$ is the usual critical argument, given by $\varphi=k\lambda-
k_{p}\lambda_{p}+(k_{p}-k)\varpi$, and
$k_{\max}=\max{\left\\{k_{p},k\right\\}}$ (it holds $k_{\max}=k_{p}$ for the
inner resonances and $k_{\max}=k$ for those exterior resonances). In
particular, we have
$\sigma=\frac{1}{2}\left({{\theta_{1}}-2{\theta_{2}}}\right)$ for the
retrograde 2:1 resonance,
$\sigma=\frac{1}{2}\left({2{\theta_{1}}-{\theta_{2}}}\right)$ for the
retrograde 1:2 resonance, and the like. Note that the resonant angle defined
by equation (11) is equal to the one defined by Lei & Li (2020) for the inner
resonances, but it is different for the outer resonances.

According to the explicit expansion of planetary disturbing function given by
equation (2), the coefficient of the cosine term associated with the resonant
angle, $\cos{(k_{\max}\sigma)}$, holds,

${\cal D}_{0}(\alpha)X_{k}^{m,-k_{p}}(e),$

which stands for the magnitude of the resonant term and it is proportional to
the power of eccentricity. In particular, at the lowest order, the magnitude
of the resonant term is proportional to $e^{k_{p}+k}$, as pointed out by
Morais & Giuppone (2012). For the prograde counterpart, it is known that the
magnitude of the resonant term, at the lowest order, is proportional to
$e^{|k_{p}-k|}$ (Murray & Dermott, 1999). Thus, for a certain $k_{p}$:$k$
resonance with a given eccentricity, the retrograde case has much weaker
resonant strength (or, equivalently, smaller resonant width) compared to the
prograde case, showing that the planetary perturbations for retrograde MMRs
are generally weaker (Morais & Namouni, 2013b). Based on this fact, Morais &
Namouni (2013b) explained that retrograde resonances are more stable than
prograde resonances is due to the weak perturbations in retrograde
configurations where the encounter of two objects occurs at a higher relative
velocity.

Observing the expression of $\varphi$, we can find the usual critical
argument, given by $\varphi=k\lambda-k_{p}\lambda_{p}+(k_{p}-k)\varpi$, keeps
the same expression as that in the prograde configuration (Murray & Dermott,
1999; Morbidelli, 2002). This is the third advantage for the definition of
$\varpi(=\Omega-\omega)$ adopted in this work. Based on this advantage, it is
possible to unify the resonant models of prograde and retrograde resonances
from theoretical point of view (this point is very important to understand the
retrograde resonances).

To formulate the resonant model, let us further introduce the following
canonical transformations:

$\displaystyle{\Gamma_{1}}=\frac{k_{\max}}{k}\Theta_{1},\quad{\sigma_{1}}=\frac{1}{{{k_{\max}}}}\left(k\theta_{1}-k_{p}\theta_{2}\right)=\sigma,$
(12)
$\displaystyle{\Gamma_{2}}=\Theta_{2}+\frac{k_{p}}{k}\Theta_{1},\quad{\sigma_{2}}=\theta_{2}$

with the generating function as

${\cal
S}=\frac{\Gamma_{1}}{k_{\max}}\left(k\theta_{1}-k_{p}\theta_{2}\right)+\theta_{2}\Gamma_{2}.$

Using the defined canonical variables, the Hamiltonian given by equation (7)
can be written as

$\displaystyle{{\cal H}}=$
$\displaystyle-\frac{{{\mu^{2}}}}{{2{{\left({\frac{k}{{{k_{\max}}}}{\Gamma_{1}}}\right)}^{2}}}}-\frac{{{k_{p}}}}{{{k_{\max}}}}{\Gamma_{1}}$
(13) $\displaystyle-\sum\limits_{{k_{1}},{k_{2}}}{{\cal
D}_{{k_{1}},{k_{2}}}}\left({\Gamma_{1},\Gamma_{2}}\right)\cos\left[{k_{\max}}\frac{k_{1}}{k}{\sigma_{1}}+\left(k_{2}+\frac{k_{1}}{k}\right){\sigma_{2}}\right].$

When the test particle is located inside a retrograde MMR, the resonant angle
$\sigma_{1}(=\sigma)$ becomes a ‘long-period’ variable, while the angle
$\sigma_{2}$ is a ‘short-period’ variable. Thus, in the expansion of
disturbing function those terms involving the angle $\sigma_{2}$ have short-
period effects upon the evolution of test particles. In order to study
resonant dynamics, it is usual to filter out those short-period terms from the
Hamiltonian function, producing the so-called resonant Hamiltonian as follows:

$\displaystyle{{\cal H}^{*}}=$
$\displaystyle\frac{1}{{2k\pi}}\int\limits_{0}^{2k\pi}{{\cal
H}\left({\Gamma_{1},\Gamma_{2},{\sigma_{1}},{\sigma_{2}}}\right){\rm
d}{\sigma_{2}}}$ (14) $\displaystyle=$
$\displaystyle-\frac{{{\mu^{2}}}}{{2{{\left({\frac{k}{{{k_{\max}}}}{\Gamma_{1}}}\right)}^{2}}}}-\frac{{{k_{p}}}}{{{k_{\max}}}}{\Gamma_{1}}$
$\displaystyle-\frac{1}{{2k\pi}}\int\limits_{0}^{2k\pi}{{\cal
R}\left({\Gamma_{1},\Gamma_{2},{\sigma_{1}},{\sigma_{2}}}\right){\rm
d}{\sigma_{2}}}$

which becomes

${{\cal
H}^{*}}=-\frac{{{\mu^{2}}}}{{2{{\left({\frac{k}{{{k_{\max}}}}{\Gamma_{1}}}\right)}^{2}}}}-\frac{{{k_{p}}}}{{{k_{\max}}}}{\Gamma_{1}}-\sum\limits_{m=0}^{M}{{{\cal
C}_{m}}\cos\left({m{k_{\max}}{\sigma_{1}}}\right)}$ (15)

where ${{\cal C}_{m}}$ are the coefficients related to the action variables
$\Gamma_{1}$ and $\Gamma_{2}$ and the number $M$ stands for the number of
harmonics of the angle $k_{\max}\sigma_{1}$ in the expansion of resonant
disturbing function. The last term of equation (15) corresponds to the
resonant disturbing function, whose explicit expression is provided in
Appendix A.

Figure 2: Comparisons between the resonant Hamiltonian obtained by means of
numerical-averaging technique and series expansions truncated at orders $N=2$,
$N=4$ and $N=6$ for the retrograde 2:1 resonance. In simulations, the
semimajor axis is taken as the nominal value of resonance (i.e.,
$a_{0}=0.62976016$ in normalized units), the value of eccentricity is fixed at
$e=0.3$ and the number of harmonics is assumed at $M=2$. For the retrograde
2:1 resonance, the resonant angle is defined by
$\sigma=\frac{1}{2}(\theta_{1}-2\theta_{2})=\frac{1}{2}\varphi$.

Alternatively, the exact value of resonant disturbing function arising in
equation (14) can be obtained by means of numerical-averaging approach (i.e.,
the unexpanded disturbing function is numerically averaged with respect to the
fast angle $\sigma_{2}$ over $k$ periods), which has been widely used in
previous works (Beaugé, 1994; Gallardo, 2006b, a; Morais & Namouni, 2013b; Li
et al., 2014a, b, 2020; Huang et al., 2018; Gallardo, 2019b, 2020; Li et al.,
2020). Please refer to Fig. 2 for the comparisons between the resonant
Hamiltonian obtained by means of numerical-averaging technique and series
expansions truncated at orders $N=2$, $N=4$ and $N=6$ for the retrograde 2:1
resonance with $e=0.3$ (the semimajor axis is taken as the nominal value of
resonance at $a_{0}=0.62976016$ in normalized units). It shows that the curve
corresponding to a larger $N$ is closer to the reference curve obtained by
means of numerical averaging, indicating that the series expansion truncated
at a higher order has a better accuracy. This is expected by us.

The resonant Hamiltonian given by equation (15) determines a dynamical model
with a single degree of freedom. Since the angle $\sigma_{2}$ is absent from
the Hamiltonian, its conjugate momentum becomes the motion integral in the
resonant model, given by

$\displaystyle{\Gamma_{2}}$
$\displaystyle=\Theta_{2}+\frac{k_{p}}{k}\Theta_{1}$ (16)
$\displaystyle=\sqrt{\mu
a}\left({\frac{{{k_{p}}}}{k}+\sqrt{1-{e^{2}}}}\right)={\rm const},$

which means that, in the resonant model, the semimajor axis and eccentricity
exchange with each other in the long-term evolution. For convenience, the
motion integral $\Gamma_{2}$ can be specified by the minimum semimajor axis
$a_{\min}$ in the following manner:

${\Gamma_{2}}=\sqrt{\mu{a_{\min}}}\left({\frac{{{k_{p}}}}{k}+1}\right),$

where $a_{\min}$ is the magnitude of semimajor axis when the eccentricity is
assumed at zero.

For the retrograde 2:1 and 1:2 resonances, the level curves of the motion
integral $\Gamma_{2}$ are plotted in the $(a,e)$ space, as shown in Fig. 3. In
particular, two representatives with $\Gamma_{2}=2.34$ (or
$a_{\min}=0.608981$) and $\Gamma_{2}=1.85$ (or $a_{\min}=1.522563$) are marked
by red lines.

Figure 3: Level curves of the motion integral ${\Gamma_{2}}$ shown in the
$(a,e)$ space for the retrograde 2:1 resonance (_upper panel_) and the
retrograde 1:2 resonance (_bottom panel_). The motion integral is defined by
$\Gamma_{2}=\sqrt{\mu a}(\frac{k_{p}}{k}+\sqrt{1-e^{2}})$, which can be
equivalently represented by $a_{\min}$ in the form of $\Gamma_{2}=\sqrt{\mu
a_{\min}}(\frac{k_{p}}{k}+1)$ ($a_{\min}$ is the semimajor axis when the
eccentricity is assumed at zero). The red lines marked in both panels
correspond to $\Gamma_{2}=2.34$ (or $a_{\min}=0.608981$) for the retrograde
2:1 resonance and $\Gamma_{2}=1.85$ (or $a_{\min}=1.522563$) for the
retrograde 1:2 resonance, which are to be used as typical examples in the
following simulations. The nominal location of resonance is marked by a short
and black vertical line. Figure 4: Schematic diagram for illustrating elliptic
points (i.e., resonant centres), saddle points, dynamical separatrices, curves
of libration and circulation, and resonant width measured from the phase
portrait in the $(\sigma_{1},\Gamma_{1})$ space. Here, we take the 2:1
resonance as an example to show the diagram.

##  4 Analytical study

In the previous section, we have formulated the Hamiltonian model for
retrograde MMRs by introducing the resonant angle $\sigma(=\sigma_{1})$ and by
performing a series of canonical transformations. The resulting resonant model
is of a single degree of freedom, and thus it is integrable. In this section,
we will introduce the analytical method on how to determine the resonant
centre and the associated resonant width, and then apply it to some retrograde
MMRs (in practice, the retrograde 2:1 and 1:2 resonances are taken as examples
in simulations).

###  4.1 Resonant centre and width

In the resonant model (remind that in the resonant model $\Gamma_{2}$ is a
motion integral), the equilibrium points should satisfy the following
stationary conditions:

$\displaystyle{{\dot{\sigma}}_{1}}$ $\displaystyle=\frac{{\partial{{\cal
H}^{*}}}}{{\partial{\Gamma_{1}}}}=0,$ (17) $\displaystyle{{\dot{\Gamma}}_{1}}$
$\displaystyle=-\frac{{\partial{{\cal
H}^{*}}}}{{\partial{\sigma_{1}}}}=\sum\limits_{m=1}^{M}{m{k_{\max}}{{\cal
C}_{m}}\sin\left({m{k_{\max}}{\sigma_{1}}}\right)}=0.$

The second condition of equation (17) shows that the equilibrium points are
located at

${\sigma_{1}}=\frac{q\pi}{k_{\max}},\quad q\in\mathbb{Z}.$

Replacing the resonant Hamiltonian $\cal H^{*}$ in the first condition of
equation (17), we can obtain the equilibrium points satisfy

$k{\mu^{2}}{\left({\frac{k}{{{k_{\max}}}}{\Gamma_{1}}}\right)^{-3}}={k_{p}}+{k_{\max}}\sum\limits_{m=0}^{M}{\frac{{\partial{{\cal
C}_{m}}}}{{\partial{\Gamma_{1}}}}\cos\left({mq\pi}\right)}$ (18)

By numerically solving this equation, the positions of equilibrium points,
denoted by $(\sigma_{1}=\frac{q\pi}{k_{\max}},\Gamma_{1})$, can be determined.

In particular, if the second term in the right-hand side of equation (18) is
ignored, it is possible to obtain the nominal resonant centre at

$\Gamma_{1}=\left(\frac{\mu^{2}k}{k_{p}}\right)^{1/3}\frac{k_{\max}}{k}$

which is

$a={\mu}^{1/3}(\frac{k}{k_{p}})^{2/3}.$

Thus, the deviation of the exact location of equilibrium points from the
nominal resonance centre is due to the rightmost term in equation (18).

For a given motion integral $\Gamma_{2}$, let us denote the stable equilibrium
points (i.e., elliptic points or resonant centres) by
$({\sigma_{1,s}},{\Gamma_{1,s}})$ and the unstable equilibrium points (i.e.,
saddle points) as $(\sigma_{1,u},\Gamma_{1,u})$. According to the expression
of $\Gamma_{1}(=\Gamma_{1,s})$ and $\Gamma_{2}$, we can identify the location
of resonant centres in the $(a,e)$ space, denoted by $(a_{0},e_{0})$. As
discussed by (Morbidelli, 2002), the level curves of resonant Hamiltonian
passing through saddle points play the role of dynamical separatrices,
dividing the total phase space into regions of libration and circulation. In
addition, the distance between adjacent separatrices evaluated at the resonant
centre stands for the resonant width (please refer to Fig. 4 for the
definition of resonant width).

Let’s denote the states on the separatrices evaluated at the resonant centre
with $\sigma_{1}=\sigma_{1,s}$ by $(\sigma_{1,s},\Gamma_{1}^{L})$ and
$(\sigma_{1,s},\Gamma_{1}^{R})$ (see Fig. 4 for details). Here, the
separatrices stem from the saddle point at $(\sigma_{1,u},\Gamma_{1,u})$.
According to the definition of resonant width, we can obtain the following
relation:

${\cal H}^{*}(\Gamma_{2};\sigma_{1,u},\Gamma_{1,u})={\cal
H}^{*}(\Gamma_{2};\sigma_{1,s},\Gamma_{1}^{L})={\cal
H}^{*}(\Gamma_{2};\sigma_{1,s},\Gamma_{1}^{R}).$ (19)

For a given motion integral $\Gamma_{2}$, it is possible to obtain the action
variables $\Gamma_{1}^{L}$ and $\Gamma_{1}^{R}$ on the separatrices by
numerically solving equation (19). Naturally, the resonant width can be
represented by $\Delta\Gamma_{1}=\Gamma_{1}^{R}-\Gamma_{1}^{L}$.

In our study, we would like to introduce the variations of semimajor axis and
eccentricity ($\Delta a$ and $\Delta e$) to measure the resonant width. To
this end, we need to determine the points on the separatrices, shown by
$(a_{L},e_{L})$ and $(a_{R},e_{R})$. According to the relationship between
$\Gamma_{1}$ and the semimajor axis $a$,

$\Gamma_{1}=\frac{k_{\max}}{k}\sqrt{\mu a},$

it is possible to determine the semimajor axes of the points at the
boundaries, denoted by $a_{L}$ and $a_{R}$. Furthermore, the conservation of
motion integral $\Gamma_{2}$ leads to the following relations (Lei & Li,
2020):

$\displaystyle\Gamma_{2}$ $\displaystyle=\sqrt{\mu
a_{0}}\left({\frac{{{k_{p}}}}{k}+\sqrt{1-{e_{0}^{2}}}}\right)$ (20)
$\displaystyle=\sqrt{\mu
a_{L}}\left({\frac{{{k_{p}}}}{k}+\sqrt{1-{e_{L}^{2}}}}\right)$
$\displaystyle=\sqrt{\mu
a_{R}}\left({\frac{{{k_{p}}}}{k}+\sqrt{1-{e_{R}^{2}}}}\right)$

where $(a_{0},e_{0})$ stands for the position of the resonant centre. Solving
equation (20), we can identify the points on the boundaries in the $(a,e)$
space as $(a_{L},e_{L})$ and $(a_{R},e_{R})$. Thus, the resonant width in
terms of the variations of semimajor axis and eccentricity can be expressed by

$\Delta a=a_{R}-a_{L},\quad\Delta e=e_{R}-e_{L}.$

###  4.2 Analytical results

In the resonant model, there is a motion integral, denoted by $\Gamma_{2}$,
which can be specified by $a_{\min}$. For such an integral system, the global
dynamics in the phase space can be revealed by phase portraits (i.e., level
curves of resonant Hamiltonian) with a given motion integral $\Gamma_{2}$ (or
$a_{\min}$). For intuition, in this work we show the phase portraits in the
$(e\cos{\sigma},e\sin{\sigma})$ and $(\sigma,a)$ spaces.

Figure 5: The level curves of resonant Hamiltonian (i.e., phase portraits)
shown in the $(e\cos{\sigma},e\sin{\sigma})$ space (_left panel_) and in the
$(\sigma,a)$ space (_right panel_) for the retrograde 2:1 resonance specified
by $a_{\min}=0.608981$ (or $\Gamma_{2}=2.34$). For the considered resonance,
the critical argument is defined as
$\sigma=\frac{1}{2}\left(\lambda-2\lambda_{p}+\varpi\right)=\frac{1}{2}\left(\theta_{1}-2\theta_{2}\right)$.
The resonant centres are located at $\sigma=0$ and $\sigma=\pi$ (corresponding
to the usual argument at $\varphi=0$), and the saddle points are located at
$\sigma=\pm\pi/2$ (corresponding to $\varphi=\pi$). The level curves passing
through the saddle points are marked in red lines, which play the role of
dynamical separatrices, dividing the phase space into libration and
circulation regions. The distance between two nearby separatrices, evaluated
at the resonant centre, measures the resonant width. In particular, the
resonant width in terms of the variation of eccentricity $\Delta e$ is marked
in the _left panel_ , and the one in terms of the variation of semimajor axis
$\Delta a$ is marked in the _right panel_.

Figure 5 shows the phase portraits of the retrograde 2:1 resonance specified
by $a_{\min}=0.608981$ (or $\Gamma_{2}=2.34$). The isoline of
$\Gamma_{2}=2.34$ is shown in the upper panel of Fig. 3 and, for the
retrograde 2:1 resonance, the resonant angle is defined by
$\sigma=\frac{1}{2}\left(\theta_{1}-2\theta_{2}\right)=\frac{1}{2}\left(\lambda-2\lambda_{p}+\varpi\right)=\frac{1}{2}\varphi$.
From Fig. 5, it is observed that (a) the resonant centres are located at
$\sigma=0$ and $\sigma=\pi$ (corresponding to the usual argument at
$\varphi=0$), (b) the saddle points are placed at $\sigma=\pi/2$ and
$\sigma=-\pi/2$ (corresponding to $\varphi=\pi$), (c) the level curves passing
through saddle points (marked in red lines) play the role of dynamical
separatrices, which divide the entire phase space into libration and
circulation regions, and (d) the distance between neighboring separatrices
evaluated at the resonant centre can measure the resonant width and, in
particular, the resonant width in terms of the variation of eccentricity
$\Delta e$ is marked in the _left panel_ and the one in terms of the variation
of semimajor axis $\Delta a$ is marked in the _right panel_.

In the polar coordinate plane (see the _left panel_ of Fig. 5), the coordinate
origin with zero eccentricity is no longer a stationary point. This is
different from the prograde case (it is known from Malhotra & Zhang (2020) and
Lei & Li (2020) that in the prograde case the zero-eccentricity point is
always a saddle point in the resonant model).

Figure 6: Same as Fig. 5 but for the retrograde 1:2 resonance specified by
$a_{\min}=1.522563$ (or $\Gamma_{2}=1.85$). For the considered resonance, the
critical argument is defined as
$\sigma=\frac{1}{2}\left(2\lambda-\lambda_{p}-\varpi\right)=\frac{1}{2}\left(2\theta_{1}-\theta_{2}\right)$.
The resonant centres are located at $\sigma=\pi/2$ and $\sigma=3\pi/2$
(corresponding to $\varphi=\pi$), and the saddle points are located at
$\sigma=0$ and $\sigma=\pi$ (corresponding to $\varphi=0$).

The phase portraits of the retrograde 1:2 resonance specified by
$a_{\min}=1.522563$ (or $\Gamma_{2}=1.85$) are shown in Fig. 6, where the
level curves passing through saddle points are marked in red lines. For the
considered resonance, the isoline of $\Gamma_{2}=1.85$ has been marked in the
bottom panel of Fig. 3 and the resonant angle is defined by
$\sigma=\frac{1}{2}\left(2\theta_{1}-\theta_{2}\right)=\frac{1}{2}\left(2\lambda-\lambda_{p}-\varpi\right)=\frac{1}{2}\varphi$.
From Fig. 6, it is observed that the resonant centres are located at
$\sigma=\pi/2$ and $\sigma=3\pi/2$ (corresponding to the usual argument at
$\varphi=\pi$) and the saddle points are located at $\sigma=0$ and
$\sigma=\pi$ (corresponding to the usual argument at $\varphi=0$). The
distance between neighboring separatrices evaluated at the resonant centre
also stands for the resonant width and, in particular, the resonant width in
terms of $\Delta a$ and $\Delta e$ is indicated in the phase portraits.

It is known that the prograde 1:$n$ resonances hold asymmetric libration
centres with the usual critical argument $\varphi(=k_{\max}\sigma)$ different
from zero or $\pi$ (Beaugé, 1994; Morbidelli, 2002). The appearance of the
asymmetric libration centres is because the second-order harmonics dominates
the resonant Hamiltonian (Morbidelli, 2002). However, there is no asymmetric
libration centres in the phase portraits of the retrograde 1:2 resonance. The
difference about the symmetric properties of libration centres between the
prograde and retrograde 1:$n$ resonances has been noticed by Li et al. (2020).

From Figs 5 and 6, we can also see that the phase portraits shown in the
$(e\cos{\sigma},e\sin{\sigma})$ space are symmetric with respect to the lines
corresponding to $\sin{\sigma}=0$ and $\cos{\sigma}=0$. These symmetric
structures indicate that the resonant centres in the analytical model have the
same dynamics. Thus, in the following discussions, we take one of them into
consideration to discuss the dynamics.

Figure 7: Analytical results about the distributions of resonant centre
$(a_{0},e_{0})$ and resonant width $(\Delta a=a_{R}-a_{L},\Delta
e=e_{R}-e_{L})$ for the retrograde 2:1 resonance (_left panel_) and for the
retrograde 1:2 resonance (_right panel_). For convenience, the level curves of
the motion integral $\Gamma_{2}$ are plotted and, in particular, the isolines
of $\Gamma_{2}=2.34$ and $\Gamma_{2}=1.85$ are marked in blue dotted lines.
The red lines stand for the location of resonant centres, and the shaded areas
represent the libration regions. Note that the resonant width is measured
along the isoline of the motion integral $\Gamma_{2}$, as shown in the _left
panel_. The nominal location of resonance is marked by a short and black
vertical line at the bottom of each figure.

According to the method presented in Section 4.1, the location of resonant
centre and the boundaries that separate the libration regions from circulation
regions can be analytically determined. In Fig. 7, the resonant centre and
resonant width are reported in the $(a,e)$ space for the retrograde 2:1
resonance in the _left panel_ and for the retrograde 1:2 resonance in the
_right panel_. For convenience, in Fig. 7, the level curves of the motion
integral $\Gamma_{2}$ are also provided and the nominal location of resonance
is marked by a short and black vertical line. The curve of $(a_{0},e_{0})$
stands for the location of resonant centre, the curve of $(a_{L},e_{L})$ for
the left separatrix, and the curve of $(a_{R},e_{R})$ for the right
separatrix. The shaded areas bounded by the left and right separatrices
represent the libration regions in the $(a,e)$ space, and the resonant width
is measured along the isoline of $\Gamma_{2}$. It is observed from Fig. 7
that, in both cases, (a) the left and right separatrices are not symmetric
with respect to the nominal resonance location, (b) the curve of resonant
centre cannot extend to the extremely low-eccentricity region, (c) the
resonant centre deviates from the nominal location of resonance and (d) the
resonant width ($\Delta a$ or $\Delta e$) is an increasing function of the
eccentricity. In particular, for the retrograde inner 2:1 resonance, the
resonant centres are on the left-hand side of the nominal location of
resonance and, for the retrograde outer 1:2 resonance, the resonant centres
are on the right-hand side (please refer to the short and vertical black lines
in Fig. 7 for the nominal locations of resonance).

##  5 Numerical study

To validate the analytical results (including the phase portraits and resonant
width), in this section we turn to explore the dynamical structures of
retrograde MMRs by analyzing Poincaré sections.

###  5.1 Poincaré surfaces of section

As discussed in Section 2, the full model represented by equation (8) is of
two degrees of freedom with $\theta_{1}$ and $\theta_{2}$ as angular
coordinates. For such a two-degree of freedom system, Poincaré surface of
section is an efficient technique to explore its global dynamics in the phase
space.

Figure 8: Schematic diagrams for the definition of Poincaré sections for the
inner and outer cases of retrograde MMRs with a given motion integral
$\Gamma_{2}$. In both diagrams, the central body with mass $m_{0}$ stands for
the Sun and the perturber with mass $m_{p}$ stands for the Jupiter-mass
planet. It is known that, in the planar retrograde problem, there are two
independent angular variables including $\theta_{1}=\lambda-\varpi$ and
$\theta_{2}=\lambda_{p}-\varpi$. For the inner case, the angular variable
$\theta_{1}=\lambda-\varpi$ is faster than $\theta_{2}=\lambda_{p}-\varpi$, so
that the section is defined by $\theta_{1}=\lambda-\varpi=0$ (i.e., the
particle is at its perihelion), where the angular separation of the perturber
from the test particle’s perihelion is ${\theta_{2}}=\lambda_{p}-\varpi$. For
the outer case, the angular variable $\theta_{2}=\lambda_{p}-\varpi$ is faster
than $\theta_{1}=\lambda-\varpi$, so that the section is defined by
$\theta_{2}=\lambda_{p}-\varpi=0$ (i.e., the perturber is on the line
directing from the central body towards the particle’s perihelion), where the
angular separation between the test particle and its perihelion is
$\theta_{1}=\lambda-\varpi$. It is noted that the definition in the inner case
is similar to the one defined by Malhotra & Zhang (2020), but the one in the
outer case is different.

In this part, we adopt the numerical technique based on Poincaré sections to
explore the dynamical structures of retrograde MMRs. In particular, we make
some modifications for producing Poincaré sections in comparison to the
traditional versions used by Malhotra (1996); Winter & Murray (1997a, b);
Morais & Giuppone (2012); Morais & Namouni (2013b); Wang & Malhotra (2017);
Malhotra et al. (2018) and Malhotra & Zhang (2020) in terms of the following
two points.

  * •

A modified condition is adopted for computing Poincaré sections. It is known
from the Hamiltonian given by equation (7) that, in the non-averaged model,
the particles move in a four-dimensional phase space
$(a,e,\theta_{1},\theta_{2})$, where $\theta_{1}=\lambda-\varpi$ and
$\theta_{2}=\lambda_{p}-\varpi$. In general, for the inner resonances (i.e.,
$a<a_{p}$), the angle $\theta_{1}$ is faster than $\theta_{2}$ and, for those
outer resonances (i.e., $a>a_{p}$), the angle $\theta_{2}$ is faster than
$\theta_{1}$. Based on this, we produce Poincaré sections by recording those
states of test particles when the short-period variable ($\theta_{1}$ for the
inner case and $\theta_{2}$ for the outer case) is equal to zero, as
illustrated in Fig. 8. In particular, for the inner case, the section is
defined by $\theta_{1}=\lambda-\varpi=0$, meaning that the test particle is
located at its perihelion (see the _left panel_ of Fig. 8) and, for the outer
case, the section is defined by $\theta_{2}=\lambda_{p}-\varpi=0$, where the
perturber is in the same direction of test particle’s eccentricity vector (see
the _right panel_ of Fig. 8). It is noted that the choice of Poincaré section
in the inner case is similar to that taken by Wang & Malhotra (2017); Malhotra
et al. (2018) and Malhotra & Zhang (2020), but the one in the outer case is
different.

  * •

The motion integral $\Gamma_{2}$ (rather than the approximate Jacobi constant)
is taken as the conserved quantity. In the production of Poincaré sections
adopted by Wang & Malhotra (2017); Malhotra et al. (2018) and Malhotra & Zhang
(2020), the approximate Jacobi constant, given by

$C_{J}\approx\frac{\mu}{a}+2\sqrt{\mu a(1-e^{2})},$

is taken as the conserved quantity. In this work, in order to compare with the
analytical developments presented in Section 3, we adopt the motion integral,
given by

$\Gamma_{2}=\sqrt{\mu a}\left(\frac{k_{p}}{k}+\sqrt{1-e^{2}}\right),$

as the conserved quantity in producing Poincaré sections. It is noted that the
motion integral $\Gamma_{2}$ is an approximate constant in the non-averaged
model and, in the vicinity of resonance, the periodic variation is on the
order of the mass ratio between the perturber and the Sun. For a given motion
integral $\Gamma_{2}$, the four-dimensional phase space
$(a,e,\theta_{1},\theta_{2})$ is restrained to a three-dimensional space,
i.e., $(a,\theta_{1},\theta_{2})$ or $(e,\theta_{1},\theta_{2})$. Furthermore,
due to the choice of Poincaré section, the three-dimensional space is further
reduced to a two-dimensional phase space, i.e., $(a,\theta_{2})$ or
$(e,\theta_{2})$ for the inner case and $(a,\theta_{1})$ or $(e,\theta_{1})$
for the outer case. Consequently, the Poincaré sections can be shown in a two-
dimensional parameter space.

For the inner case, we can further observe that, in magnitude, the angular
separation $\theta_{2}=\lambda_{p}-\varpi$ is equal to $\lambda_{p}-\lambda$
because it holds $\theta_{1}=\lambda-\varpi=0$ on the Poincaré sections.
Similarly, for the outer case we can see that the angle
$\theta_{1}=\lambda-\varpi$ is equal to $\lambda-\lambda_{p}$ because it holds
$\theta_{2}=\lambda_{p}-\varpi=0$ on the Poincaré sections. In other words,
for both cases, the angular separation ($\theta_{2}$ for the outer case and
$\theta_{1}$ for the inner case) on the Poincaré section stands for the
synodic angle between the test particle and the perturber.

In the production of Poincaré sections, the equations of motion of the non-
averaged model (i.e., the planar circular restricted three-body problem) are
numerically integrated over a long enough time and the states are recorded at
every time when test particles pass through the defined sections. The
numerical integrator we used in this study is an eighth-order Runge–Kutta
algorithm with step-size control at the seventh order (Fehlberg, 1968), and
the relative and absolute error tolerances are controlled to be smaller than
$1.0\times 10^{-12}$.

Next, let us consider the relationship between the resonant angle $\sigma$
defined in Section 3 and the angular separation $\theta_{2}$ (or $\theta_{1}$)
for the inner (or outer) resonances. For the retrograde inner resonances (in
this case, it holds $k_{p}>k$), the resonant angle is given by
$\sigma=\frac{1}{k_{p}}\left(k\theta_{1}-k_{p}\theta_{2}\right)$, which is
equal to $-\theta_{2}$ on the Poincaré section defined by $\theta_{1}=0$ (see
the _left panel_ of Fig. 8). Similarly, for the retrograde outer resonances
(in this case, it holds $k_{p}<k$), the resonant angle
$\sigma=\frac{1}{k}\left(k\theta_{1}-k_{p}\theta_{2}\right)$ is equal to
$\theta_{1}$ on the Poincaré section defined by $\theta_{2}=0$ (see the _right
panel_ of Fig. 8). Thus, for a given motion integral $\Gamma_{2}$, it is
possible for us to directly compare the phase portraits obtained in the
previous section with the associated Poincaré sections for the retrograde
$k_{p}$:$k$ resonance. This comparison is very helpful to understand the
structures arising in the Poincaré sections and also to validate the
analytical developments formulated in Sections 3 and 4.

Figure 9: Poincaré surfaces of section defined by
$\theta_{1}=\lambda-\varpi=0$ for the retrograde 2:1 resonance specified by
$a_{\min}=0.608981$ (or $\Gamma_{2}=2.34$). At the Poincaré section, the
critical argument
$\sigma=\frac{1}{2}\left(\lambda-2\lambda_{p}+\varpi\right)=\frac{1}{2}\left(\theta_{1}-2\theta_{2}\right)$
becomes $\sigma=\varpi-\lambda_{p}=-\theta_{2}$ (because of
$\theta_{1}=0^{\circ}$ on the section), measuring the angular separation
between the test particle and the perturber when the former is located at its
perihelion (i.e., $\theta_{1}=\lambda-\varpi=0$). In the sections, the curves
separating libration regions from circulation regions are marked by red lines.
The resonant width in terms of variations of semimajor axis and eccentricity
($\Delta a$ and $\Delta e$) is explicitly marked.

Figure 9 presents the Poincaré sections defined by
$\theta_{1}=\lambda-\varpi=0$ in the $(e\cos{\theta_{2}},e\sin{\theta_{2}})$
and $(\theta_{2},a)$ spaces for the retrograde 2:1 resonance specified by
$a_{\min}=0.608981$ (or $\Gamma_{2}=2.34$). The motion integral adopted here
is equal to that used in Fig. 5. In the Poincaré sections, it is known that
the scattered points stands for the chaotic motion and the smooth curves
represent regular orbits in the full model. In particular, inside the island,
the motions are of libration and, outside the island, the motions are of
circulation. From Fig. 9, it is observed that, in the considered phase space,
all the orbits are regular (i.e., no chaotic motions are found in the
sections) and there are two islands of libration centred at $\theta_{2}=\pi/2$
and $\theta_{2}=-\pi/2$ and two saddle points located at $\theta_{2}=0$ and
$\theta_{2}=\pi$. According to the relationship between the resonant angle
$\sigma$ and the angular separation $\theta_{2}$, we have $\theta_{2}=-\sigma$
on the sections defined by $\theta_{1}=0$. Thus, we can get that the islands
arising in the Poincaré sections are centred at $\sigma=-\pi/2$ and
$\sigma=\pi/2$ (corresponding to the usual critical argument at $\varphi=\pi$)
and the saddle points are located at $\sigma=0$ and $\sigma=\pi$
(corresponding to $\varphi=0$). The red curves bounding the islands of
libration play the role of dynamical separatrix, which divides the entire
phase space into regions of libration and circulation. In particular, in the
area bounded by the separatrix, the motions are of libration (i.e., they are
quasi-periodic orbits in the full model). In the Poincaré sections, the size
of libration island can be measured by the distance between two neighboring
separatrices evaluated at the centre of libration island (Malhotra & Zhang,
2020), as shown by $\Delta e$ in the _left panel_ and $\Delta a$ in the _right
panel_.

Comparing the Poincaré sections shown in Fig. 9 and the phase portraits shown
in Fig. 5 (both of them have the same motion integral at $\Gamma_{2}=2.34$),
we can observe an excellent correspondence between the structures appearing in
the Poincaré sections and the ones arising in the phase portraits: (a) the
resonant centres and saddle points in both the Poincaré sections and phase
portraits are in perfect agreement, and (b) the resonant width, measuring the
size of libration island, can be determined by analyzing both the Poincaré
sections and phase portraits (they are called analytical and numerical widths
of resonance, respectively).

Figure 10: Poincaré surfaces of section defined by
$\theta_{2}=\lambda_{p}-\varpi=0$ for the retrograde 1:2 resonance specified
by $a_{\min}=1.522563$ (or $\Gamma_{2}=1.85$). Regarding the points on the
Poincaré section, the critical argument
$\sigma=\frac{1}{2}\left(2\lambda-\lambda_{p}-\varpi\right)=\frac{1}{2}\left(2\theta_{1}-\theta_{2}\right)$
becomes $\sigma=\lambda-\varpi=\theta_{1}$ (because of $\theta_{2}=0^{\circ}$
on the section), measuring the angular separation between the test particle
and its perihelion when the perturber is in the same direction of the
particle’s perihelion (i.e., $\theta_{2}=\lambda_{p}-\varpi=0$). Similar to
Fig. 9, the curves separating the libration regions from circulation regions
are marked by red lines, and the resonant width in terms of the variations of
semimajor axis and eccentricity ($\Delta a$ and $\Delta e$) is marked.

Similarly, for the retrograde 1:2 resonance specified by $a_{\min}=1.522563$
(or $\Gamma_{2}=1.85$), the Poincaré surfaces of section are reported in Fig.
10. The motion integral adopted here is equal to that taken in Fig. 6. For the
considered retrograde outer resonance, we have $\sigma=\theta_{1}$ on the
sections defined by $\theta_{2}=0$. Observing Fig. 10, we can see that (a) the
islands of libration are centred at $\theta_{1}=\sigma=\pi/2$ and
$\theta_{1}=\sigma=3\pi/2$ (corresponding to $\varphi=\pi$) and (b) the saddle
points are located at $\theta_{1}=\sigma=0$ and $\theta_{1}=\sigma=\pi$
(corresponding to $\varphi=0$). The curves passing through the saddle points
are marked in red lines and the size of libration island is measured by
$\Delta e$ and $\Delta a$, as shown by the left and right panels,
respectively. Comparing Fig. 10 with Fig. 6 leads us to the conclusion that,
for the retrograde 1:2 resonance, the Poincaré sections and phase portraits
are in good agreement in terms of the dynamical structures arising in the
phase space.

###  5.2 Numerical results

As discussed in the previous subsection, the dynamical separatrices can be
found in the Poincaré sections. Thus, it is possible for us to evaluate their
distance at the centre of libration island, which stands for the numerical
width of resonance. As shown in Figs 9 and 10, the width of resonance is
described by the variations of semimajor axis and eccentricity, denoted by
$\Delta a$ and $\Delta e$. Note that the determination of resonant width by
analyzing Poincaré sections can be found in a series of publications, e.g.
Malhotra (1996); Winter & Murray (1997a); Wang & Malhotra (2017); Lan &
Malhotra (2019); Malhotra & Zhang (2020).

Figure 11: Numerical results about the resonant centre $(a_{0},e_{0})$ and
resonant width $(\Delta a=a_{R}-a_{L},\Delta e=e_{R}-e_{L})$ produced by
analyzing the Poincaré sections for the retrograde 2:1 resonance (_left
panel_) and the retrograde 1:2 resonance (_right panel_). For convenience, the
level curves of the motion integral $\Gamma_{2}$ are still plotted and, in
particular, the isolines of $\Gamma_{2}=2.34$ and $\Gamma_{2}=1.85$ are marked
in blue dotted lines. The resonant width is measured along the isoline of
$\Gamma_{2}$, as shown in the _left panel_. For the retrograde 2:1 resonance,
the width is evaluated at $\theta_{2}=0$ and $\theta_{2}=\pi$ and, for the
retrograde 1:2 resonance, the width is evaluated at $\theta_{1}=\pi/2$ and
$\theta_{1}=3\pi/2$. It is observed that, for the retrograde 1:2 resonance,
the results evaluated at $\theta_{1}=\pi/2$ and $\theta_{1}=3\pi/2$ are
coincident (this is due to the symmetry between the islands of libration
centred at $\theta_{1}=\pi/2$ and $\theta_{1}=3\pi/2$). The nominal location
of resonance is marked by a short and black vertical line.

For a given motion integral $\Gamma_{2}$, we can analyze the associated
Poincaré sections to determine the location of libration centre, denoted by
$(a_{0},e_{0})$, the point on the left boundary, denoted by $(a_{L},e_{L})$,
and the point on the right boundary, denoted by $(a_{R},e_{R})$. When the
motion integral varies in a certain interval, the curve of the libration
centre as well as the points on the boundaries can be determined.

In Fig. 11, the numerical width of resonance and distribution of libration
centre are reported for the retrograde 2:1 resonance in the _left panel_ and
for the retrograde 1:2 resonance in the _right panel_. The shaded areas stand
for the libration regions in the $(a,e)$ space. For convenience, the nominal
resonance location is marked by a short and vertical black line and the level
curves of the motion integral are also presented. Similar to the analytical
results shown in Fig. 7, the resonant width should also be measured along the
isoline of the motion integral, i.e., $\Delta a=a_{R}-a_{L}$ and $\Delta
e=e_{R}-e_{L}$, as shown in the _left panel_. Observing Fig. 11, we can see
that (a) the location of libration centre is different from the nominal
location of resonance, (b) the left and right separatrices are not symmetric
with respect the nominal resonance location and (c) the resonant width is an
increasing function of the eccentricity.

In particular, for the retrograde 2:1 resonance, the resonant centres are
located on the left-hand side of the nominal resonance location and the
resonant widths evaluated at $\theta_{2}=0$ and $\theta_{2}=\pi$ have slight
difference. The difference can be understood by the asymmetry between the
island of libration centred at $\theta_{2}=0$ and the one at $\theta_{2}=\pi$
(please see Fig. 9 for the asymmetric structures). For the retrograde 1:2
resonance, the resonant centres are placed on the right-hand side of the
nominal location of resonance and it is found that the resonant widths
evaluated at $\theta_{1}=\pi/2$ and $\theta_{1}=3\pi/2$ are perfectly
coincident. This coincidence is due to the symmetry between the islands of
libration centred at $\theta_{1}=\pi/2$ and $\theta_{1}=3\pi/2$ (please see
Fig. 10 for the symmetric structures).

From now on, we denote the outcomes produced from the resonant model as
analytical results and the ones obtained by analyzing Poincaré sections as
numerical results.

##  6 Comparison between analytical and numerical results

In this section, let us compare the analytical results obtained from the
resonant model discussed in Sections 3 and 4 with the numerical results
produced by analyzing Poincaré sections described in Section 5 in order to
validate the analytical developments.

Figure 12: Comparisons between the analytical and numerical results of the
resonant centre and resonant width. Analytical results are obtained from the
resonant Hamiltonian model, while the numerical results are produced by
analyzing Poincaré sections. The analytical results are the same as the ones
shown in Fig. 7 and the numerical results are the same as the ones given in
Fig. 11. In each panel, the nominal location of resonance is indicated by a
short and black vertical line.

Figure 12 reports the analytical and numerical results in the $(a,e)$ space
for the distribution of resonant centre and resonant width for the retrograde
2:1 resonance in the first two panels and for the retrograde 1:2 resonance in
the last panel. The level curves of the motion integral $\Gamma_{2}$ are
plotted and the resonant width should be evaluated along the isoline of
$\Gamma_{2}$. Observing Fig. 12, we can see that the analytical results are in
quite good agreement with the numerical results in terms of the location of
resonant centre and resonant width, indicating that our analytical model
formulated in Section 3 is accurate to describe the dynamics of retrograde 2:1
and 1:2 resonances.

Figure 13: Phase portraits (_left panels_) as well as the Poincaré sections
(_right panels_) for the retrograde 3:1 and 1:3 resonances. For the retrograde
3:1 resonance, the resonant argument is
$\sigma=\frac{1}{3}(\lambda-3\lambda_{p}+2\varpi)=\frac{1}{3}\left(\theta_{1}-3\theta_{2}\right)$,
the motion integral is taken as $\Gamma_{2}=2.7$ (or $a_{\min}=0.456060$ in
normalized units) and the Poincaré section is defined by
$\theta_{1}=\lambda-\varpi=0^{\circ}$ (on the sections it holds
$\sigma=-\theta_{2}$). For the retrograde 1:3 resonance, the resonant argument
is
$\sigma=\frac{1}{3}(3\lambda-\lambda_{p}-2\varpi)=\frac{1}{3}\left(3\theta_{1}-\theta_{2}\right)$,
the motion integral is taken as $\Gamma_{2}=1.8$ (or $a_{\min}=1.824240$ in
normalized units) and the Poincaré section is defined by
$\theta_{2}=\lambda_{p}-\varpi=0^{\circ}$ (on the sections it holds
$\sigma=\theta_{1}$). The resonant width in terms of the variation of
semimajor axis, $\Delta a$, is marked in both the phase portraits and Poincaré
sections. It is interesting to observe that there is a perfect correspondence
between the phase portraits and the associated Poincaré sections.

Next, let us apply both the analytical and numerical approaches discussed
before to the retrograde 3:1 and 1:3 resonances. In Fig. 13, their phase
portraits and Poincaré sections with the same motion integral (or $a_{\min}$)
are reported. The _left panels_ are for the phase portraits and the right ones
are for the Poincaré sections. The curves passing through the saddle points
are shown in red lines, and the distance $\Delta a$ stands for the resonant
width in terms of the variation of semimajor axis.

For the retrograde 3:1 resonance, the resonant angle is
$\sigma=\frac{1}{3}(\theta_{1}-3\theta_{2})$, which becomes
$\sigma=-\theta_{2}$ on the sections defined by $\theta_{1}=\lambda-\varpi=0$.
Observing panels (a) and (b) of Fig. 13 (the motion integral is
$\Gamma_{2}=2.7$ or $a_{\min}=0.456060$), we can see that there are similar
structures arising in the phase portrait and in the Poincaré section: (i) the
resonant centres are located at $\sigma=\pm\pi/3,\pm\pi$ (corresponding to
$\varphi=\pi$) in the phase portrait or, equivalently, at
$\theta_{2}=\pm\pi/3,\pm\pi$ on the Poincaré section, (ii) the saddle points
are located at $\sigma=0,\pm\pi$ (corresponding to $\varphi=0$) in the phase
portrait or, equivalently, at $\theta_{2}=0,\pm\pi$ on the Poincaré section,
(iii) in both plots the red curves passing through saddle points play the role
of dynamical separatrices, dividing the entire phase space into regions of
libration and circulation, and (iv) the distance between the nearby
separatrices evaluated at the resonant centre stands for the resonant width,
denoted by $\Delta a$, measuring the size of libration region.

For the retrograde 1:3 resonance, the resonant angle is
$\sigma=\frac{1}{3}(3\theta_{1}-\theta_{2})$, which is equal to
$\sigma=\theta_{1}$ on the section defined by
$\theta_{2}=\lambda_{p}-\varpi=0$. Similarly, comparing panels (c) and (d) of
Fig. 13 (the motion integral is $\Gamma_{2}=1.8$ or $a_{\min}=1.824240$), we
can observe similar structures in both the phase portrait and the Poincaré
section: (i) the resonant centres are located at $\sigma=\pi/3,\pi,5\pi/3$
(corresponding to $\varphi=\pi$) in the phase portrait or, equivalently,
$\theta_{1}=\pi/3,\pi,5\pi/3$ on the section, and (ii) the saddle points are
placed at $\sigma=0,2\pi/3,4\pi/3$ (corresponding to $\varphi=0$) in the phase
portrait or, equivalently, $\theta_{1}=0,2\pi/3,4\pi/3$ on the section.

Figure 14: Comparisons between the analytical and numerical results in terms
of the resonant centre and resonant width. The _upper panels_ are for the
retrograde 3:1 resonance, and the _bottom panels_ are for the retrograde 1:3
resonance. For convenience, the level curves of $\Gamma_{2}$ are plotted and,
in particular, the isolines of $\Gamma_{2}=2.7$ and $\Gamma_{2}=1.8$ are
marked (the associated phase portraits and Poincaré sections are provided in
Fig. 13). For the retrograde 3:1 resonance, the numerical width of resonance
is evaluated at $\theta_{2}=\pm\pi/3$ (panel ‘a’) and $\theta_{2}=\pi$ (panel
‘b’) and, for the retrograde 1:3 resonance, the numerical width of resonance
is evaluated at $\theta_{1}=\pm\pi/3$ (panel ‘c’) and $\theta_{1}=\pi$ (panel
‘d’). The nominal location of resonance is indicated by a short and black
vertical line.

Figure 14 reports the analytical and numerical results in terms of the
distribution of resonant centre and resonant width in the _left panel_ for the
retrograde 3:1 resonance and in the _right panel_ for the retrograde 1:3
resonance. For convenience, the level curves of the motion integral
$\Gamma_{2}$ are plotted in both panels. It is noted that the numerical widths
of resonance are evaluated at $\theta_{2}=\pm\pi/3,\pi$ for the retrograde 3:1
resonance (see panels ‘a’ and ‘b’ of Fig. 14) and at $\theta_{1}=\pm\pi/3,\pi$
for the retrograde 1:3 resonance (see panels ‘c’ and ‘d’ of Fig. 14). For the
analytical width of resonance, the results have no change when the widths are
evaluated at different centres (due to the symmetry of phase portraits).
However, the numerical width of resonance will shift slightly if it is
evaluated at a different resonant centre, because the islands of libration in
the Poincaré section are not exactly symmetric (please refer to the Poincaré
sections shown in Fig. 13 for detailed structures).

According to Fig. 14, we can observe a good agreement between the analytical
and numerical results: (a) both the analytical and numerical widths of
resonance ($\Delta a$ or $\Delta e$) increase with the eccentricity, and (b)
the location of resonant centre is different from the nominal location. In
particular, it is observed from Fig. 14 that the resonant centres are located
on the left-hand (right-hand) side of the nominal location of resonance for
the retrograde 3:1 (1:3) resonance and their deviation is dependent on the
eccentricity. It is not difficult to understand that the slight deviation
between analytical and numerical results is caused by the ‘short-term’ effects
filtered in the analytical model.

In summary, the comparisons made in Fig. 12 for the retrograde 2:1 and 1:2
resonances and in Figs 13 and 14 for the retrograde 3:1 and 1:3 resonances
show that the analytical results could match well with the numerical results
in terms of the location of resonant centre and resonant width, indicating
that our analytical model is accurate and applicable in predicting the
dynamics of retrograde MMRs.

##  7 Numerical widths over the full range of eccentricity

In previous sections, we concentrate on the dynamics of retrograde MMRs with
eccentricities smaller than that of the planet-crossing orbit (i.e.,
$e<e_{c}$) considering the convergence of series expansion of disturbing
function as well as the availability of perturbation theory. In addition, the
analytical approach based on series expansion cannot deal with the retrograde
co-orbital resonance due to the divergence of Laplace coefficients with
$\alpha$ close to unity (Murray & Dermott, 1999). However, we know that the
non-perturbative technique based on computing Poincaré sections is not limited
by the planet-crossing condition and the co-orbital condition (Wang &
Malhotra, 2017; Malhotra et al., 2018; Lan & Malhotra, 2019), thus it is
possible to explore the dynamics of retrograde MMRs (including the interior,
co-orbital and exterior resonances) over the entire range of eccentricity
$e\in(0,1)$ in order to provide global pictures in the phase space. To this
end, we apply the numerical approach described in Section 5 to the retrograde
2:1, 1:1, 1:2, 1:3 and 1:4 resonances.

According to the traditional notations (Malhotra & Zhang, 2020), the resonant
centres with the usual critical argument at $\varphi=0$ belong to the
pericentric branch and the ones at $\varphi=\pi$ belong to the apocentric
branch. From equation (11), we have $\varphi=k\theta_{1}-k_{p}\theta_{2}$.
Thus, on the Poincaré sections defined by $\theta_{1}=0$ (inner and co-orbital
resonances), the resonance centres at $\theta_{2}=0$ (or $\theta_{2}=\pm
2\pi/k_{p}$) belong to the pericentric branch and the ones at
$\theta_{2}=\pm\pi/k_{p}$ belong to the apocentric branch. On the sections
defined by $\theta_{2}=0$ (outer resonances), the resonance centres at
$\theta_{1}=0$ (or $\theta_{1}=\pm 2\pi/k$) belong to the pericentric branch
and the ones at $\theta_{1}=\pm\pi/k$ belong to the apocentric branch.

Figure 15: Poincaré sections of the retrograde 2:1 resonance specified by
$\Gamma_{2}=2.22$ (_left panel_) and the numerical widths of resonance over
the full range of eccentricity (_right panel_). The numerical widths are
evaluated at $\theta_{2}=0$ and $\theta_{2}=\pm\pi/2$. In the _right panel_ ,
the level curves of the motion integral $\Gamma_{2}$ are presented, and the
isoline of $\Gamma_{2}=2.22$ and the planet-crossing line corresponding to
$a(1+e)=1$ are also plotted. The nominal location of resonance is indicated by
a short and vertical black line.

In the _left panel_ of Fig. 15, the Poincaré sections for the retrograde 2:1
resonance specified by $\Gamma_{2}=2.22$ are presented in the $(\theta_{2},e)$
and $(\theta_{2},a)$ spaces. According to the Poincaré sections, it is
observed that the islands of resonance are centred at $\theta_{2}=0$,
$\theta_{2}=\pm\pi/2$ and $\theta_{2}=\pm\pi$ and those libration islands
centred at $\theta_{2}=\pm\pi/2$ (or $\theta_{2}=\pm\pi$) are symmetric with
respect to the line of $\theta_{2}=0$. In addition, the eccentricity
($e$$\sim$0.6) is close to that of planet-crossing orbit, leading to strong
perturbation, and thus chaotic layers appear between neighboring islands.
According to the Poincaré sections, we can see that chaotic layers replace the
role of dynamical separatrices. By analyzing Poincaré sections, we can
numerically identify the resonant widths by evaluating the lower and upper
boundaries inside which resonance can always take place. In particular, the
resonant width is characterized by the variations of semimajor axis and
eccentricity ($\Delta a$ and $\Delta e$), as shown in the _left panel_ of Fig.
15.

The numerical widths of resonance over the full range of eccentricity are
reported in the _right panel_ of Fig. 15. It should be noted that the
numerical widths are evaluated at $\theta_{2}=0$ and $\theta_{2}=\pm\pi/2$
(the width evaluated at $\theta_{2}=\pi$ has a similar behavior to the one
evaluated at $\theta_{2}=0$, thus it is not considered here). For convenience,
the isoline of $\Gamma_{2}=2.22$ and the planet-crossing critical line
corresponding to $a(1+e)=1$ are provided. In the following discussions, we use
$e_{c}(a)=\frac{1}{a}-1$ to stand for the critical eccentricity of planet-
crossing orbit. In the region above the line of $e_{c}(a)$, close encounter
(or even collision) with planet may occur and, in the region below the
critical line, no collision can happen. In addition, the nominal resonance
location at $a=a_{0}$ is marked by a short and vertical black line (for the
retrograde 2:1 resonance, it holds $a_{0}=0.62976016$). From the _right panel_
of Fig. 15, it is observed that (a) the pericentric libration (with
$\theta_{2}=0$) occurs over the entire range of eccentricity while the
apocentric libration (with $\theta_{2}=\pm\pi/2$) appears only in the region
above the line of $e_{c}(a)$, (b) the resonant centres with
$\theta_{2}=\pm\pi/2$ (apocentric libration) are located on the right hand of
the nominal resonance location (i.e., $a=a_{0}$) and they diverge away from
$a=a_{0}$ when the eccentricity decreases close to $e_{c}(a)$, (c) the
resonant width evaluated at $\theta_{2}=\pm\pi/2$ shrinks as the eccentricity
approaches $e_{c}(a)$, (d) the resonant centres with $\theta_{2}=0$ are
located close to $a=a_{0}$ but on the left-hand side of it, and (e) the
resonant width evaluated at $\theta_{2}=0$ reaches the maximum in the vicinity
of $e_{c}(a)$.

Figure 16: Poincaré sections of the retrograde 1:1 resonance specified by
$\Gamma_{2}=1.48$ (_left panel_) and the numerical widths of resonance in the
full range of eccentricity (_right panel_). The location of asteroid 2015
BZ509 is marked by a black square.

For the retrograde co-orbital resonance, the angular variables
$\theta_{1}=\lambda-\varpi$ and $\theta_{2}=\lambda_{p}-\varpi$ have
comparable periods thus, in this case, both conditions illustrated in Fig. 8
can be utilized to produce Poincaré sections. In practice, we take the
definition of the inner case, i.e., the sections are defined by
$\theta_{1}=\lambda-\varpi=0$. As an example, the motion integral is taken at
$\Gamma_{2}=1.48$, and the associated Poincaré sections are reported in the
$(\theta_{2},e)$ and $(\theta_{2},a)$ spaces, as shown in the _left panel_ of
Fig. 16. From the Poincaré sections, it is observed that (a) the islands of
resonance are centred at $\theta_{2}=0$ (pericentric libration) and
$\theta_{2}=\pi$ (apocentric libration), (b) the island centred at
$\theta_{2}=\pi$ is much smaller than that at $\theta_{2}=0$, and (c) chaotic
motions can be found in the sections and they fill the space outside the
islands.

Similarly, the numerical widths of resonance ($\Delta a$ and $\Delta e$) can
be measured from the sections, and they are reported in the _right panel_ of
Fig. 16. For convenience, the isoline of $\Gamma_{2}=1.48$ and the nominal
resonance location (i.e., $a=a_{0}$) are marked. It is observed that (a) the
islands of resonance centred at $\theta_{2}=\pi$ (apocentric libration) occupy
two distinct regions (in particular, they disappear in the medium-eccentricity
regions), one region with low eccentricities is located on the right-hand side
of $a=a_{0}$ and the other region with high eccentricities occupies on the
left-hand side of $a=a_{0}$, (b) the resonant centres with $\theta_{2}=0$
(pericentric libration) are located on the right-hand side and they diverge
away from $a=a_{0}$ when the eccentricity approaches zero, and (c) the
resonant width evaluate at $\theta_{2}=0$ shrinks with eccentricity
approaching zero.

As discussed in the introduction, asteroid 2015 BZ509 is the first identified
asteroid located inside retrograde co-orbital resonance with Jupiter (Wiegert
et al., 2017). For convenience, the location of 2015 BZ509 is marked in the
$(a,e)$ space, as shown in Fig. 16. It is observed that the asteroid is
located near the centre of the pericentric libration region.

Figure 17: Poincaré sections of the retrograde 1:2 resonance specified by
$\Gamma_{2}=1.72$ (_left panel_) and the numerical widths of resonance in the
full range of eccentricity (_right panel_). The locations of asteroids 2008
SO218 and 1999 LE31 are marked by black squares.

Figure 17 presents the Poincaré sections and numerical widths for the
retrograde 1:2 resonance. For convenience, the nominal resonance location
($a=a_{0}$) is marked, and the isoline of $\Gamma_{2}=1.72$ and the planet-
crossing critical line determined by $a(1-e)=1$ are presented. In the
following discussions, we use $e_{c}(a)=1-\frac{1}{a}$ to stand for the
critical eccentricity of planet-crossing orbit. From the Poincaré sections, we
can observe that (a) the islands of resonance are located at $\theta_{1}=0$,
$\theta_{1}=\pm\pi/2$ and $\theta_{1}=\pm\pi$, (b) the islands centred at
$\theta_{1}=\pm\pi/2$ (or $\theta_{1}=\pm\pi$) are symmetric with respect to
the line of $\theta_{1}=0$ and (c) the separatrices are replaced by chaotic
layers. The numerical widths ($\Delta a$ and $\Delta e$) can be measured from
the Poincaré sections and, in particular, we consider the resonant widths
evaluated at $\theta_{1}=0$ and $\theta_{1}=\pm\pi/2$, as shown in the _right
panel_ of Fig. 17. It is noted that the width evaluated at $\theta_{1}=\pm\pi$
has a similar behavior to the one evaluated at $\theta_{1}=0$, thus it is not
considered here. According to the numerical widths of resonance, we can
observe that (a) the islands of resonance centred at $\theta_{1}=0$
(pericentric libration) appear in the region above the line of $e_{c}(a)$, (b)
as the eccentricity decreases close to $e_{c}(a)$, the resonant centres at
$\theta_{1}=0$ diverge away from $a=a_{0}$, (c) the resonant centres at
$\theta_{1}=\pm\pi/2$ (apocentric libration) stay close to $a=a_{0}$ but on
the right-hand side of it, and (d) the resonant width evaluated at
$\theta_{1}=\pm\pi/2$ reaches the maximum in the vicinity of $e_{c}(a)$.

According to Morais & Namouni (2013a) and Li et al. (2019), it is known that
asteroids 2008 SO218 and 1999 LE31 are currently inside the retrograde 1:2
resonance with Jupiter. For convenience, we mark their locations in the
$(a,e)$ space, as shown in Fig. 17. Evidently, they are located inside the
libration regions.

Figure 18: The numerical widths of resonance over the full range of
eccentricity for the retrograde 1:3 resonance (_left panel_) and for the
retrograde 1:4 resonance (_right panel_). The locations of retrograde
asteroids 2000 DG8 and 2014 AT28 (in retrograde 1:3 resonance) and 2016 LS and
2016 JK24 (in retrograde 1:4 resonance) are marked by black squares.

In Fig. 18, the numerical widths of the retrograde 1:3 and 1:4 resonances are
reported in the $(a,e)$ space. About these two resonances, only those islands
centred at $\theta_{1}=0$ and $\theta_{1}=\pm\pi/k$ are taken into
consideration. In both panels, the nominal resonant location ($a=a_{0}$) and
the planet-crossing line corresponding to $a(1-e)=1$ are provided (the
critical eccentricity is denoted by $e_{c}(a)=1-\frac{1}{a}$). Evidently, the
numerical results about the resonant widths for the 1:3 and 1:4 resonances are
similar to that of the retrograde 1:2 resonance shown in the _right panel_ of
Fig. 17, in terms of the following points: (a) the islands of resonance
centred at $\theta_{1}=0$ (pericentric libration) appear in the region above
the line of $e_{c}(a)$, (b) the resonant centres at $\theta_{1}=0$ diverge
away from $a=a_{0}$ when the eccentricity approaches $e_{c}(a)$, (c) the
resonant centres at $\theta_{1}=\pm\pi/k$ (apocentric libration) are located
close to $a=a_{0}$, and (d) the resonant width evaluated at
$\theta_{1}=\pm\pi/k$ reaches the maximum in the vicinity of $e_{c}(a)$.

According to Li et al. (2019), it is known that asteroids 2000 DG8 and 2014
AT28 (or 2016 LS and 2016 JK24) are potentially located inside the retrograde
1:3 (or 1:4) resonance with Jupiter. For convenience, their locations are
presented in the $(a,e)$ space, as shown in Fig. 18. It is observed that
asteroids 2000 DG8 and 2016 JK24 sit near the centre of the libration regions,
so they will be protected by the resonant configurations for a long enough
time. However, asteroid 2016 LS and 2014 AT28 lie close to the boundaries of
libration zones, resulting in weak protection effects in the long-term
evolution, so the asteroids may exit from the current resonant states. Li et
al. (2019) showed that 2016 LS would be possibly captured in the retrograde
3:5 resonance with Saturn (with possibility of 13.6%).

In recent years, Huang et al. (2018) discussed the resonant widths by means of
semianalytical approach regarding the retrograde 1:1 resonance, and Li et al.
(2020) applied the semianalytical method to the dynamics of retrograde 1:2,
1:3, 1:4 and 1:5 resonances. By analyzing the phase portraits, Huang et al.
(2018) and Li et al. (2020) measured resonant widths by evaluating the
distance of adjacent separatrices that bound libration islands when $e<e_{c}$
and of the largest boundaries that do not cross the collision curve for
$e>e_{c}$ (Morbidelli, 2002). In comparison to our numerical results produced
in this section, some discussions are made as follows.

  * •

For the retrograde 1:1 resonance, Huang et al. (2018) noticed that the width
evaluated at the usual critical argument $\varphi=0$ (the pericentric branch)
keeps growing with increasing eccentricity (see fig. 2 in their paper).
Qualitatively, this behavior of width is in quite good agreement with our
numerical results (see the right panel of Fig. 16). For the apocentric branch,
Huang et al. (2018) observed that no apocentric libration with $\varphi=\pi$
exists when the motion integral changes from $-1.94$ to $-1.65$ (corresponding
to $\Gamma_{2}$ used in this work ranging from 1.65 to 1.94). The authors
pointed out that the resonant centres in the apocentric branch survives when
the eccentricity is small or extremely high. This behavior is reproduced in
our numerical results. Please refer to the _right panel_ of Fig. 16 for
details, which shows that the apocentric libration exists in two distinct
regions in the $(a,e)$ space: one with small eccentricities and the other one
with high eccentricities (in the medium-eccentricity region the apocentric
libration vanishes).

  * •

Concerning the retrograde 1:2, 1:3, 1:4 and 1:5 resonances, Li et al. (2020)
reported their resonant widths in the $(a,e)$ space and they noticed that all
the retrograde 1:$n$ resonances have similar structures (see figs 7 and 10 for
details). The authors assumed that the resonance occurs at the nominal
resonance location ($a=a_{0}$), and their results show that the boundaries are
symmetric with respect to $a=a_{0}$. From the numerical viewpoint, we confirm
that all the retrograde 1:2, 1:3 and 1:4 resonances considered in this work
hold similar structures about resonant widths. However, some discrepancies are
observed. Firstly, our numerical results (see Figs 17 and 18) show that the
resonant centres are not located at the nominal resonance locations and, in
particular, the libration centres in the pericentric branch ($\theta_{1}=0$)
diverge away from the nominal resonance location when the eccentricity is
close to $e_{c}$. Secondly, our results indicate that the boundaries in the
$(a,e)$ are not symmetric with respect to the line of $a=a_{0}$. At last, our
results show that the apocentric libration (with $\theta_{1}=\pm\pi/k$)
vanishes when the eccentricity is higher than a threshold value.

##  8 Summary and discussion

In this work, we performed both analytical and numerical studies about the
dynamics of retrograde MMRs in the planar circular restricted three-body
problem and, in particular, we made direct comparisons between the analytical
and numerical results, including the dynamical structures in the phase space,
location of resonant centre and resonant width.

Regarding the analytical study, we first presented the explicit expansion of
disturbing function for test particles moving on retrograde co-planar orbits
and then formulated the Hamiltonian model by (a) introducing a new critical
argument as $\sigma=\frac{1}{k_{\max}}\left[k\lambda-
k_{p}\lambda_{p}+(k_{p}-k)\varpi\right]$ with
$k_{\max}={\max}\left\\{k_{p},k\right\\}$ and (b) performing a series of
canonical transformations. The resulting resonant Hamiltonian determines a
single-degree of freedom dynamical model, and the global dynamics of
retrograde MMRs can be explored by analyzing phase portraits, where the
resonant centre and resonant width can be analytically identified
(corresponding to analytical results).

Our numerical investigation is based on the non-perturbative technique by
analyzing Poincaré sections. We produced Poincaré surfaces of section by
recording those states of test particles when the ‘short-period’ angular
variable is equal to zero under the condition that the motion integral
$\Gamma_{2}=\sqrt{\mu a}(\frac{k_{p}}{k}+\sqrt{1-e^{2}})$ is provided. In
particular, the section is defined by $\theta_{1}=\lambda-\varpi=0$ for inner
resonances and by $\theta_{2}=\lambda_{p}-\varpi=0$ for outer resonances. It
is found that, on the Poincaré sections, the angular separation ($\theta_{2}$
for the inner resonances or $\theta_{1}$ for the outer resonances) is equal to
the synodic angle between the test particle and the perturber in magnitude. By
analyzing the structures arising in the sections, it is possible to
numerically determine the location of resonant centre and resonant width
(corresponding to numerical results).

According to the definition of sections, there is a relationship between the
resonant angle $\sigma$ and the angular separation on the sections:
$\sigma=-\theta_{2}$ for the inner resonances and $\sigma=\theta_{1}$ for the
outer resonances. This relationship makes it possible to compare the phase
portraits and Poincaré sections for a certain retrograde MMR with a given
motion integral $\Gamma_{2}$. Naturally, we can make a comparison between the
analytical and numerical widths of resonance. The main results of the
comparative study are summarized as follows.

  * •

As for a certain retrograde MMR, there is a perfect correspondence between the
phase portraits and Poincaré sections with the same motion integral
$\Gamma_{2}$. This correspondence is helpful to understand the structures
arising in the Poincaré sections and also to validate our analytical
developments.

  * •

The structures arising in the phase portrait are in good agreement with the
ones appearing in the Poincaré section, including the islands of libration,
dynamical separatrices, resonant centres and saddle points.

  * •

For the retrograde MMRs considered in this study, the analytical results
(including the resonant width and the location of resonant centre) obtained
from the resonant model agree well with the associated numerical results
produced by analyzing Poincaré sections. It shows that our analytical model
formulated in this work is applicable to describe the dynamics of retrograde
MMRs in the non-crossing regions.

  * •

For the retrograde 2:1 resonance, it is found that the zero-eccentricity point
(i.e., the origin in the polar coordinate plane) is no longer a stationary
point in both the phase portraits and Poincaré sections. This is different
from the prograde 2:1 resonance (please refer to Malhotra & Zhang (2020) and
Lei & Li (2020) for discussions on the prograde 2:1 resonance).

  * •

For the outer resonances (i.e., the retrograde 1:2 and 1:3 resonances), both
the analytical and numerical results indicate that there are no asymmetric
libration centres. This is in agreement with the conclusion given by Li et al.
(2020).

  * •

Both the analytical and numerical results show that the width of resonance
increases with the eccentricity, the location of resonant centre is dependent
on the eccentricity and it deviates from its nominal resonance location (i.e.,
the Law of Structure exists for retrograde MMRs).

The numerical approach based on Poincaré sections is not limited to the non-
crossing and non-coorbital conditions. Thus, in the last part of this work, we
applied the numerical technique to identifying dynamical structures of
retrograde MMRs (including the interior, co-orbital and exterior resonances)
over the full range of eccentricity. For the retrograde co-orbital resonances,
it is found that (i) the apocentric libration occurs in two distinctive
regions in the phase space (one with low eccentricities and the other one with
high eccentricities) and (ii) the libration centres in the pericentric branch
diverge away from the nominal resonance location with eccentricity deceasing
close to zero. About the retrograde interior and exterior MMRs, it is found
that (a) in the non-crossing regions there is one branch of libration centres
while in the planet-crossing regions an additional branch of libration centres
appears, (b) chaotic layers replace the role of dynamical separatrices in the
planet-crossing regions, making that the boundaries are no longer smooth
curves, and (c) the libration centres located in the planet-crossing regions
diverge away from the nominal resonance location when the eccentricity
decreases close to the critical value $e_{c}(a)$. In particular, the
retrograde asteroids 2015 BZ509 (in retrograde coorbital resonance), 2008
SO218 and 1999 LE31 (in retrograde 1:2 resonance), 2000 DG8 and 2014 AT28 (in
retrograde 1:3 resonance), and 2016 LS and 2016 JK24 (in retrograde 1:4
resonance) are located inside the libration regions (please refer to Figs
16–18 for details), indicating that these seven asteroids are inside the
retrograde resonances with Jupiter.

Figure 19: The numerical widths of retrograde 1:1, 1:2, 1:3 and 1:4 resonances
with Jupiter, together with all the known retrograde asteroids with semimajor
axes $a\in(3.12\;\mathrm{au},14.04\;\mathrm{au})$. For convenience, seven
potential examples of asteroids inside retrograde MMRs, which are confirmed in
previous works, are indicated.

Figure 19 reports the numerical widths of resonance together with the
retrograde asteroids222https://minorplanetcenter.net//iau/MPCORB.html,
retrieved 3 February 2020 with semimajor axes
$a\in(3.12\;\mathrm{au},14.04\;\mathrm{au})$. Evidently, besides the seven
known examples of asteroids potentially inside retrograde MMRs with Jupiter,
there are additional asteroids which are located inside the libration zones
(for example, it is observed that asteroid 2018 TL6 is located at the centre
of the libration zone of the retrograde 1:2 resonance). The asteroids inside
resonance zones could be considered as new potential candidates to be trapped
inside retrograde MMRs. However, we need to notice two points: (a) the PCRTBP
adopted in this work is a greatly simplified model for the Solar system, so
that the results obtained in such a simplified model could only provide
preliminary boundaries for retrograde MMRs, and (b) the inclinations of
practical retrograde asteroids are not equal to $180^{\circ}$. Thus, whether
these asteroids are really trapped inside resonances or not, it requires
numerical integrations to confirm them in the Solar system model. We leave the
confirmation to a future study.

## Acknowledgments

This work is supported by the National Natural Science Foundation of China
(Nos 12073011, 11973027, 41774038). This research has made use of data and/or
services provided by the International Astronomical Union’s Minor Planet
Centre.

## Data availability

The analysis and codes are available upon request.

## References

  * Beaugé (1994) Beaugé C., 1994, Celest. Mech. Dyn. Astron., 60, 225
  * Chen et al. (2016) Chen Y.-T., Lin H. W., Holman M. J., Payne M. J., Fraser W. C., Lacerda P., Ip W.-H., Chen W.-P., Kudritzki R.-P., Jedicke R., et al., 2016, ApJ, 827, L24
  * Connors & Wiegert (2018) Connors M., Wiegert P., 2018, Planet. Space Sci., 151, 71
  * Ellis & Murray (2000) Ellis K. M., Murray C. D., 2000, Icarus, 147, 129
  * Fehlberg (1968) Fehlberg E., 1968, Classical fifth-, sixth-, seventh-, and eighth-order Runge-Kutta formulas with step-size control. NASA
  * Gallardo (2006a) Gallardo T., 2006a, Icarus, 184, 29
  * Gallardo (2006b) Gallardo T., 2006b, Icarus, 181, 205
  * Gallardo (2019a) Gallardo T., 2019a, MNRAS, 487, 1709
  * Gallardo (2019b) Gallardo T., 2019b, Icarus, 317, 121
  * Gallardo (2020) Gallardo T., 2020, Celest. Mech. Dyn. Astron., 132, 9
  * Gladman et al. (2009) Gladman B., Kavelaars J., Petit J.-M., Ashby M., Parker J., Coffey J., Jones R., Rousselot P., Mousis O., 2009, ApJ, 697, L91
  * Huang et al. (2018) Huang Y., Li M., Li J., Gong S., 2018, AJ, 155, 262
  * Hughes (1981) Hughes S., 1981, Celest. Mech., 25, 101
  * Kotoulas & Voyatzis (2020a) Kotoulas T., Voyatzis G., 2020a, Planet. Space Sci., 182, 104846
  * Kotoulas & Voyatzis (2020b) Kotoulas T., Voyatzis G., 2020b, Celest. Mech. Dyn. Astron., 132, 1
  * Lan & Malhotra (2019) Lan L., Malhotra R., 2019, Celest. Mech. Dyn. Astron., 131, 39
  * Lei & Li (2020) Lei H., Li J., 2020, MNRAS, 499, 4887
  * Li et al. (2020) Li J., Lawler S., Zhou L.-Y., Sun Y.-S., 2020, MNRAS, 492, 3566
  * Li et al. (2014a) Li J., Zhou L.-Y., Sun Y.-S., 2014a, MNRAS, 437, 215
  * Li et al. (2014b) Li J., Zhou L.-Y., Sun Y.-S., 2014b, MNRAS, 443, 1346
  * Li et al. (2018) Li M., Huang Y., Gong S., 2018, A&A, 617, A114
  * Li et al. (2019) Li M., Huang Y., Gong S., 2019, A&A, 630, A60
  * Li et al. (2020) Li M., Huang Y., Gong S., 2020, ApSS, 365, 1
  * Malhotra (1996) Malhotra R., 1996, AJ, 111, 504
  * Malhotra et al. (2018) Malhotra R., Lan L., Volk K., Wang X., 2018, AJ, 156, 55
  * Malhotra & Zhang (2020) Malhotra R., Zhang N., 2020, MNRAS, 496, 3152
  * Morais & Giuppone (2012) Morais M., Giuppone C., 2012, MNRAS, 424, 52
  * Morais & Namouni (2013a) Morais M., Namouni F., 2013a, MNRAS, 436, L30
  * Morais & Namouni (2013b) Morais M., Namouni F., 2013b, Celest. Mech. Dyn. Astron., 117, 405
  * Morais & Namouni (2017) Morais M., Namouni F., 2017, MNRAS, 472, L1
  * Morais & Namouni (2019) Morais M., Namouni F., 2019, MNRAS, 490, 3799
  * Morbidelli (2002) Morbidelli A., 2002, Modern celestial mechanics: aspects of solar system dynamics. Taylor & Francis, London and New York
  * Murray & Dermott (1999) Murray C. D., Dermott S. F., 1999, Solar system dynamics. Cambridge university press
  * Namouni & Morais (2018) Namouni F., Morais H., 2018, Comput. Appl. Math., 37, 65
  * Shevchenko (2016) Shevchenko I. I., 2016, The Lidov-Kozai effect-applications in exoplanet research and dynamical astronomy. Vol. 441, Springer
  * Wang & Malhotra (2017) Wang X., Malhotra R., 2017, AJ, 154, 20
  * Wiegert et al. (2017) Wiegert P., Connors M., Veillet C., 2017, Nature, 543, 687
  * Winter & Murray (1997a) Winter O., Murray C., 1997a, A&A, pp 290–304
  * Winter & Murray (1997b) Winter O., Murray C., 1997b, A&A, pp 399–408

##  Appendix A Resonant disturbing function

For the $k_{p}$:$k$ resonance, the critical argument is given by
$\sigma=\frac{1}{k_{\max}}\left[k\lambda-
k_{p}\lambda_{p}+(k_{p}-k)\varpi\right]$. The resonant disturbing function can
be written as

$\displaystyle{{\cal R}^{*}}$ $\displaystyle={\cal
G}{m_{p}}\sum\limits_{n=0}^{N}{\sum\limits_{\scriptstyle
j=-\infty\hfill\atop\scriptstyle\bmod(j,{k_{p}})=0\hfill}^{\infty}{\sum\limits_{m=0}^{n}{{{\left({-1}\right)}^{n-m}}{A_{n,j}}\left(\alpha\right){n\choose
m}}}}$ (21) $\displaystyle\times
X_{-jk/{k_{p}}}^{m,j}\left(e\right)\cos\left({\frac{j}{{{k_{p}}}}{k_{\max}}\sigma}\right)\;$
$\displaystyle-{\cal
G}{m_{p}}aX_{-k}^{1,1}\left(e\right){\delta_{{k_{p}},1}}\cos\left({{k_{\max}}\sigma}\right)$

where ${\delta_{k_{p},1}}=1$ when $k_{p}=1$ and ${\delta_{k_{p},1}}=0$ when
$k_{p}\neq 1$. For convenience, the resonant disturbing function is denoted by

${{\cal R}^{*}}=\sum\limits_{q=0}^{Q}{{{\cal
C}_{q}}\cos\left({q{k_{\max}}\sigma}\right)},$ (22)

where the coefficients ${\cal C}_{q}$ are related to the semimajor axis and
eccentricity. Based on equation (21), it is not difficult to obtain

$\displaystyle{{\cal C}_{q}}\left({a,e}\right)$ $\displaystyle={\cal
G}{m_{p}}\sum\limits_{n=0}^{N}{\sum\limits_{m=0}^{n}{{{\left({-1}\right)}^{n-m}}{A_{n,q{k_{p}}}}\left(\alpha\right){n\choose
m}X_{-qk}^{m,q{k_{p}}}\left(e\right)}}$ $\displaystyle-{\cal
G}{m_{p}}aX_{-k}^{1,1}\left(e\right){\delta_{k_{p},1}}{\delta_{q,1}}.$