# Dynamics around non-spherical symmetric bodies - I.
The case of a spherical body with mass anomaly

G. Madeira, S. M. Giuliatti Winter, T. Ribeiro, O. C. Winter

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

► To cite this version:

G. Madeira, S. M. Giuliatti Winter, T. Ribeiro, O. C. Winter. Dynamics around non-spherical symmetric bodies - I. The case of a spherical body with mass anomaly. Monthly Notices of the Royal Astronomical Society, 2022, 510, pp.1450-1469. 10.1093/mnras/stab3552. insu-03644352

HAL Id: insu-03644352

https://insu.hal.science/insu-03644352

Submitted on 20 Apr 2022

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

# Dynamics around non-spherical symmetric bodies – I. The case of a spherical body with mass anomaly

G. Madeira<sup>*</sup>,<sup>1,2★</sup> S. M. Giuliatti Winter,<sup>1★</sup> T. Ribeiro<sup>1</sup> and O. C. Winter<sup>1</sup>

Accepted 2021 December 2. Received 2021 November 10; in original form 2021 August 6

## ABSTRACT

The space missions designed to visit small bodies of the Solar system boosted the study of the dynamics around non-spherical bodies. In this vein, we study the dynamics around a class of objects classified by us as non-spherical symmetric bodies, including contact binaries, triaxial ellipsoids, and spherical bodies with a mass anomaly, among others. In this work, we address the results for a body with a mass anomaly. We apply the pendulum model to obtain the width of the spin-orbit resonances raised by non-asymmetric gravitational terms of the central object. The Poincaré surface of section technique is adopted to confront our analytical results and to study the system’s dynamics by varying the parameters of the central object. We verify the existence of two distinct regions around an object with a mass anomaly: a chaotic inner region that extends beyond the corotation radius and a stable outer region. In the latter, we identify structures remarkably similar to those of the classical restrict and planar three-body problem in the Poincaré surface of sections, including asymmetric periodic orbits associated with 1:1+p resonances. We apply our results to a Chariklo with a mass anomaly, obtaining that Chariklo rings are probably related to first kind periodic orbits and not with 1:3 spin-orbit resonance, as proposed in the literature. We believe that our work presents the first tools for studying mass anomaly systems.

**Key words:** celestial mechanics – Kuiper belt objects: individual: (10199) Chariklo – minor planets, asteroids: general – planets and satellites: dynamical evolution and stability – planets and satellites: rings.

## 1 INTRODUCTION

In the last three decades, the acquisition of data on the shape of small heliocentric bodies, by ground- and space-based observations (Hudson & Ostro 1995; Hanuš, Marchis & Ďurech 2013; Hanuš et al. 2017) and by space-mission explorations – such as OSIRIS-REx and Hayabusa spacecraft (Yoshikawa et al. 2015; Lauretta et al. 2017) – fostered the study of the dynamics around these bodies. This class of objects, which involves asteroids, trans-Neptunian objects, Centaurs, and comets, characteristically has diameters of less than 1000 km (Jorda et al. 2016). Due to their small sizes, these bodies do not have enough mass to reach hydrostatic equilibrium, showing irregular and asymmetric shapes.

The development of space missions was a strong motivation for the search of equilibrium regions around irregular bodies, as accomplished, e.g. by Scheeres, Williams & Miller (2000), which obtained stable orbits around the asteroid 433 Eros for the spacecraft NEAR-Shoemaker (Prockter et al. 2002). Some other works with such purpose are Yu & Baoyin (2012), Shang, Wu & Cui (2015), Wang, Li & Gong (2016), Winter et al. (2020), and Moura et al. (2020). The discovery of satellites and rings around this class of objects was also justifications for the interest in the stability of

irregular bodies systems (Chapman et al. 1995; Merline et al. 2002; Braga-Ribas et al. 2014; Ortiz et al. 2017).

When investigating the motion around irregular bodies, it is essential to consider the gravitational field generated by their odd shape. One method used for this is to approximate the irregular shape to a symmetric one – such as a MacLaurin spheroid or a triaxial ellipsoid – which allows studying the system theoretically or through low-cost simulations. Another course of action is to decompose the irregular body into a set of regular polyhedra (Polyhedron Shape Model; Werner 1994) or mass points (Mascon Model, Geissler et al. 1996). Despite the high level of accuracy, this methodology has a higher computational cost.

In the current and subsequent works, we study the dynamics around a class of objects classified by us as non-spherical symmetric bodies (NSSBs): contact binaries, triaxial ellipsoids with uniform density, and spherical bodies with a mass anomaly. The motion around NSSBs has already been studied in some articles, such as Lages, Shepelyansky & Shevchenko (2017) that analysed the stability around contact binaries through a generalized Kepler map technique (Meiss 1992; Shevchenko 2011), obtaining chaotic gravitational zones around the central body, similar to those found for symmetrical elongated bodies (Mysen, Olsen & Aksnes 2006; Mysen & Aksnes 2007). Their results are applicable to the asteroids 243 Ida and 25143 Itokawa (Lages et al. 2017).

Lages, Shevchenko & Rollin (2018) also use the Kepler map technique to study the chaotic region around cometary nuclei of

<footnote n="1"><sup>1</sup>Grupo de Dinâmica Orbital and Planetologia, São Paulo State University – UNESP, Avenida Ariberto Pereira da Cunha, 333, Guaratinguetá, SP 12516-410, Brazil</footnote>

<footnote n="2"><sup>2</sup>Institut de Physique du Globe de Paris, CNRS, Université de Paris, F-75005 Paris, France</footnote>

<footnote>* E-mail: gustavo.o.madeira@unesp.br (GM); giuliatti.winter@unesp.br (SMGW)</footnote>
dumb-bell shape, obtaining that such region is responsible for engulfing most of the Hill sphere of Comet 1P/Halley. Amarante & Winter (2020) studied the dynamics around 486958 Arroloth, an object similar to a contact binary, using a Polyhedron Shape Model and found an unstable zone in the equatorial region of the asteroid. Rollin, Shevchenko & Lages (2021) obtain that the particles in the equatorial plane of 486958 Arroloth are lost due to the chaotic diffusion of the orbits, which results in collisions or particle ejection. Interestingly, Rollin et al. (2021) also obtain theoretical dumb-bell-shaped objects with certain combinations of mass and spin period that host not a complete chaotic zone, but a chaotic ring.

The dynamics around triaxial ellipsoids were previously studied by Scheeres (1994) and Vantieghem (2014), and in particular by Winter, Borderes-Motta & Ribeiro (2019), who analysed the motion around 136108 Haumea, an ellipsoidal-shape object. This dwarf planet is particularly interesting due to its complex system that includes a pair of satellites, Hi'iaka and Namaka, and a ring (Ragozzine & Brown 2009; Ortiz et al. 2017). The non-asymmetric terms of the gravitational field of the NSSBs create strong resonances between the orbital period of the ring particles and the spin of the central body. Ortiz et al. (2017) propose that the Haumea ring would be associated with the 1:3 resonance. However, Winter et al. (2019) using the Poincaré surface of sections showed that this resonance is doubled, generating a large chaotic region in the resonance separatrix. Consequently, the ring is not associated with the 1:3 resonance but probably with first kind periodic orbits.

10199 Chariklo is another irregular body with a complex system involving a pair of narrow rings and possibly small satellites (Braga-Ribas et al. 2014; Bérard et al. 2017). The shape of Chariklo is still not well known. Observational data suggest triaxial and Jacobi ellipsoid shapes for the object (Leiva et al. 2017). Sicardy (2020) discusses the possibility of Chariklo to be a sphere with topographical features of a few kilometres, i.e. an object with a mass anomaly. Assuming a spherical Chariklo with a mass anomaly, Sicardy et al. (2019, 2020) show that particles inside the corotation radius migrate on to the body, and the outer material is pushed beyond the 1:2 resonance.

Here, we apply some well-known techniques to study the dynamics around a spherical body with a mass anomaly. Relations for the width and location of the spin-orbit resonances, also known as sectoral resonances, are presented. The dependence of resonances on the central body parameters is analysed. The Poincaré surface of section technique is also applied to the system for analysing the stability of the particles. We advance the reader the existence of a chaotic region near the object with a mass anomaly. This region extension is measured, and an adjusted equation is obtained as a function of the system parameters.

In Section 2, we present the disturbing function of our case of interest. In Section 3, we follow the prescription of the pendulum model developed by Winter & Murray (1997a) and Murray & Dermott (1999) for the restricted planar three-body problem (RP3BP) to obtain an analytical recipe for the location and width of the spin-orbit resonances. Section 4 presents the the Poincaré surface of section technique (Hénon 1965a, b, 1966a, b, 1969; Jefferys 1971). In Section 5, we identify, through numerical simulations, stable regions and give an overview of the system. In Section 6, we use the Poincaré surface of section technique to confront our analytical model and study the spin-orbit resonances in details. We apply our results to Chariklo in Section 7, exploring the dynamics around the object, in particular in the region of the rings. We address our final comments in Section 8.

<figure><img src="image_4.png" /><figcaption>**Figure 1.** Schematic diagram of the trajectory of a particle around a spherical object with a mass anomaly at its equator. The trajectory is fixed in the rotating frame with the central body's angular velocity $\omega$. $x_0$ is the initial position of the particle, and the red arrow indicates the initial velocity.</figcaption></figure>

## 2 DYNAMICAL SYSTEM

In this work, we analyse the dynamics of particles orbiting a hypothetical spherical object of mass $M$ and radius $R$, with a mass anomaly $m_a$ at its equator (Fig. 1). We assume the object with a uniform mass distribution, where the masses $M$ and $m_a$ have the same bulk density ($\rho = 1$ g cm$^{-3}$). The object is also assumed to rotate with constant angular velocity $\omega$ ($\omega = 2\pi/T$, where $T$ is the rotation period) without wobbling motion. For simplicity, we will express our physical quantities in the following units: $GM = 1$, while $R = 1$ is the distance between the system centre and the mass anomaly. We also define as a unit the Keplerian frequency of the mass anomaly, scaled by the density $\rho$ of the object:

$$ \omega_k = \sqrt{\frac{GM}{R^3}} = \sqrt{\frac{4\pi G\rho}{3}} = 1. \quad (1) $$

Two dimensionless parameters will define our dynamic system: the normalized mass anomaly $\mu = m_a/M$ and the rotating rate $\lambda = \omega/\omega_k$.

