# Chaotic zones around rotating small bodies

J. Lages, Dima Shepelyansky, Ivan I. Shevchenko

► To cite this version:

J. Lages, Dima Shepelyansky, Ivan I. Shevchenko. Chaotic zones around rotating small bodies. The Astronomical Journal, 2017, 153, pp.272. 10.3847/1538-3881/aa7203. hal-01519663

HAL Id: hal-01519663

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

Submitted on 1 Sep 2025

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

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

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

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

# Chaotic Zones around Rotating Small Bodies

José Lages¹, Dima L. Shepelyansky², and Ivan I. Shevchenko¹,³,⁴

¹ Institut UTINAM, Observatoire des Sciences de l'Univers THETA, CNRS, Université de Franche-Comté,
Besançon F-25030, France; jose.lages@utinam.cnrs.fr

² Laboratoire de Physique Théorique du CNRS, IRSAMC, Université de Toulouse, UPS, Toulouse F-31062, France

³ Pulkovo Observatory, RAS, 196140 Saint Petersburg, Russia

⁴ Lebedev Physical Institute, RAS, 119991 Moscow, Russia

Received 2017 March 9; revised 2017 May 2; accepted 2017 May 7; published 2017 June 1

## Abstract

Small bodies of the solar system, like asteroids, trans-Neptunian objects, cometary nuclei, and planetary satellites, with diameters smaller than 1000 km usually have irregular shapes, often resembling dumb-bells or contact binaries. The spinning of such a gravitating dumb-bell creates around it a zone of chaotic orbits. We determine its extent analytically and numerically. We find that the chaotic zone swells significantly if the rotation rate is decreased; in particular, the zone swells more than twice if the rotation rate is decreased 10 times with respect to the “centrifugal breakup” threshold. We illustrate the properties of the chaotic orbital zones in examples of the global orbital dynamics about asteroid 243 Ida (which has a moon, Dactyl, orbiting near the edge of the chaotic zone) and asteroid 25143 Itokawa.

**Key words:** celestial mechanics – chaos – comets: general – minor planets, asteroids: general – planets and satellites: dynamical evolution and stability

## 1. Introduction

The orbital dynamics around irregularly shaped bodies (having complex gravity fields) was extensively studied in the last two decades. The reason is twofold: first, satellites of small bodies such as asteroids were discovered; second, space missions were planned and accomplished to asteroids and cometary nuclei. Therefore, many aspects of the orbital dynamics in rotating complex gravity fields were studied, both theoretically and in numerical simulations; see Scheeres (1994, 2012) and references therein. Small bodies of the solar system (asteroids, trans-Neptunian objects, cometary nuclei, and planetary satellites) with diameters of less than 1000 km usually have strongly irregular shapes (Melnikov & Shevchenko 2010; Jorda et al. 2016), in many cases resembling dumb-bells or “contact binaries.” Various models for gravity fields of the “central body” were used: that of a triaxial ellipsoid with uniform density (Chauvineau et al. 1993; Mysen et al. 2006; Olsen 2006; Mysen & Aksnes 2007), a rod (Bartczak & Breiter 2003), a dumb-bell or “bilobed” model (Marchis et al. 2014; Feng et al. 2016), a collection (“molecule”) of gravitating points (Petit et al. 1997), a polyhedral model (Werner 1994; Werner & Scheeres 1996), and a truncated gravitational field derived from a shape model (Feng et al. 2017). Orbits around actual small bodies, such as asteroids Castalia, Eros, and Hektor, were extensively modeled (Scheeres et al. 1996, 2000; Yu & Baoyin 2012; Marchis et al. 2014). Concerning the dumb-bell model, it was also used in the problem of spin-spin resonances in a system of two aspherical gravitating bodies (Hut 1981; Batygin & Morbidelli 2015): the quadrupole moment of the secondary was represented as a dumb-bell of two equal masses. This model provides a setting for a qualitative description of the tidal evolution and the resulting spin-spin coupling of tight binary systems of elongated bodies (Batygin & Morbidelli 2015).

Many studies were devoted to resonant phenomena and the determination of orbital stability regions; see, in particular, Scheeres (1994), Hu & Scheeres (2004), Mysen et al. (2006), Olsen (2006), Mysen & Aksnes (2007), Scheeres (2012), and

references therein. The existence of “chaotic gravitational zones” around rotating elongated bodies was outlined by Mysen et al. (2006). A destabilizing role of resonances between a particle’s orbital motion and the rotational motion of the central elongated body was revealed by Mysen et al. (2006), Olsen (2006), and Mysen & Aksnes (2007); in particular, see Figures 1–5 in Olsen (2006) and Figures 2 and 9 in Mysen & Aksnes (2007), where the integer spin-orbit resonances form a characteristic “saw of instability” in the plane of initial values of the semimajor axis and eccentricity of the orbiting particle. Quite recently, numerical simulations of orbits around contact binaries were performed by Feng et al. (2016) emphasizing the stabilization effect of the fast rotation of the contact binary: for faster binary’s rotation, periodic orbits start to lose stability closer to the barycenter; this is explained by averaging of the perturbation.

The preceding theoretical studies were based on the analysis of perturbation functions and their expansions, in particular, resonant terms in the expansions. In our article, we present a different technique, based on the analysis of symplectic maps (see, e.g., Meiss 1992), in particular, a generalized Kepler map. The Kepler map approach allows one to understand straight-forwardly the global orbital behavior. Let us recall that the Kepler map is a two-dimensional area-preserving map, describing the eccentric circumbinary motion of a massless particle in the gravitational field of a primary and a perturber (the secondary moving around the primary in a circular orbit deeply inside the particle’s orbit). The motion is described in terms of changes in the particle’s energy and conjugated orbital phase measured at its apocenter and pericenter passages. In particular, it was shown that the Kepler map describes the dynamics of highly eccentric comets (Petrosky 1986; Malyshkin & Tremaine 1999), Comet Halley among them (Chirikov & Vecheslavov 1989). In an appropriate physical model, it explains the phenomenon of strong microwave ionization of excited hydrogen atoms (Casati et al. 1988) and autoionization of molecular Rydberg states (Benvenuto et al. 1994). A review of the Kepler map theory in a historical context is given in
<figure><img src="image_3.png" /><figcaption>**Figure 1.** We model a non-axisymmetric small body (e.g., here 25143 Itokawa Gaskell et al. 2008) by a contact binary $m_1 \ge m_2$ ($M = m_1 + m_2$). The size of the dumb-bell is *d*, the small body center of mass is marked by the cross. The axis of rotation is perpendicular to the figure plane and passes through the center of mass.</figcaption></figure>

