Chaotic Zones around Rotating Small Bodies

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 straightforwardly 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 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 loose 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 ω 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 ω, 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.

Zoom In Zoom Out Reset image size

The Kepler map description of orbits about a spinning non-axisymmetric body is achieved by introducing a parameter, ω, which is the rate of rotation of the model contact binary (see Figure 1). The value of ω 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 \lt {\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 \lt {\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.

We consider the motion of a passively gravitating particle in the planar circular restricted three-body problem m₁–m₂-particle, where the two masses m₁ and m₂ 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

$$\begin{eqnarray}&&{\omega }_{0}=\sqrt{\pi { \mathcal G }\rho /3},\end{eqnarray} \tag{ 1 }$$

where ρ is the density of the irregular body (Scheeres 2007). For a typical density $\rho =1$ g/cm³, we have ${\omega }_{0}=2.5\times {10}^{-4}$ s⁻¹ 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: ${ \mathcal G }M=1$ (where $M={m}_{1}+{m}_{2}$ is the total mass of the irregular gravitating body, we choose ${m}_{2}\leqslant {m}_{1}$ and we define $\mu ={m}_{2}/M\leqslant 0.5$), d = 1 is the size of the effective dumb-bell (Figure 1), and the Keplerian frequency ${\omega }_{0}=\sqrt{{ \mathcal 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 ω of the dumb-bell, takes the form (Casati et al. 1988; Benvenuto et al. 1994)

$$\begin{eqnarray}{E}_{i+1}={E}_{i}+{\rm{\Delta }}E({\phi }_{i}),\quad {\phi }_{i+1} & = & {\phi }_{i}+2\pi \omega /| 2{E}_{i+1}{| }^{3/2}\end{eqnarray} \tag{ 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}\geqslant 1$ (Casati et al. 1988; Benvenuto et al. 1994). We obtain the kick function ${\rm{\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 ${\boldsymbol{r}}=\cos (\omega t)\hat{{\boldsymbol{a}}}+\sin (\omega t)\hat{{\boldsymbol{b}}}$ is the position of the dumb-bell lobe of mass m₂ relative to the dumb-bell lobe of mass m₁ ($\hat{{\boldsymbol{a}}}$ and $\hat{{\boldsymbol{b}}}$ are any two orthogonal fixed directions of the plane), and ${\boldsymbol{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

$$\begin{eqnarray}&&\ddot{{\boldsymbol{R}}}=-{{\boldsymbol{\nabla }}}_{{\boldsymbol{R}}}{\rm{\Phi }}({\boldsymbol{R}},{\boldsymbol{r}},\mu ,\omega ),\end{eqnarray} \tag{ 3 }$$

where the gravitational potential reads

$$\begin{eqnarray}&&{\rm{\Phi }}({\boldsymbol{R}},{\boldsymbol{r}},\mu ,\omega )=-\displaystyle \frac{1-\mu }{\parallel {\boldsymbol{R}}+\mu {\boldsymbol{r}}\parallel }-\displaystyle \frac{\mu }{\parallel {\boldsymbol{R}}-(1-\mu ){\boldsymbol{r}}\parallel }.\end{eqnarray} \tag{ 4 }$$

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

$$\begin{eqnarray}{\rm{\Phi }}({\boldsymbol{R}},{\boldsymbol{r}},\mu ,\omega ) & = & -\,\displaystyle \frac{1}{R}-\mu (1-\mu )\displaystyle \frac{{r}^{2}}{2{R}^{3}}\\ & & \times \left(3{\left(\displaystyle \frac{{\boldsymbol{r}}\cdot {\boldsymbol{R}}}{{rR}}\right)}^{2}-1\right)-\mu (1-\mu )(2\mu -1)\\ & & \times \displaystyle \frac{{r}^{3}}{2{R}^{4}}\left(5{\left(\displaystyle \frac{{\boldsymbol{r}}\cdot {\boldsymbol{R}}}{{rR}}\right)}^{3}-3\displaystyle \frac{{\boldsymbol{r}}\cdot {\boldsymbol{R}}}{{rR}}\right)+O\left(\displaystyle \frac{{r}^{4}}{{R}^{5}}\right).\end{eqnarray} \tag{ 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

$$\begin{eqnarray}&&{\rm{\Delta }}E(\mu ,q,\omega ,\phi )=-{\int }_{-\infty }^{+\infty }\dot{{\boldsymbol{R}}}\cdot {\boldsymbol{\nabla }}\left({\rm{\Phi }}+\displaystyle \frac{1}{R}\right){dt}\end{eqnarray} \tag{ 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 ϕ 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 ω, keeping the two first leading terms for the kick function (6) we obtain

$$\begin{eqnarray}&&{\rm{\Delta }}E(\mu ,q,\omega ,\phi )\simeq {W}_{1}\sin (\phi )+{W}_{2}\sin (2\phi ).\end{eqnarray} \tag{ 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

$$\begin{eqnarray}{W}_{1} & \simeq & \mu (1-\mu )(1-2\mu ){2}^{1/4}{\pi }^{1/2}{\omega }^{5/2}{q}^{-1/4}\\ & & \times \exp \,\left(-\displaystyle \frac{{2}^{3/2}}{3}\omega {q}^{3/2}\right)\end{eqnarray} \tag{ 8 }$$

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

$$\begin{eqnarray}&&{W}_{2}\simeq -\mu (1-\mu ){2}^{15/4}{\pi }^{1/2}{\omega }^{5/2}{q}^{3/4}\exp \left(-\displaystyle \frac{{2}^{5/2}}{3}\omega {q}^{3/2}\right).\end{eqnarray} \tag{ 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 ${\rm{\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 μ, 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}$ h⁻¹ to $\omega /2\pi \sim {10}^{2}\,{{\rm{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₁ 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₁ over W₂ is even quadratically shifted farther from the small body for $\omega \lt {\omega }_{0}$. Conversely, which contribution, either W₁ or W₂, dominates is not so obvious for the region $q\lesssim {({\omega }_{0}/\omega )}^{2/3}d$, which for $\omega \lt {\omega }_{0}$ encompasses the immediate vicinity of the spinning small body. The two contributions W₁ (8) and W₂ (9) depend on the parameters μ, ω, and q; their relative importance is summarized in the $(\omega ,q)$ plot for different values of μ (Figure 2). We clearly see that below the frequency of disruption of a rubble-pile object ($\omega \lt {\omega }_{0}$), for any mass parameter μ, the quadrupole coefficient W₂ generally dominates the octopole coefficient W₁ 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 μ parameter.

Zoom In Zoom Out Reset image size

Typical amplitudes of energy kick functions ${\rm{\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 ${\rm{\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 (${\rm{\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 (${\rm{\Delta }}E\sim {10}^{-2}{d}^{2}{\omega }_{0}^{2}$). In comparison with an ordinary binary, such an energy kick increase induced by 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 ${\rm{\Delta }}E$ are presented divided by the mass factor $\mu (1-\mu )$ entering the expression of W₂ (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 μ are superimposed stressing again the fact that the second harmonic term is dominant for small spinning frequencies (see also Figure 3, right panel).

Zoom In Zoom Out Reset image size

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 μ 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.

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+\tfrac{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}$.

Zoom In Zoom Out Reset image size

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₂ clearly dominates over W₁ for small spinning frequencies ($\omega \lt {\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 μ approaches 1/2. By the substitution $E={W}_{2}\,y$ and $\phi =x/2$, the map (2) is reduced to

$$\begin{eqnarray}&&{y}_{i+1}={y}_{i}+\sin {x}_{i},\qquad {x}_{i+1}={x}_{i}+\lambda /| {y}_{i+1}{| }^{3/2}\end{eqnarray} \tag{ 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}_{\mathrm{cr}}={(3\lambda /2K)}^{2/5}$. The chaos parameter $K\,={K}_{{\rm{G}}}=0.9716...$ corresponds to the critical golden curve (Lichtenberg & Lieberman 1992). At $K\gt {K}_{{\rm{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

$$\begin{eqnarray}&&{\rm{\Delta }}{E}_{\mathrm{cr}}=| {W}_{2}\,{y}_{\mathrm{cr}}| \approx A{\omega }^{7/5}{q}^{3/10}\exp (-B\omega \,{q}^{3/2}),\end{eqnarray} \tag{ 11 }$$

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, ${\rm{\Delta }}{E}_{\mathrm{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 ${\rm{\Delta }}{E}_{\mathrm{cr}}=-{E}_{\mathrm{cr}}=1/2{a}_{\mathrm{cr}}=(1-{e}_{\mathrm{cr}})/2q$, is

$$\begin{eqnarray}&&{e}_{\mathrm{cr}}=1-2q{\rm{\Delta }}{E}_{\mathrm{cr}},\end{eqnarray} \tag{ 12 }$$

where ${\rm{\Delta }}{E}_{\mathrm{cr}}$ is given by (11).

Let us first consider $K={K}_{{\rm{G}}}$, i.e., the value from which chaos is global: orbits with $e\gtrsim {e}_{\mathrm{cr}}(\omega ,q)$ are chaotic. In Figure 4, the analytical curve ${e}_{\mathrm{cr}}(q)$, given by (11) and (12) at $K={K}_{{\rm{G}}}$, is superimposed on stability diagrams for different values of ω. One can see that the ${e}_{\mathrm{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}_{\mathrm{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}_{\mathrm{cr}}(q)$ for $K={K}_{{\rm{G}}}$ and K = 2 testifies the adequacy of the map’s theoretical model (Popova & Shevchenko 2016).

By calculating the ${e}_{\mathrm{cr}}(\omega ,q)$ dependence, given by (12) at $K={K}_{{\rm{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}\lt {q}_{2}$) are the roots of the equation ${e}_{\mathrm{cr}}(q)=0$ at ω fixed. Trajectories with ${q}_{1}\lt q\lt {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 ω 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 ∼2.8 times the farthest limit for the Keplerian frequency $\omega ={\omega }_{0}$. Conversely, the increase of ω 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).

Zoom In Zoom Out Reset image size

The swelling of the central chaotic zone can be explained analyzing the ω dependence of the kick amplitude W₂ (9) and of the width $2{\rm{\Delta }}{E}_{\mathrm{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}_{\mathrm{cr}}=1/(2{\rm{\Delta }}{E}_{\mathrm{cr}})\approx 500d$ and the lowest reachable eccentricity is ${e}_{\mathrm{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}_{\mathrm{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 \lt 1$ has a price: the exponent decreases, proportional to ω 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 ω 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\gt 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 μ tending to zero, the chaotic zone vanishes, because the perturbation from the second (smaller) lobe tends to zero. Hence for intermediary cases with $\mu \lt 1/2$, the chaotic zone is less extended, and the octopole contribution W₁, though weak for small ω'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}\simeq 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₂ (red domain in Figure 5, upper right panel), or with only the first harmonic contribution W₁ instead of W₂ (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₁ and W₂, 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₂ in μ tells us that the chaotic domain induced by W₂ is less and less wide as μ decreases from $\mu =1/2$ toward $\mu =0$. The chaotic domain induced by W₁ is the widest for $\mu =1/2-\sqrt{1/12}\,\simeq 0.211$, and is less and less wide as μ increases (decreases) from $\mu \simeq 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}\simeq 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).

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\lt 8^\circ$ 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}\simeq 1$ (Petit et al. 1997). We set the density ρ and the rotation period ${P}_{\mathrm{rot}}=2\pi /\omega$ of the asteroid, respectively, to be equal to 2.24 g cm⁻³ (Petit et al. 1997) and 4.63 hr (Vokrouhlický et al. 2003). Using formula (1), the corresponding spinning frequency for Ida is $\omega \simeq 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\simeq 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}_{\mathrm{cr}}(\omega ,q)$ dependence, given by (12) at $K={K}_{{\rm{G}}}$, one can find the radius of the central chaotic zone around the asteroid; it is given by the root of the equation ${e}_{\mathrm{cr}}(q)=0$ at $\omega \simeq 0.953{\omega }_{0}$. In the case of Ida, the root is $q\simeq 2.85d\simeq 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).

Zoom In Zoom Out Reset image size

Critical curves at $K={K}_{{\rm{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).

Zoom In Zoom Out Reset image size

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.

In the case of Ida, $\omega /{\omega }_{0}$ 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 ${\omega }_{0}$ and chaotic zone’s swelling is expected to be large.

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 × 310 × 260 m (“body”) and 230 × 200 × 180 m (“head”), and densities 1750 kg m⁻³ and 2850 kg m⁻³, 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\,{\rm{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₂.

Owing to the small mass, Itokawa’s zone of gravitational influence measured by its Hill radius R_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.

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.

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 $\omega /{\omega }_{0}\ll 1$ leads to an increase of the capture cross-section σ. Indeed, according to Lages & Shepelyansky (2013), we have $\sigma \sim \pi {r}_{\mathrm{st}}^{2}\,\sim {qd}{(d{\omega }_{0}/{v}_{\mathrm{st}})}^{2}$, where r_st and v_st 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\sim d{({\omega }_{0}/\omega )}^{2/3}$ the above estimate shows that the capture cross-section of the slowly spinning body $(\omega /{\omega }_{0}\ll 1)$ can be significantly enhanced compared to its geometric cross-section $\sim \pi {d}^{2}$. Such an effect may play an important role for dust capture by, e.g., a spinning satellite.

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).