Chaotic zones around rotating small bodies

TL;DR

This paper analyzes the extent of chaotic orbital zones around rotating small bodies in the Solar system, showing how these zones expand significantly with slower rotation rates, with implications for asteroid dynamics.

Contribution

It provides an analytical and numerical determination of chaotic zone extents around irregularly shaped, rotating small bodies, highlighting the impact of rotation rate on zone size.

Findings

Chaotic zones expand more than twice when rotation rate decreases tenfold.

The extent of chaotic zones is analytically and numerically characterized.

Examples include asteroid 243 Ida and asteroid 25143 Itokawa.

Abstract

Small bodies of the Solar system, like asteroids, trans-Neptunian objects, cometary nuclei, planetary satellites, with diameters smaller than one thousand kilometers 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 ten 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.

Equations

$${\omega }_0=\sqrt{\pi \mathcal{G}\rho /3},$$

$$\begin{array}[]{lll}E_{i+1}=E_{i}+\Delta E\left(\phi_{i}\right),\quad\phi_{i+1}&=&\phi_{i}+2\pi\omega/|2E_{i+1}|^{3/2}\end{array}$$

$$\ddot{\mathbf{R}}=-{\nabla }_{\mathbf{R}}\Phi \left(\mathbf{R},\mathbf{r},\mu ,\omega \right)$$

$$\Phi \left(\mathbf{R},\mathbf{r},\mu ,\omega \right)=-\frac{1-\mu }{\Vert \mathbf{R}+\mu \mathbf{r}\Vert }-\frac{\mu }{\Vert \mathbf{R}-(1-\mu )\mathbf{r}\Vert }.$$

$$\Phi \left(\mathbf{R},\mathbf{r},\mu ,\omega \right)=-\frac{1}{R}-\mu \left(1-\mu \right)\frac{r^2}{2R^3}\left(3{\left(\frac{\mathbf{r}\cdot \mathbf{R}}{rR}\right)}^2-1\right)-\mu \left(1-\mu \right)\left(2\mu -1\right)\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)$$

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

$$\Delta E\left(\mu ,q,\omega ,\varphi \right)\simeq W_1\sin \left(\varphi \right)+W_2\sin \left(2\varphi \right).$$

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

$$W_2\simeq -\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).$$

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

$$\Delta E_{\mathrm{cr}}=\left|W_2{y}_{\mathrm{cr}}\right|\approx A{\omega }^{7/5}{q}^{3/10}\exp \left(-B\omega q^{3/2}\right),$$

$$e_{\mathrm{cr}}=1-2q\Delta E_{\mathrm{cr}},$$

Full text

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, planetary satellites) with diameters less than one thousand kilometers usually have strongly irregular shapes (Melnikov & Shevchenko, 2010; Jorda et al., 2016, p. 270), 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), 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; Marchis et al., 2014; Yu & Baoyin, 2012). Concerning the dumb-bell model, it was also used in the problem on 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 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 particle’s orbital motion and the rotational motion of the central elongated body was revealed by Mysen et al. (2006); Olsen (2006); 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 loose 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 analysis of symplectic maps (see eg 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 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, so that 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, Figs. 4,5,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). But 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 (Chirikov & Vecheslavov, 1989; Petrosky, 1986) 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 it is shown in Fig. 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 Fig. 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 Fig. 1). The Keplerian rate of rotation of a contact binary, i.e. two tangent spheres, is

[TABLE]