Shevchenko (2011). Recent applications of the Kepler map theory along with its corresponding advancements concern processes of disintegration of three-body systems and Levý flight statistics in these processes (Shevchenko 2010), capture of dark matter by the solar system and by binary stellar systems (Lages & Shepelyansky 2013; Rollin et al. 2015), accurate symplectic map description of the long-term dynamics of Comet Halley (Rollin et al. 2015). In this article, the Kepler map is used mostly for analytical purposes, thus to provide an analytical description of resonances and borders of dynamical chaos in the stability diagrams. However, it is also used as a numerical tool, whose advantage is in the enormously high speed of computation, which allows one to construct the stability diagrams with very high spatial resolution (see Section 3, Figures 4, 5, and 7).

In our work, we consider a passively gravitating particle orbiting a gravitating dumb-bell. If the dumb-bell is fixed in space, the particle cannot gain or lose orbital energy or angular momentum for its orbital motion, because their source is absent. But if the dumb-bell rotates, the particle's energy or angular momentum may vary strongly, so that the particle may even escape or fall on the primary, depending on initial conditions. Obviously, one expects that the particles close to the primary are more prone to such disturbances than those away from it.

It is already known that a gravitating binary, such as a binary star or a binary asteroid, has a circumbinary chaotic zone, where all circumbinary orbits of the orbiting particles with any initial eccentricity are chaotic (Shevchenko 2015). However, what would be the case if one considers the motion around a rigid dumb-bell, for which the spinning frequency $\omega$ can be smaller or larger than the Keplerian frequency $\omega_0$ fixed by Kepler's third law? Here we give an answer to this question generalizing the Kepler map description (Petrosky 1986; Chirikov & Vecheslavov 1989) to describe the motion of a particle in the gravitational field of a spinning body modeled by a dumb-bell with masses $m_1, m_2$ separated by constant distance (dumb-bell size) *d*. In such a way, we model an irregular body by two contact uniform-density spheres (equivalent to two

point masses) as is shown in Figure 1 for an example of asteroid 25143 Itokawa (Gaskell et al. 2008). The dumb-bell is spinning around its center of mass with an angular frequency $\omega$, which can be different from the Keplerian frequency $\omega_0$ of revolution of masses $m_1, m_2$. The dynamics of particles orbiting the dumb-bell is considered in the plane orthogonal to the spin axis.

The Kepler map description of orbits about a spinning non-axisymmetric body is achieved by introducing a parameter, $\omega$, which is the rate of rotation of the model contact binary (see Figure 1). The value of $\omega$ is arbitrary. We derive analytical expressions for the kick function, representing the energy increment for the test particle when it passes the apocenter of its orbit. We consider the planar case, i.e., the case of the orbits lying in the plane orthogonal to the small body spin axis. We note that the Kepler map appears also for molecular Rydberg states with a rotating dipole core (Benvenuto et al. 1994). In the gravitational potential, the dipole term cancels, and in the dumb-bell case the quadrupole and octupole contributions of the central body's gravitational field provide leading terms in the kick function. However, we show that, in a wide range of spinning frequencies $\omega < \omega_0$, retaining the quadrupole term is enough to qualitatively describe the chaotic zone around the spinning body. Strikingly, such a zone swells significantly for $\omega < \omega_0$ down to a certain threshold. In our approach, we derive the kick function in a closed form, valid in the whole range of parameters' values. To connect our theoretical findings with observational data, we illustrate the properties of the chaotic orbital zones in examples of the global orbital dynamics about asteroid 243 Ida (which has a moon, Dactyl, orbiting near the edge of the chaotic zone) and asteroid 25143 Itokawa.

## 2. The Kepler Map Description

We consider the motion of a passively gravitating particle in the planar circular restricted three-body problem $m_1-m_2$-particle, where the two masses $m_1$ and $m_2$ are connected by a massless rigid rod, thus forming a dumb-bell (see Figure 1). The Keplerian rate of rotation of a contact binary, i.e., two tangent spheres, is

$$ \omega_0 = \sqrt{\pi G \rho / 3}, \qquad (1) $$

where $\rho$ is the density of the irregular body (Scheeres 2007). For a typical density $\rho = 1 \text{ g/cm}^3$, we have $\omega_0 = 2.5 \times 10^{-4} \text{ s}^{-1}$ corresponding to a period of about 7 hr. There are many observed asteroids with significantly larger rotation periods (see, e.g., Pravec et al. 2008). From now on, we express the physical quantities in the following units: $GM = 1$ (where $M = m_1 + m_2$ is the total mass of the irregular gravitating body, we choose $m_2 \le m_1$ and we define $\mu = m_2/M \le 0.5$), $d = 1$ is the size of the effective dumb-bell (Figure 1), and the Keplerian frequency $\omega_0 = \sqrt{G(m_1 + m_2)/d^3} = 1$; particle's energy per unit of mass, $E$, is then expressed in units of $d^2\omega_0^2$. We consider solely the case of prograde (with respect to the dumb-bell rotation) orbits of the particle; analysis of the retrograde case is analogous. The Kepler map for the motion around a gravitating dumb-bell, if one allows for the arbitrary rotation rate $\omega$ of the dumb-bell, takes the form (Casati et al. 1988; Benvenuto et al. 1994)

$$ E_{i+1} = E_i + \Delta E(\phi_i), \quad \phi_{i+1} = \phi_i + 2\pi\omega/|2E_{i+1}|^{3/2} \quad (2) $$
where the subscript *i* enumerates the pericenter passages with the rotation phase $\phi_i = \omega t_i$ and the corresponding particle energy $E_i$ taken at apocenter. We retrieve the original Kepler map derived in Chirikov & Vecheslavov (1989) and Petrosky (1986) by setting $\omega = \omega_0 = 1$. The equation for the rotation phase $\phi_i$ variation is given by the third Kepler law. Originally, the map has been derived for particle's orbit with a period larger than the period of the binary, i.e., $\omega/|2E|^{3/2} \ge 1$ (Casati et al. 1988; Benvenuto et al. 1994). We obtain the kick function $\Delta E$ generalizing to the case of a binary with an arbitrary non-Keplerian rotation velocity ($\omega \ne \omega_0$) the work of Roy & Haddow (2003) and Heggie (1975) devoted to energy change in hard binary due to distant encounters. Defining $\mathbf{r} = \cos(\omega t)\hat{\mathbf{a}} + \sin(\omega t)\hat{\mathbf{b}}$ is the position of the dumb-bell lobe of mass $m_2$ relative to the dumb-bell lobe of mass $m_1$ ($\hat{\mathbf{a}}$ and $\hat{\mathbf{b}}$ are any two orthogonal fixed directions of the plane), and $\mathbf{R}$ is the position of the test particle relative to the barycenter of the two lobes, the equation of motion for the test particle around the dumb-bell is given by