Equations of motion in a frame $Oxy$ rotating with the same period as the central body's spin are given by (Scheeres et al. 1996)

$$ \ddot{x} - 2\lambda\dot{y} = \lambda^2 x + U_x \quad (2) $$

and

$$ \ddot{y} + 2\lambda\dot{x} = \lambda^2 y + U_y, \quad (3) $$

where $U_x$ and $U_y$ stand for the partial derivatives of the gravitational potential.

The potential acting on a particle with position vector $\vec{r} = x\hat{x} + y\hat{y}$ ($r = |\vec{r}|$) in the rotating frame is obtained by adding the gravitational potential of the spherical portion of the object – at the centre of the system – with the gravitational potential of the mass anomaly, located at $\vec{R} = \hat{x}$ (Sicardy et al. 2019):

$$ U(r) = -\frac{1}{r} - \frac{\mu}{|\vec{r} - \hat{x}|} + \lambda^2 \mu (\vec{r} \cdot \hat{x}). \quad (4) $$

Note that the potential given in equation (4) differs from that acting on a particle in the RP3BP (Murray & Dermott 1999) by the rotating parameter $\lambda^2$ in the indirect term. While the secondary mass in RP3BP surrounds the central body with Keplerian velocity $\omega_k$, here the mass anomaly rotates with angular velocity $\lambda\omega_k$. We introduced the rotating parameter to correct this difference.

Similar to the dynamics of a particle in the RP3BP with a disturbing internal body, we obtain the expansion of the potential $U$ for the lowest order terms in eccentricity ($e$) as

$$ U = -\frac{1}{r} - \sum_{j=0}^{\infty} \sum_{m=-\infty}^{\infty} \mu e^j \left[ \alpha F_j b_{1/2}^{(m-j)}(\alpha) + \frac{\lambda^2}{\alpha} f_j \delta_{|m|,1} \right] \cos \phi, \quad (5) $$
<table><caption>**Table 1.** The linear operators $f_j$ and $F_j$ for $j \le 5$. The derivative operator $D$ is given by $D = d/d\alpha$.</caption><thead><tr><th>j</th><th>f<sub>j</sub></th><th>F<sub>j</sub></th></tr></thead><tbody><tr><td>1</td><td>-<sup>1</sup>&frasl;<sub>2</sub> [(−1 + 2m) + αD]</td><td></td></tr><tr><td>2</td><td>-<sup>3</sup>&frasl;<sub>8</sub> [(2 − 7m + 4m<sup>2</sup>) + (−2 + 4m)αD + α<sup>2</sup>D]</td><td></td></tr><tr><td>3</td><td>-<sup>1</sup>&frasl;<sub>3</sub> <sup>1</sup>&frasl;<sub>48</sub> [(−6 + 29m − 30m<sup>2</sup> + 8m<sup>3</sup>) + (6 − 21m + 12m<sup>2</sup>)αD + (−3 + 6m)α<sup>2</sup>D + α<sup>3</sup>D<sup>3</sup>]</td><td></td></tr><tr><td>4</td><td>-<sup>125</sup>&frasl;<sub>384</sub> <sup>1</sup>&frasl;<sub>384</sub> [(24 − 146m + 211m<sup>2</sup> − 104m<sup>3</sup> + 16m<sup>4</sup>) + (−24 + 116m − 120m<sup>2</sup> + 32m<sup>3</sup>)αD + (12 − 42m + 24m<sup>2</sup>)α<sup>2</sup>D + (−4 + 8m)α<sup>3</sup>D<sup>3</sup> + α<sup>4</sup>D<sup>4</sup>]</td><td></td></tr><tr><td>5</td><td>-<sup>27</sup>&frasl;<sub>80</sub> <sup>1</sup>&frasl;<sub>3840</sub> [(−120 + 874m − 1595m<sup>2</sup> + 1110m<sup>3</sup> − 320m<sup>4</sup> + 32m<sup>5</sup>) + (120 − 730m + 1055m<sup>2</sup> − 520m<sup>3</sup> + 80m<sup>4</sup>)αD + (−60 − 290m − 300m<sup>2</sup> + 80m<sup>3</sup>)α<sup>2</sup>D + (20 − 70m + 40m<sup>2</sup>)α<sup>3</sup>D<sup>3</sup> + (5 + 10m)α<sup>4</sup>D<sup>4</sup> + α<sup>5</sup>D<sup>5</sup>]</td><td></td></tr></tbody></table>

where $\alpha = 1/a < 1$ ($a$ = semimajor axis of the particle), $b_{1/2}^{(m)}$ is the Laplace coefficient, $f_j$ and $F_j$ are linear operators (Table 1; Murray & Dermott 1999; Ellis & Murray 2000), $\delta_{|m|,1}$ is the Kronecker delta, and $\phi$ is a characteristic angle of the system relating the rotation of the central object with the longitudes of the particle. The characteristic angle associated with the sectoral resonances is presented in Section 3.

In conservative systems, such as those analysed in this work, the Jacobi constant $C_j$ is a conserved quantity used to obtain the Poincaré surface of sections. It is expressed here in the units $R^2\omega_k^2$ and is given by (Scheeres et al. 1996)

$$C_j = \lambda^2(x^2 + y^2) + 2U(x, y) - \dot{x}^2 - \dot{y}^2. \quad (6)$$

### 3 SECTORAL RESONANCES

At the planar limit, a pair of fundamental frequencies describe the motion of a particle: the synodic and radial epicyclic frequencies. The first, $n - \omega$ ($n$ = angular frequency of the particle), corresponds to the frequency of the particle's return to a fixed position on the rotating frame. The second, $\kappa = n - \dot{\sigma}$, is the frequency of the particle's return to its pericentre, being $\dot{\sigma}$ the derivative of the particle's longitude of pericentre. If these frequencies are commensurable, the particle is in a sectoral resonance – spin–orbit resonance – with the central body. Once in resonance, the orbital evolution of the particle will be modelled by the energy balance provided by the resonant configuration. Sectoral resonances with real non-spherical bodies were studied in Borderes-Motta & Winter (2018) and Winter et al. (2019) for the asteroid 4179 Toutatis and the dwarf planet Haumea, respectively.

A particle at the centre of an $m:(m-j)$ resonance satisfies the resonance condition (Sicardy et al. 2019)

$$m\omega - (m - j)n - j\dot{\sigma} = 0, \quad (7)$$

where $m$ and $j$ are integers responsible for giving the commensurability of the frequencies. For $j = 0$, the particle is in corotation resonance, while for $j = m$, we have the apsidal resonances. Both cases are out of the scope of this work (for details, see Sicardy et al. 2019) and here we will focus on resonances with $j \ge 1$, where the numerical value of $j$ gives the order of the resonance.

When a particle is in an $m:(m-j)$ resonance, the characteristic angle $\phi$ – also called resonant angle – librates with an amplitude lower than 360°. The angle is given by

$$\phi = m\omega t - (m - j)\lambda_p - j\sigma, \quad (8)$$

where $\lambda_p$ is the mean longitude of the particle. For simplicity, we ignore variations in the mean longitude of epoch.

### 3.1 Resonance location

The angular and radial epicyclic frequencies are given by (Chandrasekhar 1942)

$$n^2 = \frac{1}{r} \frac{dU_0}{dr} \quad (9)$$

and

$$\kappa^2 = \frac{1}{r^3} \frac{d(r^4 n^2)}{dr}, \quad (10)$$

where $U_0$ is the axisymmetric part of the gravitational potential ($j = m = 0$).

From equation (5), we obtain

$$U_0 = -\frac{1}{r} - \frac{\mu}{2r} b_{1/2}^{(0)}(\alpha). \quad (11)$$

Expanding the Laplace coefficient up to second order in $\alpha$,

$$\frac{1}{2} b_{1/2}^{(0)}(\alpha) = 1 + \frac{1}{4}\alpha^2, \quad (12)$$

we obtain the axisymmetric part of the gravitational potential for the spherical body with a mass anomaly:

$$U_0 = -\frac{1}{r} \left( 1 + \mu + \frac{\mu}{4} \alpha^2 \right). \quad (13)$$

Keeping the lowest order terms in $\mu$ in equations (9) and (10), we obtain

$$n^2 = \frac{1}{r^3} \left( 1 + \mu + \frac{3\mu}{4} \alpha^2 \right) \quad (14)$$

and

$$\kappa^2 = \frac{1}{r^3} \left( 1 + \mu - \frac{3\mu}{4} \alpha^2 \right). \quad (15)$$

The location of the resonances can be obtained by numerical methods, such as the Newton–Raphson method (see Press et al. 1988; Renner & Sicardy 2006), by applying equations (14) and (15) in the resonance condition (equation 7). Table 2 shows the location of the resonances in the ranges $-4 \le m \le 4$ and $j \le 5$ (up to fifth-order resonances). The central body is a Chariklo-type body with $\mu = 10^{-3}$ and $\lambda = 0.471$ ($M = 6.3 \times 10^{18}$ kg and $T = 7.004$ h; Leiva et al. 2017), defined as our reference object. We assume $\mu = 10^{-3}$ as a reference value because it is small enough for the centre of the system to be approximately the physical centre of the spherical portion and large enough for the effects of the mass anomaly to be observed.
<table><caption>**Table 2.** The location of the *m*(*m* − *j*) resonances in the ranges −4 ≤ *m* ≤ 4 and *j* ≤ 5. We assumed a central body with parameters based on the centaur Chariklo, with *λ* = 0.471 and *µ* = 10⁻³ (reference object). The resonances marked ‘inside’ occur within the physical radius of the central body and, therefore, do not exist in the considered system. The resonances marked as ‘apsidal’ are out of the scope of this work.</caption><thead><tr><th>j</th><th>m →</th><th>-4</th><th>-3</th><th>-2</th><th>-1</th><th>1</th><th>2</th><th>3</th><th>4</th></tr></thead><tbody><tr><td rowspan="2">1</td><td>Resonance</td><td>5:6</td><td>4:5</td><td>3:4</td><td>2:3</td><td>1:2</td><td>1:0</td><td>2:1</td><td>3:2</td></tr><tr><td>a/R</td><td>1.909</td><td>1.993</td><td>2.156</td><td>2.612</td><td>Apsidal</td><td>Inside</td><td>1.256</td><td>1.358</td></tr><tr><td rowspan="2">2</td><td>Resonance</td><td>5:7</td><td>4:6</td><td>3:5</td><td>2:4</td><td>1:3</td><td>1:−1</td><td>2:0</td><td>3:1</td></tr><tr><td>a/R</td><td>2.156</td><td>2.313</td><td>2.612</td><td>3.423</td><td>Inside</td><td>Apsidal</td><td>Inside</td><td>Inside</td></tr><tr><td rowspan="2">3</td><td>Resonance</td><td>5:8</td><td>4:7</td><td>3:6</td><td>2:5</td><td>1:4</td><td>1:−2</td><td>2:−1</td><td>3:0</td></tr><tr><td>a/R</td><td>2.389</td><td>2.612</td><td>3.031</td><td>4.146</td><td>Inside</td><td>Inside</td><td>Apsidal</td><td>Inside</td></tr><tr><td rowspan="2">4</td><td>Resonance</td><td>4:8</td><td>3:7</td><td>2:6</td><td>1:5</td><td>1:−3</td><td>2:−2</td><td>3:−1</td><td>4:0</td></tr><tr><td>a/R</td><td>2.612</td><td>2.895</td><td>3.423</td><td>4.811</td><td>Inside</td><td>Inside</td><td>Inside</td><td>Apsidal</td></tr><tr><td rowspan="2">5</td><td>Resonance</td><td>4:9</td><td>3:8</td><td>2:7</td><td>1:6</td><td>1:−4</td><td>2:−3</td><td>3:−2</td><td>4:−1</td></tr><tr><td>a/R</td><td>2.825</td><td>3.164</td><td>3.793</td><td>5.433</td><td>Inside</td><td>Inside</td><td>Inside</td><td>Inside</td></tr></tbody></table>