where $\rho$ is the density of the irregular body (Scheeres, 2007). For a typical density $\rho=1$g$/$cm3 we have $\omega_{0}=2.5\times 10^{-4}$ s*-1* corresponding to a period of about $7$ hours. There are many observed asteroids with significantly larger rotation periods (see eg 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}\leq m_{1}$ and we define $\mu=m_{2}/M\leq 0.5$), $d=1$ is the size of the effective dumb-bell (Fig. 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 $\omega$ of the dumb-bell, takes the form (Casati et al., 1988; Benvenuto et al., 1994)

[TABLE]

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, ie $\omega/|2E|^{3/2}\leq 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\neq\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\left(\omega t\right)\mathbf{\hat{a}}+\sin\left(\omega t\right)\mathbf{\hat{b}}$ the position of the dumb-bell lobe of mass $m_{2}$ relative to the dumb-bell lobe of mass $m_{1}$ ($\mathbf{\hat{a}}$ and $\mathbf{\hat{b}}$ are any two orthogonal fixed directions of the plane), and $\mathbf{R}$ 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

[TABLE]

where the gravitational potential reads

[TABLE]

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

[TABLE]

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

[TABLE]

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 arbitrary frequency rate $\omega$, keeping the two first leading terms for the kick function (6) we obtain

[TABLE]

In equation (7), the exchange of energy between the small spinning body and the test particle is splat 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

[TABLE]

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

[TABLE]

We note that expression (8) restricted to the case $\omega=\omega_{0}=1$ and $\mu\ll 1$ has been obtained using 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 consider the Kepler map description of cometary dynamics around the Solar System modelized by the Sun and Jupiter as perturber (Chirikov & Vecheslavov, 1989; Petrosky, 1986). But 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}$h*-1* to $\omega/2\pi\sim 10^{2}$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 $2$ 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$ (Fig. 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 octopole 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 Fig. 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 (Fig. 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 ordinary binary, such an energy kick increase induced by a slow spinning frequency allows zone of chaos to extend quite far from the central body. In Fig. 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 Fig. 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 size 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$ (Fig. 4). The border delimiting chaotic domain (red) from regular domain (blue) is ragged; the most prominent teeth being associated to the integer $p$:1 and half-integer $p+\frac{1}{2}$:1 resonances between particle orbital frequency and small body spinning frequency. Here any neighboring integer and half-integer resonances are equal-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 stated the chaotic zone is defined as the region where even particles initially in circular orbits become dynamically chaotic. From Fig. 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 Eqs. (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 is greater than $10$ for $\mu=0$ and diverges as $\mu$ approaches $1/2$. By the substitution $E=W_{2}\,y$ and $\phi=x/2$ the map (2) is reduced to

[TABLE]

with $\lambda=2^{1/2}\pi\omega/|W_{2}|^{3/2}$. Following the standard procedure (Chirikov, 1979; Lichtenberg & Lieberman, 1992; Casati et al., 1988) 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_{\mathrm{G}}=0.9716\dots$ corresponds to the critical golden curve (Lichtenberg & Lieberman, 1992). At $K>K_{\mathrm{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_{\mathrm{G}}$ only slightly, relatively large islands of stability exist inside the global domain of chaos. At $K=2$ bifurcation of half-integer resonances occur. At this value the stability islands start to disappear. The chaos border in energy is consequently

[TABLE]

where $A=\mu^{2/5}\left(1-\mu\right)^{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_{\mathrm{cr}}$, and consequently the chaos border, is qualitatively well described by this Chirikov’s criterion derived formula (see Fig. 6 as an illustrative example). The particle critical eccentricity $e_{\mathrm{cr}}$, following from the relation $\Delta E_{\mathrm{cr}}=-E_{\mathrm{cr}}=1/2a_{\mathrm{cr}}=(1-e_{\mathrm{cr}})/2q$, is

[TABLE]

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

Let us first consider $K=K_{\mathrm{G}}$, i.e. the value from which chaos is global: orbits with $e\gtrsim e_{\mathrm{cr}}(\omega,q)$ are chaotic. In Fig. 4, the analytical curve $e_{\mathrm{cr}}(q)$, given by (11) and (12) at $K=K_{\mathrm{G}}$, is superimposed on stability diagrams for different values of $\omega$. 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 Fig. 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_{\mathrm{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_{\mathrm{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_{\mathrm{cr}}(q)=0$ at $\omega$ fixed. Trajectories with $q_{1}<q<q_{2}$ and any initial eccentricity are chaotic. In Fig. 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 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_{\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}\left(d\omega_{0}\right)^{2}$ (see Fig. 3, left panel), is inefficient to produce chaotic orbits at any eccentricity since the lowest reachable semi-major axis is $a_{\mathrm{cr}}=1/(2\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}\left(d\omega_{0}\right)^{2}$ (see Fig. 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<1$ has some price: the Lyapunov exponent decreases being 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 white dotted line for 1:1 resonance in Fig. 5 upper left panel). Below $\omega\simeq 0.068\omega_{0}$, no zeros of (12) exists 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$, Fig. 5, upper left panel). For the opposite case, at $\mu$ tending to zero the chaotic zone vanishes, because the 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}$. Fig. 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_{2}$ (red domain in Fig. 5, upper right panel), or with only the first harmonic contribution $W_{1}$ instead of $W_{2}$ (dark red domain in Fig. 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 represent 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}\simeq 0.211$, and is less and less wide as $\mu$ increases (decreases) from $\mu\simeq 0.211$ toward $\mu=1/2$ ($\mu=0$). Fig. 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 (Fig. 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, usually asteroids and cometary nuclei 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 (Fig. 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^{\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); 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 $\rho$ and the rotation period $P_{\mathrm{rot}}=2\pi/\omega$ of the asteroid, respectively, to be equal to $2.24$ g.cm*-3* (Petit et al., 1997) and $4.63$ h (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 (Fig. 6) we show the phase portrait $(E,\phi)$ of Dactyl’s dynamics around Ida obtained by iteration of the Kepler map (2). As discussed above, by calculating the $e_{\mathrm{cr}}(\omega,q)$ dependence, given by (12) at $K=K_{\mathrm{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).

Critical curves at $K=K_{\mathrm{G}}$ and at $K=2$ are superimposed on stability diagrams for Ida in the $(q,e)$ plane (Fig. 7, left panel) and in the $(a,e)$ plane (Fig. 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 Figs. 7 it is clear that Dactyl is chaotic, in agreement with the numerical-experimental findings in Petit et al. (1997).

Note that, in fact, short-time observations from the Galileo spacecraft gave no data on the stability of Dactyl’s orbit. It can well 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, $\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, $\omega$ 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 (Fig. 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 kg/m3 and 2850 kg/m3, respectively; the centers of the ellipsoids are separated by $d\simeq 280$ m (Lowry, S. C. et al., 2014). The period of rotation of Itokawa is 12.132 h (Kaasalainen, M. 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$, $\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 Fig. 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_{\mathrm{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 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 Fig. 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.

6 Capture cross-section

Particles flying by a non spherical spinning body can be captured. Following Lages & Shepelyansky (2013), the capture cross-section $\sigma$ 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 $\sigma$. Indeed, according to Lages & Shepelyansky (2013) we have $\sigma\sim\pi r_{st}^{2}\sim qd(d\omega_{0}/v_{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 slowly spinning body $(\omega/\omega_{0}\ll 1)$ can be significantly enhanced comparing 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.

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 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 ten 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, eg, 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 $\sim$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 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).