$$ \ddot{\mathbf{R}} = - \nabla_{\mathbf{R}} \Phi(\mathbf{R}, \mathbf{r}, \mu, \omega), \quad (3) $$

where the gravitational potential reads

$$ \Phi(\mathbf{R}, \mathbf{r}, \mu, \omega) = - \frac{1 - \mu}{\|\mathbf{R} + \mu\mathbf{r}\|} - \frac{\mu}{\|\mathbf{R} - (1 - \mu)\mathbf{r}\|}. \quad (4) $$

Defining $\mathbf{r} = \|\mathbf{r}\|$ and $\mathbf{R} = \|\mathbf{R}\|$, the multipole expansion of the gravitational potential gives

$$ \begin{aligned} \Phi(\mathbf{R}, \mathbf{r}, \mu, \omega) = & -\frac{1}{R} - \mu(1 - \mu)\frac{r^2}{2R^3} \\ & \times \left(3\left(\frac{\mathbf{r} \cdot \mathbf{R}}{rR}\right)^2 - 1\right) - \mu(1 - \mu)(2\mu - 1) \\ & \times \frac{r^3}{2R^4}\left(5\left(\frac{\mathbf{r} \cdot \mathbf{R}}{rR}\right)^3 - 3\frac{\mathbf{r} \cdot \mathbf{R}}{rR}\right) + O\left(\frac{r^4}{R^5}\right). \end{aligned} \quad (5) $$

Here, besides the $1/R$ term, the first two leading terms of the series are retained. This turns out to be well enough for the purposes of the present analysis, as comparisons of our results with previous simulations show (see Sections 4 and 5). The energy increment

$$ \Delta E(\mu, q, \omega, \phi) = - \int_{-\infty}^{+\infty} \dot{\mathbf{R}} \cdot \nabla \left( \Phi + \frac{1}{R} \right) dt \quad (6) $$

of a test particle forced to follow a parabola, the focus of which is the dumb-bell barycenter, is a function of the pericenter distance $q$, and of the phase of the dumb-bell $\phi$ when the test particle passes at pericenter. Here, the two lobes of the rotating dumb-bell form a circular binary. Following Roy & Haddow (2003) in the case of a circular binary but rotating at an arbitrary frequency rate $\omega$, keeping the two first leading terms for the kick function (6) we obtain

$$ \Delta E(\mu, q, \omega, \phi) \approx W_1 \sin(\phi) + W_2 \sin(2\phi). \quad (7) $$

In Equation (7), the exchange of energy between the small spinning body and the test particle is split in two terms: the first harmonic comes from the octupole term $(\propto r^3/R^4)$ of the

gravitational potential multipole expansion (5) with amplitude

$$ W_1 \approx \mu(1 - \mu)(1 - 2\mu)2^{1/4}\pi^{1/2}\omega^{5/2}q^{-1/4} \\ \times \exp\left(-\frac{2^{3/2}}{3}\omega q^{3/2}\right) \quad (8) $$

and the second harmonic comes from the quadrupole term $(\propto r^2/R^3)$ with amplitude

$$ W_2 \approx -\mu(1 - \mu)2^{15/4}\pi^{1/2}\omega^{5/2}q^{3/4} \exp\left(-\frac{2^{5/2}}{3}\omega q^{3/2}\right). \quad (9) $$

We note that expression (8), restricted to the case $\omega = \omega_0 = 1$ and $\mu \ll 1$, has been obtained using a different method in Shevchenko (2011). Usually, in the Kepler map ($\omega = \omega_0$) the kick function $\Delta E$ is proportional to $\sin\phi$, which is just the first most prominent term in the Fourier expansion of the energy increment, especially if $\mu \ll 1$ (Petrosky 1986; Shevchenko 2011). This is, for example, the case when one considers the Kepler map description of cometary dynamics around the solar system modelized by the Sun and Jupiter as perturber (Petrosky 1986; Chirikov & Vecheslavov 1989). However, with an increase of $\mu$, the second harmonic $(\propto \sin 2\phi)$ becomes more and more important, and even remains the sole term for the case $\mu = 1/2$ since the first harmonic $(\propto \sin \phi)$ disappears ($W_1 = 0$). Indeed, for $\mu = 1/2$, due to the equality of the mass of primaries, by symmetry, the perturbation frequency is effectively doubled.

Here, for the case of spinning small bodies, a wide range of rotation frequencies can be considered; in particular, spinning frequencies for asteroids range from $\omega/2\pi \sim 10^{-3} \text{ h}^{-1}$ to $\omega/2\pi \sim 10^2 \text{ h}^{-1}$ (Whiteley et al. 2002; Warner et al. 2009; Hergenrother & Whiteley 2011). For $q \gg (\omega_0/\omega)^{2/3}d$, the contribution $W_1$ is obviously dominant since a factor of two exists between the arguments of the exponentials entering Equations (8) and (9). This absolute prominence of $W_1$ over $W_2$ is even quadratically shifted farther from the small body for $\omega < \omega_0$. Conversely, which contribution, either $W_1$ or $W_2$, dominates is not so obvious for the region $q \lesssim (\omega_0/\omega)^{2/3}d$, which for $\omega < \omega_0$ encompasses the immediate vicinity of the spinning small body. The two contributions $W_1$ (8) and $W_2$ (9) depend on the parameters $\mu$, $\omega$, and $q$; their relative importance is summarized in the $(\omega, q)$ plot for different values of $\mu$ (Figure 2). We clearly see that below the frequency of disruption of a rubble-pile object ($\omega < \omega_0$), for any mass parameter $\mu$, the quadrupole coefficient $W_2$ generally dominates the octupole coefficient $W_1$ in the vicinity of the spinning small body. For example, $W_2 \gg W_1$ for $q \lesssim 3d$, $\omega \lesssim \omega_0$, and for any $\mu$ parameter.