## 3.2 Resonance width

In this subsection, we follow the classical approach of the pendulum model, presented in Winter & Murray (1997a) and Murray & Dermott (1999), to obtain the resonance width for our case of interest. A particle is in an *m*(*m* − *j*) sectoral resonance when its resonant angle *φ* librates, which means that the particle oscillates in the rotating frame around the central position of the resonance (equation 7). We can evaluate the maximum amplitude of a resonant particle through the temporal variations of *φ*:

$$ \dot{\phi} = m\omega - (m - j)n - j\dot{\omega} \quad (16) $$

and

$$ \ddot{\phi} = -(m - j)\dot{n} - j\ddot{\omega}. \quad (17) $$

Considering only the lowest order terms in eccentricity (*e*), we obtain using the Lagrange’s equations (Murray & Dermott 1999)

$$ \dot{n} = -3nC_r(m - j)e^j \sin \phi \quad (18) $$

and

$$ \dot{\omega} = je^{j-2}C_r \cos \phi, \quad (19) $$

where

$$ C_r = \mu \frac{n}{\alpha} \left[ \alpha F_j b_{1/2}^{(m-j)} + \frac{\lambda^2}{\alpha} f_j \delta_{|m|,1} \right]. \quad (20) $$

From equation (19), we obtain that the second derivative of *ω* is

$$ \ddot{\omega} = j(j-2)e^{j-3}\dot{e}C_r \cos\phi - je^{j-2}C_r \sin\phi\dot{\phi}, \quad (21) $$

where the time variation of eccentricity (*ė*) obtained through Lagrange’s equations is *ė* = −*je*<sup>*j*−1</sup>*C<sub>r</sub> sin *φ*.

It can be shown that

$$ \ddot{\omega} = j^2 e^{2(j-2)} C_r^2 \sin j\phi - je^{j-2}C_r(m\omega - (m - j)n)\sin\phi. \quad (22) $$

Therefore,

$$ \ddot{\phi} = -j^3 e^{2(j-2)} C_r^2 \sin 2\phi + 3nC_r(m-j)^2 e^j \sin \phi + j^2 e^{j-2} C_r(m\omega - (m-j)n) \sin \phi. \quad (23) $$

By inspection, we can evaluate the contribution of each term of equation (23). The *C*<sub>r</sub> function is proportional to *µ*, a value lower than one. In fact, for high values of mass anomaly (*µ* ≳ 10<sup>−2</sup>), we cannot assume the centre of mass of the system as the physical centre of the spherical object, and equation (5) is no longer applied – this range of *µ* defines another NSSB, the contact binary. Since *µ* < 1,

the term that depends on *C*<sub>r</sub> will dominate those dependent on *C*<sub>r</sub><sup>2</sup>, in principle. For first-order resonances (*j* = 1), the first and third terms in equation (23) are proportional to 1/*e*<sup>2</sup> and 1/*e*, respectively – *e* is a small value – and dominate over the second term, proportional to *e*.

### 3.2.1 Second and higher order resonances

For second and higher order resonances, the eccentricity exponents in equation (23) are positive, and we can approximate the equation to

$$ \ddot{\phi} + \omega_0^2 \sin \phi = 0, \quad (24) $$

where *ω*<sub>0</sub><sup>2</sup> = 3*n*|*C*<sub>r</sub>|(m − *j*)<sup>2</sup>*e*<sup>*j*</sup>. To obtain this result, we have assumed *m*ω − (m − *j*)*n* ≈ 0 since the particle is in resonance.

From equation (24), we can see that a resonant particle is confined in a pendulum motion around an equilibrium position of the resonance. The number of equilibrium positions of an *m*(*m* − *j*) sectoral resonance is *j*. Analogous to the simple pendulum problem, the particle reduced energy in the rotating frame is

$$ E = \frac{\dot{\phi}^2}{2} + 2\omega_0^2 \sin^2 \frac{\phi}{2}. \quad (25) $$

The maximum possible energy of the pendulum (*φ* = 0 deg and *φ* = 90 deg) defines the separatrix between libration and circulation of the resonant angle. That is, the separatrix corresponds to the boundary between bounded and unbounded motions. The energy of such trajectory is *E* = 2*ω*<sub>0</sub><sup>2</sup>, and the temporal variation of the resonant angle is *φ* = ±2*ω*<sub>0</sub> cos(*φ*/2).

Relating *φ* and *n*:

$$ dn = \frac{\dot{n}}{\dot{\phi}} d\phi = \pm \sqrt{3n|C_r|e^j} \sin \frac{\phi}{2} d\phi, \quad (26) $$

we obtain, by integration, the range of angular frequency in which a particle is in an *m*(*m* − *j*) sectoral resonance:

$$ n = n_0 \pm \sqrt{12n|C_r|e^j} \cos \frac{\phi}{2}, \quad (27) $$

where *n*<sub>0</sub> is the central angular frequency of the resonance.

Therefore, a particle is in a second or higher order resonance if its semimajor axis meets the relation:

$$ a = a_0 \pm \left( \frac{16 |C_r|}{3 n} e^j \right)^{1/2} a_0, \quad (28) $$

where *a*<sub>0</sub> is the central semimajor axis of the resonance (Section 3.1).
### 3.2.2 First-order resonances

For $m:(m-1)$ resonances, none of the terms in equation (23) can be disregarded, requiring a different solution than the one obtained. As ansatz, we assume a solution similar to equation (27), $n = n_0 + k\cos(\phi/2)$, where $k$ is an as-yet-unknown constant. By integrating equation (23), we obtain the kinetic energy of the system

$$ \frac{1}{2}\dot{\phi}^2 = \int \ddot{\phi} d\phi = \frac{C_r^2}{e^2} \left( 2\cos^2\frac{\phi}{2} + \cos^2\phi \right) + -6nC_r(m-1)^2e\cos^2\frac{\phi}{2} + \frac{4C_r}{3e}(m-1)k\cos^3\frac{\phi}{2}, \quad (29) $$

where the constant arising from the integration was determined considering $\phi = 0$ and 180 deg.

Applying $n = n_0 + k\cos(\phi/2)$ to equation (16) and assuming that the particle is exactly at the centre of the resonance ($\phi = 0$ and 180 deg), we find that $m\omega - (m-1)n_0 = -C_r/e$. From equation (16), we get

$$ \frac{1}{2}\dot{\phi}^2 = \frac{1}{2}\frac{C_r^2}{e^2}(1 + \cos\phi)^2 + \frac{1}{2}(m-1)^2k^2\cos^2\frac{\phi}{2} + \frac{C_r}{e}(1 + \cos\phi)(m-1)k\cos\frac{\phi}{2}. \quad (30) $$

Taking equations (29) and (30) as equal and assuming $\phi = 0$ deg:

$$ (m-1)^2k^2 + \frac{4C_r}{3e}(m-1)k + 12nC_r(m-1)^2e = 0. \quad (31) $$

Therefore, the boundaries of the angular frequency and semimajor axis in which a particle is in a first-order resonance are, respectively,

$$ n = n_0 \pm \sqrt{12|C_r|ne} \left( 1 + \frac{1}{27(m-1)^2e^3} \frac{|C_r|}{n} \right)^{1/2} - \frac{|C_r|}{3(m-1)e} \quad (32) $$

and

$$ a = a_0 \pm \left( \frac{16}{3} \frac{|C_r|}{n} e^j \right)^{1/2} \left( 1 + \frac{1}{27(m-1)^2e^3} \frac{|C_r|}{n} \right)^{1/2} a_0 + \frac{2}{9(m-1)e} \frac{|C_r|}{n} a_0. \quad (33) $$

## 4 POINCARÉ SURFACES OF SECTION

The Poincaré surface of section technique is usually applied in studies of the RP3BP (Hénon 1965a, b, 1966a, b, 1969; Jefferys 1971; Winter & Murray 1994a, b), to analyse the dynamics of the third body, providing information such as the location and size of stable and chaotic regions, including the mean motion resonance regions. In RP3BP, the problem is considered in a rotating system where the primary and secondary bodies are fixed, and only the third body describes a free motion. Some works have also adopted the Poincaré surface of section to study dynamical systems composed of two bodies, with a non-spherical central object. Scheeres et al. (1996) applied this technique to find periodic orbits around the asteroid 4769 Castalia. This technique was also applied by Borderes-Motta & Winter (2018) and Winter et al. (2019) to study the region around Toutatis and Haumea, respectively.

This work also applies the Poincaré surface of section to a two-body problem composed of a massive central body and a massless particle. Instead of the orbital motion between the primary and secondary bodies, the rotation of the central body gives the motion of the rotating frame. The Poincaré surface of