Typical amplitudes of energy kick functions $\Delta E$ are shown in Figure 3. Analytical curves (7) constructed using the first (8) and the second (9) harmonic terms of the multipole expansion of the dumb-bell gravitational potential are in good agreement with kick energy $\Delta E$ obtained by direct integration of Newton's equations (Figure 3). Globally, the decrease of the small body spinning frequency induces an increase of the energy kick. As expected for $q = 7d$ insignificant kick ($\Delta E \sim 10^{-8} d^2\omega_0^2$) is expected in the case of an ordinary binary rotating with $\omega = \omega_0$. However, in the case of a spinning small body at, e.g., $\omega = 0.1\omega_0$, the energy kick is strongly enhanced ($\Delta E \sim 10^{-2}d^2\omega_0^2$). In comparison with an ordinary binary, such an energy kick increase induced by
<figure><img src="image_4.png" /><figcaption>Figure 2. Contours of the function $W_1(\omega, q)/W_2(\omega, q)$ for $\mu \to 0$ (left panel), $\mu = 0.4$ (middle panel), and $\mu = 0.49$ (right panel).</figcaption></figure>

a slow spinning frequency allows the zone of chaos to extend quite far from the central body. In Figure 3 (left panel), amplitudes of kick functions $\Delta E$ are presented divided by the mass factor $\mu(1 - \mu)$ entering the expression of $W_2(9)$. For $q/d = 3, 5, 7$, we clearly see that below $\omega/\omega_0 \simeq 1, 0.5, 0.3$, curves for any reduced mass $\mu$ are superimposed stressing again the fact that the second harmonic term is dominant for small spinning frequencies (see also Figure 3, right panel).

It should be noted that upon a minor modification this study can be applied to a more generalized body, namely to a planar molecule representing a set of coplanar asymmetric dumb-bells of various sizes and $\mu$ with a common center of mass. In this way, the Kepler map is straightforwardly generalized by means of adding separate terms corresponding to each elementary dumb-bell's contribution in the equation for the energy increment; each added term has its own amplitude and constant phase shift in the body's orientation.

In the frame of 3D atoms in a monochromatic field in 3D a symplectic map was shown to give a correct description of real dynamics (Casati et al. 1988). However, for a rotating gravitating body, the generalization of our dumb-bell Kepler map to the 3D case is an analytically complicated task, as a 3D generalization of the classical Kepler map by Emelyanenko (1990) shows. We reserve this for a future study.

### 3. Stability Diagrams and Central Chaotic Zone

Stability diagrams are constructed by computing Lyapunov exponents on a fine grid of initial data, $(e, q)$ or $(e, a)$. Lyapunov exponents are computed iterating concurrently the dumb-bell Kepler map (2) and its tangent map (as, e.g., described by Chirikov 1979 in application to the standard map). The motion is regarded as chaotic, if the maximum Lyapunov exponent is non-zero and positive. Such diagrams are presented in the $(q, e)$ plane for $\mu = 1/2$ and for different values of $\omega/\omega_0 = 0.068, 0.4$, and 1 (Figure 4). The border delimiting chaotic domain (red) from regular domain (blue) is ragged; the most prominent teeth being associated with the integer $p:1$ and half-integer $p + 1/2:1$ resonances between particle orbital frequency and small body spinning frequency. Here any neighboring integer and half-integer resonances are equally sized due to the symmetry of the dumb-bell for $\mu = 1/2$, indeed half-period and full-period rotations of the symmetric dumb-bell both result in configurations identical to

the initial one. The stability diagram graphically demonstrates how the integer and half-integer resonances overlap. Let us define the central chaotic zone as the zone in $q$ such as at any initial eccentricity the particle's dynamics is chaotic. Otherwise, the chaotic zone is defined as the region where even particles initially in circular orbits become dynamically chaotic. From Figure 4, we clearly see that the central chaotic zone swells significantly as the small body spinning frequency decreases, since its farthest extent varies from $q \simeq 2.8d$ for $\omega = \omega_0$ to $q \simeq 7d$ for $\omega \simeq 0.068\omega_0$.

Based on the concept of the chaotic layer around the separatrix and using analytical expressions for the classical Kepler map parameter, a strictly analytical expression for the size of the central chaotic zone around a gravitating binary can be derived (Shevchenko 2015). In a similar way, the size of the central chaotic zone around a rotating gravitating dumb-bell can be analytically estimated. Let us retain in (7) only the second harmonic contribution, since $W_2$ clearly dominates over $W_1$ for small spinning frequencies ($\omega < \omega_0$), indeed from Equations (8) and (9), for $\omega \ll (q/d)^{-3/2}\omega_0 \ll \omega_0$, we obtain $W_2/W_1 \approx 2^{7/2}q/(1 - 2\mu)$, which diverges as $\mu$ approaches $1/2$. By the substitution $E = W_2 y$ and $\phi = x/2$, the map (2) is reduced to

$$y_{i+1} = y_i + \sin x_i, \quad x_{i+1} = x_i + \lambda/|y_{i+1}|^{3/2} \quad (10)$$

with $\lambda = 2^{1/2}\pi\omega/|W_2|^{3/2}$. Following the standard procedure (Chirikov 1979; Casati et al. 1988; Lichtenberg & Lieberman 1992), the phase Equation in (10) can be linearized in $y$ in a vicinity of resonant phases $x = 2\pi j$ with integer $j$ describing the local dynamics by the Chirikov standard map with the chaos border $y_{cr} = (3\lambda/2K)^{2/5}$. The chaos parameter $K = K_G = 0.9716...$ corresponds to the critical golden curve (Lichtenberg & Lieberman 1992). At $K > K_G$, the dynamical chaos is global, and the chaotic diffusion from resonance to resonance becomes possible (Chirikov 1979; Lichtenberg & Lieberman 1992). However, at $K$ exceeding $K_G$ only slightly, relatively large islands of stability exist inside the global domain of chaos. At $K = 2$, the bifurcation of half-integer resonances occur. At this value, the stability islands start to disappear. The chaos border in energy is consequently

$$\Delta E_{cr} = |W_2 y_{cr}| \approx A\omega^{7/5}q^{3/10}\exp(-B\omega q^{3/2}), \quad (11)$$
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Amplitude of the energy kick $\Delta E$ as a function of the small body rotation frequency $\omega$, computed by direct integration of the dynamics of a massless particle around a rotating dumb-bell, for $q = 3d$, $q = 5d$, and $q = 7d$, and for different reduced masses $\mu = 0.01$ (□), $\mu = 0.1$ (○), $\mu = 0.4$ (Δ), and $\mu = 0.5$ (⊤). Plain lines give amplitudes of the analytically determined kick function $\Delta E$ (7) using (8) and (9). Left panel: for the sake of clarity of the figure, amplitudes of the energy kick, $\max|\Delta E|$, are presented divided by the parameter $\mu(1 - \mu)$. Inset: example of kick function $\Delta E(\phi)$ for $q = 5d$, $\mu = 0.1$, and $\omega = 0.3\omega_0$ computed from direct integration of the dynamics of a massless particle around a rotating dumb-bell (×). The green solid line gives the kick function $\Delta E(\phi)$ (7). Right panel: ratio $\max|\Delta E/W_2|$ with the same data as in the left panel.</figcaption></figure>

where $A = \mu^{2/5}(1 - \mu)^{2/5}2^{13/10}3^{2/5}\pi^{3/5}K^{-2/5}$ and $B = 2^{7/2}/15$. The half-width of the chaotic layer, $\Delta E_{cr}$, and consequently the chaos border, is qualitatively well described by this Chirikov's criterion derived formula (see Figure 6 as an illustrative example). The particle critical eccentricity $e_{cr}$, following from the relation $\Delta E_{cr} = -E_{cr} = 1/2a_{cr} = (1 - e_{cr})/2q$, is

$$e_{cr} = 1 - 2q\Delta E_{cr}, \qquad (12)$$

where $\Delta E_{cr}$ is given by (11).

Let us first consider $K = K_G$, i.e., the value from which chaos is global: orbits with $e \gtrsim e_{cr}(\omega, q)$ are chaotic. In Figure 4, the analytical curve $e_{cr}(q)$, given by (11) and (12) at $K = K_G$, is superimposed on stability diagrams for different values of $\omega$. One can see that the $e_{cr}(q)$ curve (black solid line) approximately describes the ragged border of the chaotic zone. At $K = 2$, i.e., the value at which bifurcation of half-integer resonances of the standard map occurs, the $e_{cr}(q)$ curve is shown by black dashed line in Figure 4. This curve gives the location where regular islands are no more distinguishable. The good performance of the analytical expression of $e_{cr}(q)$ for $K = K_G$ and $K = 2$ testifies the adequacy of the map's theoretical model (Popova & Shevchenko 2016).

By calculating the $e_{cr}(\omega, q)$ dependence, given by (12) at $K = K_G$, one can find the limits $q_1(\omega)$ and $q_2(\omega)$ of the central chaotic zone around the spinning irregular body; these limits ($q_1 < q_2$) are the roots of the equation $e_{cr}(q) = 0$ at $\omega$ fixed. Trajectories with $q_1 < q < q_2$ and any initial eccentricity are chaotic. In Figure 5, upper left panel, the central chaotic zone around a spinning symmetric dumb-bell ($\mu = 1/2$) is represented by the red domain. This global picture confirms that the central chaotic zone swells significantly as $\omega$ decreases. For $\mu = 1/2$, the farthest limit of the central chaotic zone, $q \simeq 7.8d$, occurs for $\omega \simeq 0.08\omega_0$. This is $\sim 2.8$ times the

farthest limit for the Keplerian frequency $\omega = \omega_0$. Conversely, the increase of $\omega$ beyond $\omega_0$ leads to a shrinking of the central chaotic zone in agreement with the stabilization effect around the fast rotating contact binary (Feng et al. 2016).

The swelling of the central chaotic zone can be explained analyzing the $\omega$ dependence of the kick amplitude $W_2$ (9) and of the width $2\Delta E_{cr}$ (11) of the chaotic layer around the separatrix ($E = 0$). Taking the example of a symmetric dumb-bell ($\mu = 1/2$), for $q = 5d$ and a spinning rate $\omega = \omega_0$, the kick amplitude, $W_2 \approx 10^{-8}(d\omega_0)^2$ (see Figure 3, left panel), is inefficient to produce chaotic orbits at any eccentricity since the lowest reachable semimajor axis is $a_{cr} = 1/(2\Delta E_{cr}) \approx 500d$ and the lowest reachable eccentricity is $e_{cr} \approx 0.99$. For $q = 5d$, but with a much slower dumb-bell spinning rate, e.g., $\omega = 0.3\omega_0$, the kick amplitude is switched on, $W_2 \approx 2 \cdot 10^{-3}(d\omega_0)^2$ (see Figure 3, left panel), in comparison to the $\omega = \omega_0$ case, giving $a_{cr} \approx q$, and thus creating a chaotic layer with orbits of any eccentricity. As a remark, we note that the swelling of the chaotic zone at $\omega < 1$ has a price: the exponent decreases, proportional to $\omega$ at $\omega/\omega_0 \ll 1$.

For $\omega \lesssim 0.24\omega_0$, a central regular zone appears in the immediate vicinity of the irregular small body. This central regular zone is surrounded by the chaotic zone and increases as $\omega$ is decreased from $\omega \simeq 0.24\omega_0$ down to $\omega \simeq 0.068\omega_0$. This central regular zone appears in a region where test particles with circular orbits have a period smaller than the rotation period of the irregular small body (see the white dotted line for 1:1 resonance in Figure 5, upper left panel). Below $\omega \simeq 0.068\omega_0$, no zeros of (12) exist and consequently no central chaotic zone exists around the irregular spinning body (i.e., for $q > 1$).

The most extended chaotic zone is provided by the symmetric case ($\mu = 1/2$, Figure 5, upper left panel). For the opposite case, at $\mu$ tending to zero, the chaotic zone vanishes, because the
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Stability diagrams for $\mu = 1/2$ and for $\omega/\omega_0 = 0.068$ (left panel), 0.4 (middle panel), and 1 (right panel). Chaotic (regular) domains are shown by reddish (blueish) areas. Chaos is determined by computing the Lyapunov exponent $\Lambda$ for a trajectory with initial orbital elements $(q, e)$. Here $10^6$ iterations of the Kepler map for dumb-bell (2) have been computed for each initial orbital elements $(q, e)$. The solid line gives the chaos border given by the analytical formula (12) with $K = K_G$. The dashed line gives the border of the bifurcation of half-integer resonances given by (12) with $K = 2$. Symbols ●, ▼, and ▲ mark the limit of the central chaotic zone (see Figure 5).</figcaption></figure>