<figure><img src="image_5.png" /><figcaption>**Figure 2.** Poincaré surface of section for $C_J = 2.032 R^2\omega_k^2$ around an object with $\mu = 10^{-3}$ and $\lambda = 0.471$. The black islands are quasi-periodic orbits associated with the periodic orbit of first kind. Blue islands are associated with the 1:3 resonance and the green ones with the 2:7 resonance. The red points are chaotic orbits that cross the phase plane irregularly.</figcaption></figure>

section applied to the two-body problem with a mass anomaly provides information about stability and resonances. However, in this case, there are spin-orbit resonances instead of mean motion resonances.

The Poincaré surface of sections are maps generated in the phase space through the intersection points of the particle orbits with a fixed section in the system. These maps are generated for fixed values of the Jacobi constant (equation 6). In Fig. 2, we see an example of this map for a system composed of a massive central body with a mass anomaly. The Poincaré surface of section was defined in the plane $y = 0$ around our reference object and for the fixed value of the Jacobi constant $C_J = 2.032 R^2\omega_k^2$. We distributed the initial conditions on the x-axis.

In Fig. 2, the different sets of closed curves, called stability islands, delimit the stable regions of the system. Each stability island is formed by a single quasi-periodic orbit that is named because it does not have a defined orbital period. At the centre of the stability islands, we have periodic orbits. The latter crosses the Poincaré surface of section always at the same points and can be classified into two kinds (Poincaré 1895): those not associated with resonances are the first kind, and those associated with resonance are the second kind.

The point in the centre of all black closed curves is a first kind periodic orbit. In contrast, the points in the centres of the blue and green islands are the second kind orbits associated with the 1:3 and 2:7 resonances, respectively (Fig. 2). A single stability island identifies periodic orbits of first kind, while one or more stability islands can identify the orbits of second one. The number of islands for the second kind orbits is related to the order of the resonance (Winter & Murray 1997b). For example, the pair of blue islands in Fig. 2 is formed by quasi-periodic orbits that librate around the periodic orbit associated with the 1:3 resonance, a second-order resonance. In the same vein, each particle in the 2:7 resonance – a fifth-order resonance – generates five islands on the surface of section.
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Boundary curves between the chaotic (on the left) and stable (on the right) regions. The solid black line corresponds to the reference object, while the coloured solid and dashed lines are the cases in which we varied the parameters $\lambda$ and $\mu$, respectively.</figcaption></figure>

A set of first and second kind orbits belonging to the same resonance usually appears in a continuous Jacobi constant range and defines a family of orbits. The families present evolution in structure and position in the Poincaré surface of section during the variation of the Jacobi constant. In addition to the stability region delimited by the islands, there are also unstable regions filled with scattered red points in the figure, created by chaotic orbits. These chaotic regions are seen around the stability islands associated with the periodic orbits of 1:3 and 2:7 resonances. They are associated with the resonance separatrix and do not enter the stable regions, as we can see in Fig. 2. A stable region bounded by quasi-periodic orbits (black curves) separates the two chaotic regions.

In the following sections, we use the Poincaré surface of section to explore the stability around bodies with a mass anomaly by varying the central body parameters.

## 5 SYSTEM OVERVIEW

We studied the dynamics around our object by simulating a set of particles with pericentre distance $q$ and eccentricity in the ranges $1 < q/R \le q_f$ and $0 \le e \le 0.5$, respectively ($\Delta e = 0.05$ and $\Delta q/R = 0.01$). $q_f$ is a given value of $q$ for which all particles survive. The particles were simulated for 10 000 orbits. We assumed the parameters $\lambda$ and $\mu$ in the ranges $0.01 \le \lambda \le 1$ and $10^{-6} \le \mu \le 5 \times 10^{-3}$, respectively. Except for the near-Earth asteroids, we have that the vast majority of small heliocentric bodies have $\omega > \omega_k$ (Warner, Harris & Pravec 2009), justifying the fact that we do not focus on cases with $\lambda > 1$.

We verified in all numerical simulations the existence of a chaotic region just outside the central body in which particles collide or are ejected from the system. Beyond the chaotic region, there is a stable region, and the boundaries between them are shown in Fig. 3. Particles with semimajor axis and eccentricity in the region bounded by the curve (on the left-hand side of the figure) will be lost, while those outside the boundary will survive for at least 10 000 orbits. The solid black line corresponds to our reference object, while the solid coloured and dashed lines provide the boundary curves for systems where we vary $\lambda$ and $\mu$, respectively.

The successive close encounters of the particle with the mass anomaly are responsible for exchanges of energy and angular momentum, resulting in the variation of the particles' orbital elements. Particles with sufficiently small semimajor axis show orbital evolution with chaotic diffusive character (for details, see Rollin et al. 2021). In general, eccentricities in the chaotic region tend to increase, resulting in occasional collisions or until the orbit becomes hyperbolic.

Fig. 4 shows the trajectory in the rotating frame (a) and eccentricity (b) of a pair of particles initially in circular orbits around a central

<figure><img src="image_7.png" /><figcaption>**Figure 4.** (a) Trajectory in the rotating frame and (b) temporal evolution of the eccentricity. The innermost particle (red line) is at $1.74 R$ and the outermost one (blue line) at $3.48 R$, and both are initially in circular orbits. The parameters of the central body are $\lambda = 0.471$ and $\mu = 5 \times 10^{-3}$.</figcaption></figure>

body with parameters $\lambda = 0.471$ and $\mu = 5 \times 10^{-3}$. The semimajor axes of the innermost (red line) and outermost (blue line) particles are $1.74 R$ and $3.48 R$, respectively. We observe that the eccentricity shows a secular increase for the innermost particle, reaching values up to 0.15. The particle collides with the central body after about $3.5\omega_k^{-1}$ or $\sim 12$ spin periods. The eccentricity shows periodic variations for the outermost particle, and the particle remains stable around the central body.

The boundary curves are robust against the final simulation period and are preserved when we extend the simulations to 100 000 orbits. Lages et al. (2017) analyse through the Lyapunov exponent the stability of particles around a contact binary, obtaining a chaotic region coherent with ours. Our boundary curve is also coherent with the region where the particles are lost in the numerical simulations for a Chariklo with a mass anomaly performed by Sicardy et al. (2019).

The boundary between chaotic and stable regions has only a slight dependence on the relative mass anomaly. Although the increase of $\mu$ generates only a small swell of the chaotic region, it produces larger increments in eccentricity and the particles are lost more quickly. The extension of the chaotic region is mainly affected by the parameter $\lambda$. Decreasing in 10 times the rotating rate, we obtain that the chaotic region is more than doubled, a result also obtained by Lages et al. (2017).

In order to crudely evaluate the extension of the chaotic region, we calculate for a set of systems the semimajor axis at which a particle in circular orbit will survive for up to 10 000 orbits – the threshold
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Threshold semimajor axis obtained in selected numerical simulations (markers) and through equation (34) (solid lines). The x-axis gives the normalized mass anomalies, and the different colours and markers give the rotating rates.</figcaption></figure>

semimajor axis ($a_t$). The curve fitted from the numerical results is given by

$$ \frac{a_t}{R} = [1.298 - 0.007M + 0.006M^2 + 0.674\lambda^{-0.75}], \quad (34) $$

where $M = -\log \mu$. Physically, we can interpret the threshold semimajor axis as the minimum semimajor axis, beyond which rings and satellites can exist around a body with a mass anomaly.

Fig. 5 shows the threshold semimajor axis obtained in the numerical simulations. The x-axis gives the normalized mass anomaly of each simulation, and different colours and markers show the different rotating rates. The solid lines correspond to the curves given by equation (34) (the colour of the lines matches the colour of the markers for the same $\lambda$).

Fig. 6 shows the position and width of the sectoral resonances, obtained theoretically (Section 3), where each panel corresponds to a different normalized mass anomaly, while the rotating rate is fixed as $\lambda = 0.471$. The vertical black line in each panel corresponds to the corotation radius $a_c$ of the system ($a_c = \lambda^{-2/3}$), and the coloured lines give the boundaries of the resonances. The white region on the left provides the initial conditions of particles with pericentre within the central body, and the white area on the right is the stable region. The grey area places the chaotic region.

Since the sectoral resonances are spin-orbit resonances, we have that $\mu$ has a minor effect on their locations, as seen in the figure. However, the resonance width will depend on $\mu$, as an increase in the mass anomaly will enhance the gravitational perturbation felt by the particles, allowing larger regions to be connected to the resonance equilibrium points. As we increase the numerical value of $m$, the resonances approach the corotation radius. The first-order resonances with $|m| > 4$ overlap for the case with $\mu = 10^{-4}$ (Fig. 6a). For $\mu = 10^{-3}$ (Fig. 6b), we see that additional first-order resonances, such as 4:5 and 3:4, overlap for high eccentricities, while for $\mu = 5 \times 10^{-3}$ (Fig. 6c) the overlap is intensified, covering the 2:3 resonance.

The overlap of first-order resonances is generally responsible for eliminating stable regions associated with the resonance (Wisdom 1980; Winter & Murray 1997a, b). Thus, they should contribute to the chaotic behaviour verified in the systems. It is not by chance that the region surrounding the corotation radius is always chaotic. However, while the overlap helps to carve the chaotic region around the central object, it is not the primary source of chaoticity for the system. Such a fact can be seen in 7a, where the chaotic region covers a region where there is no overlap of first-order resonances. As already mentioned, encounters with the mass anomaly produce chaotic diffusion of the orbits, clearing an entire region that extends beyond the corotation radius. In our numerical simulations, we did not find stable particles in internal resonances.

Analogously to Fig. 6, we present in Fig. 7 the resonances and chaotic and stable regions around an object with a mass anomaly,

now keeping $\mu$ constant and varying $\lambda$. As we can see in the figure, the parameter affects both the location and the width of the resonances (equations 7 and 20). When $\lambda$ is increased by one order, resonances move more than four times closer to the body, while chaotic regions approach it by only a factor of 2. So, by changing $\lambda$, we change which resonances will be in the stable region.

In the case shown in Fig. 7(c), the rotation frequency is equal to the Keplerian one, which places the corotation radius on the surface of the spherical portion of the body. Consequently, the internal resonances and some external ones reside within the central body, with only a few resonances in the stable region. Assuming objects with even faster rotation, we get a narrower chaotic region with fewer resonances outside the object, corresponding to unattractive cases. In the hypothetical case where the spin frequency tends to infinity, it would be not sectoral resonances or chaotic region since it falls in the case of a non-rotating spherical object with a ridge at its equator.

In Section 6, we analyse the evolution of the stable region and the external resonances using the Poincaré surfaces of section.

## 6 STABLE REGION

In Section 5, we have shown the existence of two distinct regions around a spherical body with a mass anomaly: a chaotic region where particles collide or are ejected, and a stable region, which will be our focus in this section. First, we compare the resonance widths obtained by numerical simulations with those predicted by the analytical model described in Section 3. We then analyse the motion of test particles in the vicinity of external resonances.

We put our analytical model to the test using the following methodology: (i) for a given resonance, we theoretically calculate its central position (Section 3.1) and the Jacobi constant (equation 6), initially assuming a circular orbit; (ii) we generate the Poincaré surface of section of a broad region around the central position; (iii) by visual inspection, we obtain the position of the stable fixed point of the resonance – which corresponds to the real central position of the resonance – and the limits of the widest island surrounding the point – the width of the resonance; and (iv) then, we successively increase the eccentricity by $10^{-2}$ and repeat the previous steps until the islands disappear or until we reach $e = 0.5$.

Fig. 8 shows the resonance widths obtained theoretically and through the Poincaré surfaces of sections, for our reference object and an object with $\lambda = 0.157$ and $\mu = 10^{-3}$. We found that numerical data agree reasonably well with the analytical model, indicating that the pendulum model with necessary adaptations applies to our system. In general, we obtain that the largest divergences occur for larger eccentricities ($e \gtrsim 0.2$). It is expected, since we assumed first-order approximations in eccentricity in the development of the pendulum model.

We verify that the innermost resonances present the largest displacements in the central position for the reference case. These same displacements are verified for the case with $\lambda = 0.157$ (Fig. 8b) in which resonances are at least twice as far from the central body. As a rule, we obtain that displacements depend on the distance $a_t/a_c$ from the resonance to the corotation radius. The central positions we obtained differ by less than 5 per cent from those theoretically obtained, demonstrating the robustness of the analytical method.

After the validity of the analytical model is attested, we turn our attention to the resonance dynamics. Fig. 9(a) shows the Poincaré surface of section of the region around the 2:5 resonance for the reference object and $C_1 = 1.964 R^2 \omega_k^2$. In the figure, we can identify four types of motion in the stable region: periodic motion of first kind, quasi-periodic motion associated with the latter, resonant and chaotic
<figure><img src="image_9.png" /><figcaption>**Figure 6.** Semimajor axis ($a/R$) versus eccentricity ($e$) for systems, with $\lambda = 0.471$ and (a) $\mu = 10^{-4}$, (b) $\mu = 10^{-3}$, and (c) $\mu = 5 \times 10^{-3}$. Particles with initial $a/R$ and initial $e$ in the left white region have pericentre within the central body and collide. Particles in the grey area collide with the central body or are ejected, and those in the right white one remain in the system for more than 10 000 orbits. The dashed black lines correspond to the corotation radius, and the coloured lines provide the theoretical boundaries of the resonances. Coloured lines not referenced on the label and close to the corotation radius correspond to first-order resonances with $|m| > 4$.</figcaption></figure>

motion – which is not stable despite being within the region defined by us as stable. The single black dot at $x_0/R = 3$ corresponds to the orbit classified as periodic of first kind. Periodic orbits in a Poincaré surface of section divide the x-axis positions into pericentric positions – at smaller $x$ – and apocentric positions – at larger $x$. Seeing the right part of the figure, we have particles with higher initial eccentricity, forming black closed curves surrounding the periodic orbit. They are quasi-periodic orbits and define regions where particles remain indefinitely in stable motion without other effects.

The orange islands correspond to orbits associated with the 2:5 resonance, where every single dot in the centre of an island is a

stable fixed point of the resonance. All three orange dots in the figure correspond to a single periodic orbit of second kind. Due to energy exchanges between the central body and resonant particles, the latter can remain stable even in the presence of other effects, depending on the system conditions (Madeira & Giuliatti Winter, 2020).

The red dots show the chaotic zone between the resonant and quasi-periodic orbits. Appendix A shows the Poincaré surfaces of section of the resonances found around the reference object. As one can see, the chaotic zones at the resonances separatrices are always narrow, showing that the region we named the ‘stable region’ has in fact a few very small strips of confined chaotic motion. Moving
<figure><img src="image_10.png" /><figcaption>**Figure 7.** Semimajor axis ($a/R$) versus eccentricity ($e$) of systems with $\mu = 10^{-3}$ and rotating rate $\lambda$ given in the caption of each panel. Particles with initial conditions in the white region on the left have pericentre within the central body and collide, while those in the grey area show chaotic behaviour. The white region on the right is the stable region. The dashed black line provides the corotation radius, and the coloured lines give the theoretical boundaries of the resonances. Coloured lines not referenced on the label and close to the corotation radius correspond to first-order resonances with $|m| > 4$.</figcaption></figure>

to the right in Fig. 9(a), there is another region with quasi-periodic orbits that extends up to $x_0/R = 3.57$. After this limit, we have the chaotic region, where the red dot at $x_0/R = 3.65$ corresponds to a particle that collides with the central body.

In Fig. 9(b), we present the evolution of the 2:5 resonance, showing the largest stable orbit of the resonance, an intermediate one, and the central orbit, for different values of the Jacobi constant. Since the 2:5 resonance is a third-order resonance ($j = 3$), each initial condition produces three distinct islands in the Poincaré surfaces of section. Except for the case with $C_J = 1.959 R^2 \omega_k^2$, we see that the three islands shrink and get closer as the value of $C_J$ increases. The Jacobi constant and the eccentricity are inversely proportional, so the

latter decreases from right to left in the figure. The resonance width decreases with the eccentricity, and the resonant orbits tend towards the periodic orbit of the first kind, explaining why the islands shrink until they disappear.

To understand why the largest red island is smaller than the largest blue island (Fig. 9b), we present in Fig. 10(a) the Poincaré surfaces of section for $C_J = 1.959 R^2 \omega_k^2$, and in Fig. 10(b), the theoretical and numerical boundaries of the 2:5 resonance. The red dashed line places the case with $C_J = 1.959 R^2 \omega_k^2$ and the grey area is the chaotic region. For this value of Jacobi constant, the resonance is at the edge between the stable and chaotic regions. As a result, the particles most bounded to the resonance – closer to the stable fixed point – remain
<figure><img src="image_11.png" /><figcaption>**Figure 8.** The width of the external sectoral resonances in the stable region for (a) an object with $\lambda = 0.471$ and $\mu = 10^{-3}$ and (b) for an object with $\lambda = 0.157$ and $\mu = 10^{-3}$. The solid and dashed lines give the widths predicted by the analytical model, and the coloured filled regions delimit the obtained numerically widths. The grey region corresponds to the chaotic region near the central body.</figcaption></figure>

stable (in orange), while particles closer to the resonance boundaries initially follow the pattern expected for resonant particles. However, they are showing the stickiness phenomenon behaviour mimicking the resonant behaviour, but they are lost from the system at some point.

The eccentricity of one of these less bounded particles is shown in Fig. 10(c) by the solid red line, while the eccentricity of the central resonant orbit is the solid orange line. Both particles show a periodic variation in eccentricity. However, the eccentricity of the less bounded particle also shows an increase, reaching $e \sim 0.145$.

Fig. 9(c) shows, in the rotating frame, the periodic orbit of the second kind seen in Fig. 9(a), where the colour coding gives the velocity. Since the particle is at the stable fixed point of the resonance, the orbit is closed. Also, the orbit is retrograde ($v < 0$) because the resonance is beyond the corotation radius. As the central body is symmetric, there will always be at least one axis that divides the orbit into two symmetric parts. For example, for the orbit shown in Fig. 9(c), this axis corresponds to $y = 0$.

Sicardy (2020) discusses some additional symmetries expected for the trajectory of a particle in an $m:(m-j)$ sectoral resonance. The orbit is invariant by a rotation of $360/|m|$ deg, and it has a total of $|m|(j-1)$ self-crossing. For the 2:5 case ($j=3$ and $m=-2$), we see that the orbit is invariant by a rotation of 180 deg and has four self-crossings.

A peculiarity of our system is the positions of the particle pericentre and apocentre. In RP3BP, in which the disturbing body is

at $x_0/R = 1$ and the particle is initially at $x_0/R < 1$, the gravitational effect felt by the particle is weaker (stronger) when it is on the x-axis at $x/R > 0$ ($x/R < 0$). Consequently, the particle starts at the pericentre, the apocentre being in the opposite direction. In our case, we have the opposite scenario. The orbit position where a particle feels the strongest gravitational effect is on the x-axis at $x/R > 0$ – where the modulus of gravitational force is the sum of the forces of the two portions of the central body. Thus, a particle initially on the x-axis ($x/R > 0$) starts at its apocentre (minimum velocity), as we can see from the dot labelled ‘1’ in Fig. 9.

Fig. 11 shows, from top to bottom, the Poincaré surface of section for $C_j = 2.087 R^2 \omega_c^2$, the whole evolution of the islands of the 1:4 resonance, and the trajectory of a particle at a stable fixed point of the resonance. As shown in Fig. 11(a), the overview of the resonance neighbourhood is similar to the 2:5 resonance, with a narrow, chaotic region at the resonance boundaries, surrounded by a region with periodic/quasi-periodic orbits. A crucial difference, however, is obtained in the resonance islands. While the 2:5 resonance has three stable fixed points, we obtained in Fig. 11(a) six stable points for the 1:4 resonance.