perturbation from the second (smaller) lobe tends to zero. Hence for intermediary cases with $\mu < 1/2$, the chaotic zone is less extended, and the octopole contribution $W_1$, though weak for small $\omega$'s, is not negligible around and beyond $\omega = \omega_0$. Figure 5, upper right panel, gives the example of the central chaotic zone for non-symmetric dumb-bell with $\mu = 1/2 - \sqrt{1/12} \approx 0.211$. We have computed the analytical border of the central chaotic zone using (12) with either, as explained above, only the second harmonic contribution $W_2$ (red domain in Figure 5, upper right panel), or with only the first harmonic contribution $W_1$ instead of $W_2$ (dark red domain in Figure 5, upper right panel). We observe quite a continuous overlap between the chaotic zones induced by the two contributions $W_1$ and $W_2$, which give together a qualitative global picture of the central chaotic zone around an irregular spinning small body. The white dashed line represents the contour where $W_1 = W_2$. The dependence of $W_2$ in $\mu$ tells us that the chaotic domain induced by $W_2$ is less and less wide as $\mu$ decreases from $\mu = 1/2$ toward $\mu = 0$. The chaotic domain induced by $W_1$ is the widest for $\mu = 1/2 - \sqrt{1/12} \approx 0.211$, and is less and less wide as $\mu$ increases (decreases) from $\mu \approx 0.211$ toward $\mu = 1/2$ ($\mu = 0$). Figure 5, bottom panels, show stability diagrams of test particles initially in circular orbits ($e=0$) for the symmetric case $\mu = 1/2$ (left panel) and for the $\mu = 1/2 - \sqrt{1/12} \approx 0.211$ non-symmetric case (right panel). The fractal contour (Figure 5 bottom panels) of the central chaotic zone around the small body is well approximate by analytically obtained contours (12).

## 4. Ida and Dactyl

We now apply the Kepler map approach to real celestial bodies. Among the solar system bodies, there exists quite a marked size border line between the close-to-spherical large bodies and the essentially ellipsoidal (potato-like) small bodies. This border lies at $R = 300–500$ km, where $R$ is the characteristic radius of the body (see Figures 1–2 in Melnikov & Shevchenko 2010).

Moreover, asteroids and cometary nuclei usually resemble dumb-bells, i.e., they are more like dumb-bells than ellipsoids. A well-known example is the nucleus of comet 67P/Churyumov–Gerasimenko, the target of the Rosetta mission

(Jorda et al. 2016). Another example is asteroid 25143 Itokawa (Figure 1), the target of the Hayabusa mission (Fujiwara et al. 2006). In fact, several asteroids are observed to have a bilobed shape; in particular, 243 Ida among them, is famous to have a small natural satellite. The satellite, named Dactyl, moves in an orbit prograde with the rotation of Ida, with a very small inclination ($i < 8°$ Petit et al. 1997) with respect to the equatorial plane of Ida.

Asteroid 243 Ida can be approximately described as a symmetric dumb-bell ($\mu = 1/2$). As follows from data presented in Belton et al. (1995, 1996) and Petit et al. (1997), Ida resembles an aggregate of two merged bodies with the ratio of masses $m_2/m_1 \approx 1$ (Petit et al. 1997). We set the density $\rho$ and the rotation period $P_{rot} = 2\pi/\omega$ of the asteroid, respectively, to be equal to $2.24 \text{ g cm}^{-3}$ (Petit et al. 1997) and 4.63 hr (Vokrouhlický et al. 2003). Using formula (1), the corresponding spinning frequency for Ida is $\omega \approx 0.953\omega_0$. Besides, for the twin binary, consisting of two tangent spherical masses $m$, one has $\rho\pi d^3/3 = 2m = M$, where $M$ and $d$ are, respectively, the total mass and size of the dumb-bell. Therefore, for Ida, one has $d \approx 24.9$ km.

As an illustration (Figure 6), we show the phase portrait $(E, \phi)$ of Dactyl's dynamics around Ida obtained by iterating the Kepler map (2). As discussed above, by calculating the $e_{cr}(\omega, q)$ dependence, given by (12) at $K = K_G$, one can find the radius of the central chaotic zone around the asteroid; it is given by the root of the equation $e_{cr}(q) = 0$ at $\omega \approx 0.953\omega_0$. In the case of Ida, the root is $q \approx 2.85d \approx 71$ km. This estimate for the chaotic zone extent is in good qualitative agreement with the numerical-experimental findings on the stability limit for Dactyl's orbit size found in Petit et al. (1997).

Critical curves at $K = K_G$ and at $K = 2$ are superimposed on stability diagrams for Ida in the $(q, e)$ plane (Figure 7, left panel) and in the $(a, e)$ plane (Figure 7, middle panel); the location of Dactyl is shown by a black dot. The $51/1$ and $52/1$ resonant teeth engulf the cell where Dactyl is located. The resonances densely accumulate higher in the diagram, on approaching the parabolic separatrix. From Figure 7, it is clear that Dactyl is chaotic, in agreement with the numerical-experimental findings in Petit et al. (1997).
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Extent of the central chaotic zone around a small body as a function of the its spinning frequency ω. Upper left panel: case of a symmetric dumb-bell (μ = 1/2): the analytically obtained central chaotic zone is represented by the red domain. The blue area represents the complementary zone of stable orbits. The symbols ●, ▽, and ▲ mark the limit of the central chaotic zone for ω ≃ 0.068, 0.4, and 1, respectively (see the corresponding symbols in Figure 4). Upper right panel: case of a non-symmetric dumb-bell (here μ = 1/2 - √1/12 ≃ 0.211): the central chaotic zone in red (dark red) is obtained analytically assuming that the second (first) harmonic term in (7) with amplitude W₂ (W₁) is dominant. The white dashed line represent the curve on which W₁(q, ω) = W₂(q, ω). Bottom left and right panels: stability diagrams in the (q, ω) plane for e = 0. The reduced mass is μ = 0.5 (bottom left panel) and μ = 1/2 - √1/12 ≃ 0.211 (bottom right panel). Chaos is determined by computing the Lyapunov exponent Λ. Here 10⁶ iterations of the Kepler map for dumb-bell (2) have been computed for each couple of initial parameters (q, ω) with e = 0. Solid white lines delimit central chaotic zones obtained analytically (see upper panels). In each of the panels, white dotted lines represent p:1 and 1:p resonances. For the sake of clarity, all of the resonances, marked by dotted lines, are labeled only in the bottom left panel. The location and distribution of resonances are determined by the ratio of orbital period to dumb-bell spinning period. The resonance p': p is given by the curve q/d = (ωp/ω₀p')⁻²/³.</figcaption></figure>