To understand the dynamics of the resonance, we colour green (Fig. 11a) the trajectory of a particle near one of these points. The particle is responsible for forming three islands around three of the stable fixed points (in green in the figure). This fact leads us to conclude that the resonance is the 1:4 (of third-order) and not the 2:8 (of sixth-order) as we would obtain in the ellipsoidal problem
<figure><img src="image_12.png" /><figcaption>**Figure 9.** (a) Poincaré surface of section for $C_J = 1.964 R^2 \omega_k^2$, with $\lambda = 0.471$ and $\mu = 10^{-3}$. We assume initial conditions with $3.15 \le x_0/R \le 3.84$. The black curves are the periodic and quasi-periodic orbits of first kind, and the orange curves are orbits associated with the 2:5 resonance. Red dots correspond to chaotic orbits. (b) Evolution of the 2:5 resonance islands, where the colours of the dots correspond to values of $C_J$ given on the figure's label. (c) Central orbit of the 2:5 resonance for $C_J = 1.964 R^2 \omega_k^2$ in the rotating frame. The temporal evolution of the orbit is given by numbers and dots equally spaced in time, while the colour coding gives the velocity in the rotating frame.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Figure 10.** (a) Poincaré surface of section for $C_J = 1.959 R^2 \omega_k^2$ where the periodic/quasi-periodic orbits of first kind are in black, the 2:5 resonance islands are in orange, and the particles in the chaotic region are in red. (b) Theoretical boundaries of the 2:5 resonance are shown by the solid orange lines. In contrast, the filled orange and grey regions are regions numerically obtained for the 2:5 resonance and the chaotic region, respectively. The red dashed line gives the initial conditions of the simulation with $C_J = 1.959 R^2 \omega_k^2$. (c) Eccentricity of a pair of particles: the one that remains in the system is orange, and the unstable one is red.</figcaption></figure>

– which we will address in a following publication. The islands produced by the particle have the particularity of being asymmetric in relation to the *x*-axis – we say that the particle is in asymmetric libration (Beaugé 1994; Winter & Murray 1997b).

In Fig. 11(b), we highlight the Poincaré surface of section islands of some particles by plotting them in black, intending to show the asymmetric libration. Each island produced by a particle in asymmetric libration has a mirror image obtained from the motion of
<figure><img src="image_14.png" /><figcaption>**Figure 11.** (a) Poincaré surface of section for $C_J = 2.087 R^2 \omega_k^2$, with $\lambda = 0.471$ and $\mu = 10^{-3}$. We assumed initial conditions with $3.70 \le x_0/R \le 5.97$ and separated the distinct types of orbits by colour: the periodic/quasi-periodic orbits of first kind are in black, the 1:4 resonance orbits are in purple and green, and chaotic ones in red. (b) Resonance islands for different values of $C_J$. The label on the panel gives the colour of the largest island for each value of $C_J$. (c) Central orbit in the rotating frame of one of the families associated with the 1:4 resonance (in green in the top panel) for $C_J = 2.087 R^2 \omega_k^2$. The numbers and colours on the panel provide time evolution and velocity in the rotating frame, respectively.</figcaption></figure>

a different particle in asymmetric libration. Closer to the resonance boundaries, we also get ‘horseshoe fashion’ orbits encompassing pairs of fixed points of two different trajectories.

When we refer to asymmetric resonance or libration, we refer to the symmetry of the trajectory in the Poincaré surface of section and not in the *xy*-plane. As already mentioned, the trajectory in the *xy*-plane of the resonant particles has a symmetry axis due to the symmetric mass distribution in the central body. For example, in Fig. 11(c), the axis of symmetry would correspond to the axis connecting the point ‘6’ to the centre of the system.

Several works such as Message (1970), Frangakis (1973a, b), Message & Taylor (1978), Bruno (1994), Beaugé (1994), and Winter & Murray (1997b) have studied asymmetric periodic orbits in the context of RP3BP, showing that these orbits are characteristics of 1:1+p resonances and are obtained only for particles with eccentricities above a critical value. Similar to particles in *m*(m-j) resonance with m≠-1, the ones with eccentricity lower than this threshold value present symmetric libration in the Poincaré surface of section. We obtained these same results for the case with mass

anomaly. In Fig. 11(b), the critical eccentricity is reached somewhere between the Jacobi constants $2.122 R^2 \omega_k^2$ and $2.152 R^2 \omega_k^2$. Carrying out a set of Poincaré surface of section in this interval, we obtain that the critical eccentricity for the 1:4 resonance is $e \sim 0.167 (2.136 R^2 \omega_k^2)$.

Fig. 12(a) shows one island of the 1:4 resonance for the critical eccentricity ($2.136 R^2 \omega_k^2$, in green), for $2.133 R^2 \omega_k^2$ (in purple), and $2.139 R^2 \omega_k^2$ (in blue). For the highest value of Jacobi constant (smallest eccentricity), we see a single stable point in the figure related to a single family of resonant orbits. The critical eccentricity is reached by decreasing the Jacobi constant, and the stable point bifurcates into two points (the stars in the figure). Each of the points gives rise to an independent family of resonant orbits. The *x*-axis, which previously allocated the single stable point, now allocates the unstable equilibrium point after the bifurcation, corresponding to the inflexion position of the ‘horseshoe fashion’ orbits.

Fig. 12(b) shows the trajectories of the stable points given by stars in the top panel. The orbits are mirror versions of each other. The same is obtained for eccentricities in Fig. 12(c), in which the
<figure><img src="image_15.png" /></figure>

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

<figure><img src="image_17.png" /><figcaption>**Figure 12.** (a) Poincaré surface of section of one island of the 1:4 resonance for $C_J = 2.133, 2.136$, and $2.139 R^2\omega_k^2$ (in purple, green, and blue, respectively). The black and red stars are the stable points obtained after bifurcation. (b) Trajectories and (c) eccentricities of the stable points given by stars in the top panel, where the colour of the solid lines coincides with the colour of the star for the same stable point.</figcaption></figure>

red curve is the mirror version of the black one with respect to time $t \approx 4.2\omega_k^{-1}$ (pericentre passage time). As discussed in Bruno (1994), the bifurcation of the stable points is related to the indirect term of the disturbing function, which differs from zero only for $1: 1 + p$ resonances (equation 5).

We show another example of particles in asymmetric libration in Fig. 13. From top to bottom, this figure shows the Poincaré surface of section of the region of 1:2 resonance for a central body with $\lambda = 0.157$ and $\mu = 10^{-3}$, the whole evolution of the resonance, and the

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

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

<figure><img src="image_20.png" /><figcaption>**Figure 13.** (a) Poincaré surface of section for $C_J = 0.915 R^2\omega_k^2$, with $\lambda = 0.157$ and $\mu = 10^{-3}$. The non-resonant orbits are in black. Particles in 1:2 resonance and chaotic orbits are in orange and green and red, respectively. (b) Poincaré surface of section for some particles in 1:2 resonance, with $C_J = 0.907, 0.912, 0.915$, and $0.918 R^2\omega_k^2$. Different colours of the islands involved by the same 'horseshoe fashion' orbit correspond to different particles. (c) Trajectory of a stable fixed point shown in orange in the top panel, where the colour coding gives the velocity and the numbers and dots, the time evolution of the orbit.</figcaption></figure>
<figure><img src="image_21.png" /><figcaption>**Figure 14.** Semimajor axis versus eccentricity for Chariklo system, where coloured lines place the sectoral resonances, and the grey area corresponds to the chaotic region. Coloured lines not referenced on the label, between 1.4 and 1.8, correspond to first-order resonances with $|m| > 4$. A vertical dashed line at $a/R \approx 1.7$ gives the corotation radius, while vertical dotted lines give the central location of the rings.</figcaption></figure>

trajectory of a particle in a stable fixed point of the resonance. As in the reference case, chaotic behaviour is seen only in a narrow region in the separatrices, with a large regular region of periodic/quasi-periodic orbits of first kind around the resonances. For the 1:2 resonance, we have a low critical eccentricity ($e \sim 10^{-2}$), with symmetric libration only in the cases where the resonance islands are tiny.

Trajectories of particles in 1:2 resonance are the only ones without self-crossings in the rotating frame, as we can see in the bottom panel. Such fact has implications for the temporal evolution of a ring of particles, as self-crosses increase collisions between particles. In this context, a ring with particles into the 1:2 resonance or in periodic/quasi-periodic orbits of first kind – which do not show self-crossing either – should have a lower rate of collisions than a ring with particles in other resonances, disregarding other external effects. The particles shown in Figs 11(c) and 13(c) do not start at the apocentre of the orbit because they are not initially with $\dot{x} = 0$. A particle around a body with mass anomaly will start at its apocentre only when that condition is met.

## 7 APPLICATION TO THE CHARIKLO SYSTEM

Leiva et al. (2017), using stellar occultation data, investigated the shape of Chariklo, obtaining four distinct shape models for the object: a sphere, a MacLaurin spheroid, a triaxial ellipsoid, and a Jacobi ellipsoid. According to Sicardy et al. (2019), observational data suggest the presence of topographic features of typical heights of 5 km in the spherical solution. This fact places Chariklo as a possible body with a mass anomaly. In this section, we briefly study the dynamics around Chariklo, in particular in the region of the ring. The rings have orbital radii of 391 and 405 km, with radial widths of 7 and 3 km, respectively (Bérard et al. 2017).

We performed numerical simulations adopting the spherical Chariklo given by Leiva et al. (2017), $\lambda = 0.471$, with a mass anomaly of $\mu = 7 \times 10^{-6} = [5 \text{ km}/(2 \times 129 \text{ km})]^3$. Fig. 14 shows the width of the resonances and the location of the chaotic region. The vertical dashed line gives the corotation radius and the central location of the rings by the vertical dotted lines. We obtained a threshold semimajor axis of $a_t/R = 2.5$ in the numerical simulation. This result is in good agreement with our adjusted equation (equation 34) that returns $a_t/R = 2.6$. The 1:2 resonance is the only first-order resonance beyond the chaotic region. The region beyond the chaotic one is essentially stable, hosting the rings and possibly moons, depending on their eccentricity.

Leiva et al. (2017) propose that the inner ring is associated with the 1:3 spin-orbit resonance. Therefore, we studied this resonance in detail, as it is close to both rings. Fig. 15(a) shows the Poincaré surface of section for the largest Jacobi constant obtained by us for the 1:3 resonance ($C_J = 2.038 R^2 \omega_k^2$). For this value of $C_J$, the resonance has not reached the critical eccentricity, and we obtain only a single symmetric periodic orbit. A narrow, chaotic region surrounds the islands of resonance, but the whole set is surrounded by a stable region associated with periodic/quasi-periodic orbits of first kind.

Fig. 15(b) shows the motion in the rotating frame of the periodic orbits given in Fig. 15(a), in which the colour of the orbits matches those given in the top panel, and the green regions correspond to the location of the rings. For clarity, we only show the portions of the orbits with $y > 0$. We show the radial variation of the orbits in Fig. 15(c).