Note that, in fact, short-time observations from the *Galileo*
spacecraft gave no data on the stability of Dactyl’s orbit. It can
be chaotic and thus short-lived. On the other hand, the
determination of Dactyl’s orbit may have also suffered
inaccuracies (again due to the shortness of the observations),
occasionally placing Dactyl in the chaotic region of the
diagram.

**5. Itokawa and Hayabusa**

In the case of Ida, ω/ω₀ is not far from unity, therefore, the
found central chaotic zone is analogous to the one existing
usual Keplerian binary. In our second example, 25143 Itokawa,
ω is much less than ω₀ and chaotic zone’s swelling is expected
to be large.
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Poincaré section $(E, \phi)$ for Dactyl's dynamics around Ida computed from (2) with $\mu = 1/2$. Left panel: the chaotic layer in the vicinity of the separatrix (chaotic sea) is shown by the light blue color. Ida's parameters ($d \simeq 24.9$ km, $\omega \simeq 0.953\omega_0$) have been derived from physical parameters (Petit et al. 1997; Vokrouhlický et al. 2003). A possible dynamical locus of Dactyl ($a \simeq 14.3d$, $q \simeq 3.20d$), derived from Petit et al. (1997, black point on Figure 19), is shown by the dash-dotted line. The dashed line shows the analytical estimation of the chaotic sea border according to Chirikov's criterion (11) at $K = K_G$. Top right panel: close-up around Dactyl's dynamical location. Bottom right panel: close-up of the chaotic sea border. The last invariant KAM curve separating the chaotic sea (above) from the regular domain (below) is shown in red.</figcaption></figure>

Itokawa was the target of the Hayabusa mission (Fujiwara et al. 2006). Its shape is bilobed (Figure 1), and is described as a contact binary of two ellipsoids with sizes $490 \times 310 \times 260$ m (“body”) and $230 \times 200 \times 180$ m (“head”), and densities $1750 \text{ kg m}^{-3}$ and $2850 \text{ kg m}^{-3}$, respectively; the centers of the ellipsoids are separated by $d \simeq 280$ m (Lowry et al. 2014). The period of rotation of Itokawa is 12.132 hr (Kaasalainen et al. 2003), and its mass is estimated as $3.58 \times 10^{10}$ kg (Fujiwara et al. 2006). Based on these observational data one readily calculates $P_0 = 2\pi/\omega_0 = 4.54$ h, $\omega = 0.37$, $m_1/m_2 \simeq 2.9$, and $\mu \simeq 0.26$.

The stability diagram computed on the basis of these data using the Kepler map for non-symmetric dumb-bell (2) is shown in Figure 7, right panel. The radius of the central chaotic zone $q \simeq 4.6$ is almost twice the chaotic zone we would have obtained for Itokawa’s parameters but $\omega = 1$ (not shown). We also clearly see that the central chaotic zone radius is well

estimated by the critical curves derived using only the second harmonic contribution $W_2$.

Owing to the small mass, Itokawa’s zone of gravitational influence measured by its Hill radius $R_{\text{Hill}}$ is also pretty small: it can be as small as 25 km (Fuse et al. 2008). What is more, for a probe with large solar panels such as Hayabusa, due to the effect of the solar radiation pressure, the outer limits of the zone of Itokawa’s ability to sustain satellites diminish substantially to about 3 km (Zimmer et al. 2014).

On the other hand, numerical modeling in Zimmer et al. (2014) showed “that orbits below 1 km in semimajor axis were more susceptible to the complex gravity of a rotating, non-uniform body with the spacecraft either impacting or being ejected after only a few orbits.” That is why, instead of trying to orbit Itokawa, Hayabusa moved in a neighboring orbit around the Sun. From Figure 7, right panel, it is clear that indeed no stable circular orbit can be found below $q \simeq 4.6d \simeq 1.3$ km.
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Stability diagram around Ida (left and middle panels, *d* ≃ 24.9 km, *ω* ≃ 0.953*ω*₀, *μ* ≃ 1/2) and around Itokawa (right panel, *d* ≃ 280 m, *ω* ≃ 0.37*ω*₀, *μ* ≃ 0.26): the chaotic domain is shown by the reddish area. Chaos is determined by computing the Lyapunov exponent Λ for a trajectory with initial orbital elements (*q*, *e*). Here 10⁶ iterations of the Kepler map for dumb-bell (2) have been computed for each of the initial orbital elements (*q*, *e*) (left and right panels) or (*a*, *e*) (middle panel). The black dot (left and middle panels) gives a possible current dynamical position of Dactyl according to Petit et al. (1997, Figure 19). The critical curves *e*cr(*q*) (12) for overlap of integer resonances *K* = *K*G (solid line) and for bifurcation of half-integer resonances *K* = 2 (dashed line) are plotted.</figcaption></figure>

Itokawa has no satellites, as reported in Fuse et al. (2008). The formation of the extended central chaotic zone, in concert with the smallness of the Hill sphere, explains the lack of moons. This effect also explains why Hayabusa could not be put in orbit around Itokawa.

**6. Capture Cross-section**

Particles flying by a non-spherical spinning body can be captured. Following Lages & Shepelyansky (2013), the capture cross-section σ characterizes the probability that a spinning body captures a scattering particle after a passage at the pericenter. The fact that chaotic zones increases significantly at ω/ω₀ ≪ 1 leads to an increase of the capture cross-section σ. Indeed, according to Lages & Shepelyansky (2013), we have σ ~ πr<sub>st</sub>² ~ qd(dω₀/v<sub>st</sub>)², where r<sub>st</sub> and v<sub>st</sub> are the impact distance and the mean velocity of a scattering particle at infinity. Since from (8)–(9) the exchange of energy (7) is non-negligible for pericenters up to q ~ d(ω₀/ω)²/³ the above estimate shows that the capture cross-section of the slowly spinning body (ω/ω₀ ≪ 1) can be significantly enhanced compared to its geometric cross-section ~πd². Such an effect may play an important role for dust capture by, e.g., a spinning satellite.

**7. Conclusions**

We have generalized the Kepler map technique to describe the motion of a particle in the gravitational field of a rotating irregular body modeled by a dumb-bell. This has been achieved by the introduction of an additional parameter responsible for the arbitrary rate of rotation of the “central binary.” We have found that the chaotic zone around the dumb-bell swells significantly if its rotation rate is decreased; in particular, the zone swells more than twice if the rotation rate is decreased 10 times with respect to the “centrifugal breakup” threshold. We have determined the extent of the chaotic zone both analytically and numerically.