The 1:3 resonance orbit has one self-crossing at $y = 0$ and a period of almost $6.4\omega_k^{-1}$. In contrast, the trajectory of the first kind follows the Chariklo shape, with a period of almost $3.2\omega_k^{-1}$. As one can see in the figure, both orbits are initially in the inner ring – near its outer edge. However, only the first kind of periodic orbit remains within the ring throughout the simulation, while the resonant orbit crosses the ring edges and reaches the outer ring.

The difference in radial variation is due to the different nature of the orbits. Periodic orbits of first kind correspond to nearly circular orbits, while those of second kind are intrinsically eccentric, explaining why the latter has a significantly larger radial variation. Here, when we refer to eccentricity, we are referring to osculating elements defined in the context of the classical two-body problem. We refer the reader to the work of Ribeiro et al. (2021) for a detailed discussion regarding the orbital elements in the context of NSSBs.

The results discussed in the last paragraph lead us to question whether the inner ring is associated with the 1:3 resonance. To verify this, we present in Fig. 16(a) diagram with the semimajor axis versus eccentricity for a range of values corresponding to the rings (in green). We also show the orbital elements obtained for the periodic orbits of first and second kinds (black and blue lines, respectively). The region filled in blue shows the 1:3 resonance boundaries obtained in the Poincaré surface of sections.

The largest possible eccentricity for a particle to remain within the boundaries of the inner ring is $e = 9 \times 10^{-3}$, which is smaller than the smallest eccentricity obtained for the resonant orbits ($e = 10^{-2}$). In addition, the resonant periodic orbits and the ring are displaced, indicating that the ring is not confined by such resonance.
<figure><img src="image_22.png" /><figcaption>**Figure 15.** (a) Poincaré surface of section of the Chariklo system for $C_J = 2.038 R^2 \omega_k^2$. We show different orbits by different colours: the non-resonance orbits are black, the 1:3 resonant orbits are blue, and the chaotic ones are red. (b) Motion in the rotating frame for $y > 0$ and (c) radial variation of periodic orbits shown in panel (a). The orbits of the first and second kind are given by black and blue lines, respectively, and the green regions correspond to the positions of Chariklo rings.</figcaption></figure>

Meanwhile, periodic orbits of first kind cover a broad region and encompass both rings. The entire region shown in Fig. 16, which is not associated with the 1:3 resonance (blue region), is composed of periodic/quasi-periodic orbits of first kind, including the ring region. Therefore, we conclude that Chariklo rings are associated with first kind orbits and not with the 1:3 resonance, as proposed by Leiva et al. (2017). Similar results were obtained by Winter et al. (2019) for the Haumea ring.

<figure><img src="image_23.png" /><figcaption>**Figure 16.** Diagram of the semimajor axis versus eccentricity. The green regions show the range of values that corresponds to the location of the rings. The blue line shows the orbital elements obtained for the central orbit of the 1:3 resonance, and the blue filled regions give the boundaries of the resonance. The black line gives the periodic orbits of first kind.</figcaption></figure>

study, varying the parameters of the central object. We can summarize our overall results as follows:

(i) The pendulum model with the necessary adaptations and the Poincaré surface of section proved to apply to the mass anomaly problem. We verified a strong agreement between the results by comparing both techniques.

(ii) There is a chaotic region near the central object where particles collide or are ejected due to chaotic diffusion caused by successive close encounters with the mass anomaly. Mysen et al. (2006), Mysen & Aksnes (2007), and Lages et al. (2017) also obtain chaotic regions near the central object for elongated bodies and contact binaries, respectively.

(iii) For the set of parameters analysed by us, the chaotic region extends beyond the corotation radius. This fact indicates a lack of stable internal sectoral and corotation resonances in the mass anomaly system.

(iv) Resonances location is mainly affected by the mass of the spherical portion and the spin period. In contrast, the masses of the spherical and anomalous portions of the body and spin period are responsible for determining the width of the resonances.

(v) Beyond the chaotic region, there is a region where the motion of the particles is dynamically stable. In such a region, there is chaotic behaviour only in a narrow region in the separatrices of the resonances.

(vi) The behaviour of the particles in the external sectoral resonances is similar to those obtained for the mean motion resonances in the RP3BP (Winter & Murray 1997a, b). Similar to RP3BP, we verify the existence of asymmetric periodic orbits associated with $1:1+p$ resonances.

Although objects with the shape assumed in this work are unknown so far, the completely irregular shapes known for some asteroids lead us to speculate that such a class of object might exist. We emphasize that bodies with mass anomaly are perfectly reasonable outputs from a collision of a satellite that spirals towards the central body due to tidal dissipation or a collision between two objects at low velocity, with partial accretion (Leinhardt & Stewart 2011).

Sicardy et al. (2019) discuss the possibility of Chariklo having a spherical shape with topographic feature with $\mu \sim 10^{-5}$, which places the Centaur as a first candidate to integrate the class of mass anomaly objects. We studied the dynamics around a Chariklo with mass anomaly and found that 1:3 resonant particles present radial variations too large for the radial extension of the inner ring. On the other hand, particles in periodic/quasi-periodic first kind orbits show radial motions that match the extension of the two rings of Chariklo. Consequently, the ring must be associated with these orbits and not

## 8 CONCLUSIONS AND DISCUSSION

In this paper, we have attempted to perform a general analysis of the dynamics of particles around a spherical body with a mass anomaly. For this, we used well-known techniques of the three-body problem
with orbits of second kind as proposed by Leiva et al. (2017). With the constant increase in data on small heliocentric bodies, we believe that objects with shapes similar to bodies with mass anomalies may soon be detected.

It is essential to point out that in this study, we limited ourselves to analyse the dynamics of an isolated particle around a NSSB, disregarding the effects associated with the ring particles, such as collisions between them, local viscous, and self-gravity effects. We also disregard external effects that modulate the dynamics of small particles, such as solar radiation pressure and Poynting-Robertson drag. Nevertheless, the location and width of resonances and the chaotic region are general results and should remain almost unchanged in the presence of other effects. Therefore, our work presents some tools and first general results for studies on dynamics of mass anomaly systems.

## ACKNOWLEDGEMENTS

This study was financed in part by the Coordenação de Aperfeiçoamento de Pessoal de Nível Superior (CAPES, Finance Code 001), Fundação de Amparo à Pesquisa do Estado de São Paulo (FAPESP, Proc. 2016/24561-0, Proc. 2018/23568-6, and Proc. 2021/07181-7), and Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq, Proc. 305210/2018-1 and Proc. 313043/2020-5). Finally, we thank the anonymous reviewer for the comments that significantly improved our work.

## DATA AVAILABILITY

The data underlying this article will be shared on reasonable request to the corresponding author.

## REFERENCES

Amarante A., Winter O. C., 2020, MNRAS, 496, 4154

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

Bérard D. et al., 2017, AJ, 154, 144

Borderes-Motta G., Winter O. C., 2018, MNRAS, 474, 2452

Braga-Ribas F. et al., 2014, Nature, 508, 72

Bruno A. D., 1994, The Restricted 3-Body Problem: Plane Periodic Orbits, Vol. 17. Walter de Gruyter, Berlin

Chandrasekhar S., 1942, Principles of Stellar Dynamics. The University of Chicago Press, Chicago, Il

Chapman C. et al., 1995, Nature, 374, 783

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

Frangakis C., 1973a, Ap&SS, 22, 421

Frangakis C., 1973b, Ap&SS, 23, 17

Geissler P., Petit J.-M., Durda D. D., Greenberg R., Bottke W., Nolan M., Moore J., 1996, Icarus, 120, 140

Hanuš J., Marchis F., Æurech J., 2013, Icarus, 226, 1045

Hanuš J., Marchis F., Viikinkoski M., Yang B., Kaasalainen M., 2017, A&A, 599, A36

Hénon M., 1965a, Ann. Astrophys., 28, 499

Hénon M., 1965b, Ann. Astrophys., 28, 992

Hénon M., 1966a, Ann. Astrophys., 29, 49

Hénon M., 1966b, Ann. Astrophys., 29, 57

Hénon M., 1969, C. R. Acad. Sci., Paris Ser. B Sci. Phys., 269, 223

Hudson R. S., Ostro S. J., 1995, Science, 270, 84

Jefferys W. H., 1971, An Atlas of surfaces of section for the restricted problem of three bodies. Univ. Texas Austin, Austin, TX

Jorda L. et al., 2016, Icarus, 277, 257

Lages J., Shepelyansky D. L., Shevchenko I. I., 2017, AJ, 153, 272

Lages J., Shevchenko I. I., Rollin G., 2018, Icarus, 307, 391

Lauretta D. et al., 2017, Space Sci. Rev., 212, 925

Leinhardt Z. M., Stewart S. T., 2011, ApJ, 745, 79

Leiva R. et al., 2017, AJ, 154, 159

Madeira G., Giuliani Winter S.M., 2020, The European Physical Journal Special Topics, 229, 1527

Meiss J., 1992, Rev. Mod. Phys., 64, 795

Merline W. J., Weidenschilling S. J., Durda D. D., Margot J.-L., Pravec P., Storrs A. D., 2002, in Bottke W. F., Jr, Cellino A., Paolicchi P., Binzel R. P., eds, Asteroids III, Vol. 1. Univ. Arizona Press, Tucson, AZ, p. 289

Message P. J., 1970, in Giacaglia G. E. O., ed., Periodic Orbits Stability and Resonances. D. Reidel Publ. Co., Dordrecht, Holland, p. 19

Message P. J., Taylor D. B., 1978, On Asymmetric Periodic Solutions of the Plane Restricted Problem of Three Bodies, and Bifurcations of Families. D. Reidel Publ. Co., Dordrecht, Holland, p. 319

Moura T., Winter O., Amarante A., Sfair R., Borderes-Motta G., Valvano G., 2020, MNRAS, 491, 3120

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

Mysen E., Aksnes K., 2007, A&A, 470, 1193

Mysen E., Olsen Ø., Aksnes K., 2006, Planet. Space Sci., 54, 750

Ortiz J. L. et al., 2017, Nature, 550, 219

Poincaré H., 1895, Analysis situs. Gauthier-Villars et Fils, Paris, p. 3

Press W. H., Flannery B. P., Teukolsky S. A., Vetterling W. T., 1988, Numerical Recipes in C: The Art of Scientific Computing. Cambridge University Press, Cambridge, UK

Prockter L., Murchie S., Cheng A., Krimigis S., Farquhar R., Santo A., Trombka J., 2002, Acta Astronaut., 51, 491

Ragozzine D., Brown M. E., 2009, AJ, 137, 4766

Renner S., Sicardy B., 2006, Celest. Mech. Dyn. Astron., 94, 237

Ribeiro T., Winter O., Mourão D., Boldrin L., Carvalho J., 2021, MNRAS, 506, 3068

Rollin G., Shevchenko I. I., Lages J., 2021, Icarus, 357, 114178

Scheeres D. J., 1994, Icarus, 110, 225

Scheeres D. J., Ostro S. J., Hudson R. S., Werner R. A., 1996, Icarus, 121, 67

Scheeres D., Williams B. G., Miller J. K., 2000, J. Guid. Control Dyn., 23, 466

Shang H., Wu X., Cui P., 2015, Ap&SS, 355, 69

Shevchenko I. I., 2011, New Astron., 16, 94

Sicardy B., 2020, AJ, 159, 102

Sicardy B., Leiva R., Renner S., Roques F., El Moutamid M., Santos-Sanz P., Desmars J., 2019, Nat. Astron., 3, 146

Sicardy B., Renner S., Leiva R., Roques F., El Moutamid M., Santos-Sanz P., Desmars J., 2020, in Prialnik D., ed., The dynamics of rings around Centaurs and Trans-Neptunian Objects. Elsevier, Amsterdam, Netherlands, p. 249

Vantieghem S., 2014, Proc. R. Soc. A: Math. Phys. Eng. Sci., 470, 20140093

Wang X., Li J., Gong S., 2016, MNRAS, 455, 3724

Warner B. D., Harris A. W., Pravec P., 2009, Icarus, 202, 134

Werner R. A., 1994, Celest. Mech. Dyn. Astron., 59, 253

Winter O. C., Borderes-Motta G., Ribeiro T., 2019, MNRAS, 484, 3765

Winter O. C., Murray C. D., 1994a, Atlas of the Planar, Circular, Restricted Three-Body Problem. I. Internal Orbits. No. 1, QMW Maths Notes. Queen Mary and Westfield College, London

Winter O. C., Murray C. D., 1994b, Atlas of the Planar, Circular, Restricted Three-Body Problem. II. External Orbits. No. 2, QMW Maths Notes. Queen Mary and Westfield College, London

Winter O. C., Murray C. D., 1997a, A&A, 319, 290

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

Winter O. C., Valvano G., Moura T. S., Borderes-Motta G., Amarante A., Sfair R., 2020, MNRAS, 492, 4437

Wisdom J., 1980, AJ, 85, 1122

Yoshikawa M., Kawaguchi J., Fujiwara A., Tsuchiyama A., 2015, in Michel P., DeMeo F. E., Bottke W. F., eds, Asteroids IV, Vol. 1. Univ. Arizona Press, Tucson, AZ, p. 397

Yu Y., Baoyin H., 2012, AJ, 143, 62
## APPENDIX A: POINCARÉ SURFACES OF SECTION FOR THE REFERENCE OBJECT

In this appendix, we present the Poincaré surface of sections for all resonances in the stable region of the reference case – $\lambda = 0.471$ and $\mu = 10^{-3}$. All figures have three panels, following the pattern:

a) The top panel shows the Poincaré surface of section of a broad region, with the value of $C_J$ given in the upper right corner of the figure. Black closed curves are periodic/quasi-periodic orbits of first kind. In red is the chaotic motion. Closed coloured curves correspond to islands of the resonance given on the figure label. In the case of 1: 1 + *p* resonances, we plot some islands with different colours to show the asymmetric libration of the resonance.

b) The middle panel shows the entire evolution of the resonance. The label at the top gives the colour of the largest island for each $C_J$ value. In the case of 1: 1 + *p* resonances, we plotted some islands in black to highlight the asymmetric libration observed in such resonances.

c) The bottom panel shows the trajectory of the central orbit of the resonance in the rotating frame. The value of $C_J$ is the same as the top panel. We show the temporal evolution of the orbit by numbers and dots equally spaced in time, while colour coding gives the velocity in the rotating frame.

<figure><img src="image_24.png" /><figcaption>**Figure A1.** (a) Poincaré surface of section for $C_J = 1.937 R^2 \omega^2$ in which the black curves are the periodic and quasi-periodic orbits of first kind and the blue curves are orbits associated with the 4:9 resonance. (b) Central orbit of the 4:9 resonance for $C_J = 1.937 R^2 \omega^2$ in the rotating frame. The temporal evolution of the orbit is given by numbers and dots equally spaced in time, while the colour coding gives the velocity in the rotating frame.</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>**Figure A2.** (a) Poincaré surface of section for $C_J = 1.948 R^2 \omega^2$ in which the periodic/quasi-periodic orbits of first kind are in black and the 3:7 resonance orbits are in green. (b) Resonance islands for different values of $C_J$ given in the upper panel. (c) Central orbit of the 3:7 resonance for $C_J = 1.948 R^2 \omega^2$. The numbers and colours on the panel provide time evolution and the velocity in the rotating frame, respectively.</figcaption></figure>
<figure><img src="image_26.png" /><figcaption>**Figure A3.** (a) Poincaré surface of section for $C_J = 1.964 R^2\omega^2$. The black curves are periodic and quasi-periodic orbits of first kind, and the orange curves are orbits associated with the 2:5 resonance. Red dots correspond to chaotic orbits. (b) Evolution of the 2:5 resonance islands, where the colours of the dots correspond to the values of $C_J$ given on the label of the figure. (c) Central orbit of the 2:5 resonance for $C_J = 1.964 R^2\omega^2$ in the rotating frame. The temporal evolution of the orbit is given by numbers and dots equally spaced in time, while the colour coding gives the velocity in the rotating frame.</figcaption></figure>

<figure><img src="image_27.png" /><figcaption>**Figure A4.** (a) Poincaré surface of section for $C_J = 1.986 R^2\omega^2$. Black curves are periodic and quasi-periodic orbits of first kind, and blue curves are orbits associated with the 3:8 resonance. Red dots correspond to chaotic orbits. (b) Evolution of the 3:8 resonance islands, where the colours of the dots correspond to values of $C_J$ given on the label of the figure. (c) Central orbit of the 3:8 resonance for $C_J = 1.986 R^2\omega^2$ in the rotating frame. The temporal evolution of the orbit is given by numbers and dots equally spaced in time, while the colour coding gives the velocity in the rotating frame.</figcaption></figure>

<figure><img src="image_28.png" /><figcaption>**Figure A5.** (a) Poincaré surface of section for $C_J = 2.024 R^2\omega^2$. Periodic/quasi-periodic orbits of first kind are in black, the 1:3 resonance orbits are in blue and orange, and chaotic ones in red. (b) Resonance islands for different values of $C_J$. The label on the panel gives the colour of the largest island for each value of $C_J$. (c) Central orbit in the rotating frame of one of the families associated with the 1:3 resonance (in blue in the top panel) for $C_J = 2.024 R^2\omega^2$. The numbers and colours on the panel provide time evolution and the velocity in the rotating frame, respectively.</figcaption></figure>

The pattern above was not followed only for the 4:9 resonance, for which we do not show the middle panel. We got only one value of $C_J$ in the stable region for this resonance. Figures are given from

the resonance closest to the central body to the farthest. Figs A1–A9 correspond to 4:9, 3:7, 2:5, 3:8, 1:3, 2:7, 1:4, 1:5, and 1:6 resonances, respectively.
<figure><img src="image_29.png" /><figcaption>**Figure A6.** (a) Poincaré surface of section for $C_J = 2.073 R^2 \omega_k^2$. Black curves are periodic and quasi-periodic orbits of first kind, and green curves are orbits associated with the 2:7 resonance. Red dots correspond to chaotic orbits. (b) Evolution of the 2:7 resonance islands, where the colours of the dots correspond to the values of $C_J$ given on the label of the figure. (c) Central orbit of the 2:7 resonance for $C_J = 2.073 R^2 \omega_k^2$ in the rotating frame. The temporal evolution of the orbit is given by numbers and dots equally spaced in time, while the colour coding gives the velocity in the rotating frame.</figcaption></figure>

<figure><img src="image_30.png" /><figcaption>**Figure A7.** (a) Poincaré surface of section for $C_J = 2.087 R^2 \omega_k^2$. Periodic/quasi-periodic orbits of first kind are in black, the 1:4 resonance orbits are in purple and green, and chaotic ones in red. (b) Resonance islands for different values of $C_J$. The label on the panel gives the colour of the largest island for each value of $C_J$. (c) Central orbit in the rotating frame of one of the families associated with the 1:4 resonance (in green in the top panel) for $C_J = 2.087 R^2 \omega_k^2$. The numbers and colours on the panel provide time evolution and velocity in the rotating frame, respectively.</figcaption></figure>

<figure><img src="image_31.png" /><figcaption>**Figure A8.** (a) Poincaré surface of section for $C_J = 2.147 R^2 \omega_k^2$. The periodic/quasi-periodic orbits of first kind are in black, the 1:5 resonance orbits are in brown and green, and chaotic ones in red. (b) Resonance islands for different values of $C_J$. The label on the panel gives the colour of the largest island for each value of $C_J$. (c) Central orbit in the rotating frame of one of the families associated with the 1:5 resonance (in brown in the top panel) for $C_J = 2.147 R^2 \omega_k^2$. The numbers and colours on the panel provide time evolution and the velocity in the rotating frame, respectively.</figcaption></figure>
<figure><img src="image_32.png" /><figcaption>**Figure A9.** (a) Poincaré surface of section for $C_J = 2.204 R^2 \omega_k^2$. The periodic/quasi-periodic orbits of first kind are in black, the 1:6 resonance orbits are in cyan and blue, and chaotic ones in red. (b) Resonance islands for different values of $C_J$. The label on the panel gives the colour of the largest island for each value of $C_J$. (c) Central orbit in the rotating frame of one of the families associated with the 1:6 resonance (in cyan in the top panel) for $C_J = 2.204 R^2 \omega_k^2$. The numbers and colours on the panel provide time evolution and the velocity in the rotating frame, respectively.</figcaption></figure>

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