To connect our theoretical findings with observational data, we have illustrated the properties of the chaotic orbital zones in examples of the global orbital dynamics about asteroid 243 Ida

(which has a moon, Dactyl, orbiting near the edge of the chaotic zone) and asteroid 25143 Itokawa.

Possible orbital regimes of Ida's moon Dactyl have been described by means of constructing stability diagrams of its orbital motion. The qualitative dynamics of the Ida-Dactyl asteroid-satellite system has been shown to be described adequately within this approach; in particular, an agreement has been found with previous numerical-experimental data on the stability of orbits around Ida. It has been explained why Dactyl is marginally chaotic, as its orbit is situated at the fractal border of the analytically expected central chaotic zone.

For Itokawa, it has been explained why space probe Hayabusa could not be put in orbit around it, and also why Itokawa has no natural satellites. All this is due to the swelling of the chaotic zone around slowly rotating Itokawa, this enlargement being combined with the smallness of its Hill sphere.

We highlight various important implications of emerged chaos around rotating minor bodies. The dumb-bell map technique might be perspectively applied to describe orbital motions and assess the possibility and sizes of chaotic zones around elongated minor planetary satellites, e.g., minor moons in the Pluto-Charon system. Indeed, as outlined in Quillen et al. (2017), in this system, only Hydra rotates rapidly (at the rate of ~30% of the “centrifugal breakup” threshold). Therefore, the chaotic zones around the minor moons in the Pluto-Charon system may engulf their Hill spheres substantially; this issue deserves further study.

The authors are thankful to Jean-Marc Petit, Darin Ragozzine, and the anonymous referee for valuable remarks and comments. I.I.S. benefited from a grant of Bourgogne-Franche-Comté region. I.I.S. was supported in part by the Russian Foundation for Basic Research (project No. 17-02-00028).

**References**

Bartczak, P., & Breiter, S. 2003, CeMDA, 86, 131

Batygin, K., & Morbidelli, A. 2015, ApJ, 810, 110

Belton, M. J., Chapman, C. R., Klaasen, K. P., et al. 1996, Icar, 120, 1

Belton, M. J. S., Chapman, C. R., Thomas, P. C., et al. 1995, Natur, 374, 785

Benvenuto, F., Casati, G., & Shepelyansky, D. L. 1994, PhRvL, 72, 1818
Casati, G., Guarneri, I., & Shepeliansky, D. L. 1988, IQE, 24, 1420

Chauvineau, B., Farinella, P., & Mignard, F. 1993, Icar, 105, 350

Chirikov, B. V. 1979, PhR, 52, 263

Chirikov, B. V., & Vecheslavov, V. V. 1989, A&A, 221, 146

Emelyanenko, V. V. 1990, SvAL, 16, 318

Feng, J., Noomen, R., Hou, X., Visser, P., & Yuan, J. 2017, CeMDA, 127, 67

Feng, J., Noomen, R., Visser, P., & Yuan, J. 2016, AdSpR, 58, 387

Fujiwara, A., Kawaguchi, J., Yeomans, D. K., et al. 2006, Sci, 312, 1330

Fuse, T., Yoshida, F., Tholen, D., Ishiguro, M., & Saito, J. 2008, EP&S, 60, 33

Gaskell, R., Saito, J., Ishiguro, M., et al. 2008, HAY-A-AMICA-5-ITOKAWASHAPE-V1.0., NASA Planetary Data System

Heggie, D. C. 1975, MNRAS, 173, 729

Hergenrother, C. W., & Whiteley, R. J. 2011, Icar, 214, 194

Hu, W., & Scheeres, D. J. 2004, P&SS, 52, 685

Hut, P. 1981, A&A, 99, 126

Jorda, L., Gaskell, R., Capanna, C., et al. 2016, Icar, 277, 257

Kaasalainen, M., Kwiatkowski, T., Abe, M., et al. 2003, A&A, 405, L29

Lages, J., & Shepelyansky, D. L. 2013, MNRAS Letters, 430, L25

Lichtenberg, A. J., & Lieberman, M. A. 1992, Regular and Chaotic Dynamics (New York: Springer)

Lowry, S. C., Weissman, P. R., Duddy, S. R., et al. 2014, A&A, 562, A48

Malyshkin, L., & Tremaine, S. 1999, Icar, 141, 341

Marchis, F., Durech, J., Castillo-Rogez, J., et al. 2014, ApJL, 783, L37

Meiss, J. D. 1992, RvMP, 64, 795

Melnikov, A., & Shevchenko, I. 2010, Icar, 209, 786

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

Mysen, E., Olsen, Ø., & Aksnes, K. 2006, P&SS, 54, 750

Olsen, Ø. 2006, A&A, 449, 821

Petit, J., Durda, D., Greenberg, R., Hurford, T., & Geissler, P. 1997, Icar, 130, 177

Petrosky, T. Y. 1986, PhLA, 117, 328

Popova, E. A., & Shevchenko, I. I. 2016, AstL, 42, 474

Pravec, P., Harris, A., Vokrouhlický, D., et al. 2008, Icar, 197, 497

Quillen, A. C., Nichols-Fleming, F., Chen, Y.-Y., & Noyelles, B. 2017, arXiv:1701.05594

Rollin, G., Haag, P., & Lages, J. 2015, PhLA, 379, 1017

Rollin, G., Lages, J., & Shepelyansky, D. L. 2015, A&A, 576, A40

Roy, A., & Haddow, M. 2003, CeMDA, 87, 411

Scheeres, D. 2007, Icar, 189, 370

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

Scheeres, D. J. 2012, AcAau, 72, 1

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

Scheeres, D. J., Williams, B. G., & Miller, J. K. 2000, JGCD, 23, 466

Shevchenko, I. I. 2010, PhRvE, 81, 066216

Shevchenko, I. I. 2011, NewA, 16, 94

Shevchenko, I. I. 2015, ApJ, 799, 8

Vokrouhlický, D., Nesvorný, D., & Bottke, W. F. 2003, Natur, 425, 147

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

Werner, R. A. 1994, CeMDA, 59, 253

Werner, R. A., & Scheeres, D. J. 1996, CeMDA, 65, 313

Whiteley, R. J., Hergenrother, C. W., & Tholen, D. J. 2002, in ESA Special Publication 500, Asteroids, Comets, and Meteors: ACM, ed. B. Warmbein (Noordwijk: ESA), 473

Yu, Y., & Baoyin, H. 2012, AJ, 143, 62

Zimmer, M., Williams, K., Johnson, J., et al. 2014, LPSC, 45, 2226