# Kepler map

Post-publication activity

Curator: Dima Shepelyansky

The Kepler map is a two-dimensional symplectic map

$\tag{1} \bar{w} = w + J \sin \phi = w +F(\phi/2\pi)\; , \; \bar{\phi} = \phi + 2\pi/{\bar{w}}^{3/2} \;.$ Figure 1: Poincaré sections given by the generalized Kepler map with several kick harmonics for the case of Halley comet (top) and the Kepler map (1) with corresponding $$J =0.007$$ at sine kick (bottom); white area corresponds to invariant curves (from Rollin et al., 2015b; see details there). Figure 2: Comet Halley orbit: an orthographic projection (left panel) and an arbitrary view point (right panel). The red curve shows three successive orbital revolutions of comet Halley, the other near circular orbits are for the eight planets, the yellow bright spot gives the Sun position. (from Rollin et al., 2015a) Figure 3: A schematic image of orbit of third light body (3) in the plane of rotating binary with masses $$M$$ (1) and $$m_p$$ (2) with corresponding distance at perihelion $$q$$ (closest distance between light body and center of mass of binary), large semi-axis $$a$$, effective planet radius $$a_p$$ and rescaled energy $$w$$, arrows show the rotation direction.

which describes the highly-eccentric circumbinary motion of a massless (passively gravitating) particle in the restricted three-body problem; it also describes the microwave ionization of hydrogen and Rydberg atoms with initial highly excited states in the semiclassical regime. Here $$w = -2E/(m_d {v_p}^2)$$ is a rescaling of the energy $$E$$ of the particle with mass $$m_d$$, $$v_p$$ is the velocity of a planet of mass $$m_p$$ orbiting around a star with mass $$M > m_p$$; the phase $$x=\phi/2\pi=t/T_p \pmod 1$$ is given by time $$t$$ taken at the moment of the particle's passage through perihelion, the planet's orbital period and radius are $$T_p$$ and $$a_p$$; variables with overbars $$(\bar{w},\bar{x})$$ mark the new values of variables $$(w,x)$$ after one orbital period of the particle given by the third Kepler law (Kepler, 1619). A kick in energy $$\bar{w}-w$$ takes place during the passage through the perihelion, which has a minimal distance between a particle and the center of mass of binary. For the perihelion distance $$q\gg a_p$$, the kick is a sine function with an amplitude $$J$$ which is proportional to the mass ratio $$m_p/M$$ and depends also on the perihelion distance $$q$$. For $$q \sim a_p$$ one has a generalized Kepler map with a kick function $$J\sin(2\pi x) \rightarrow F(x)$$ containing several harmonics. The generalized Kepler map describes, e.g., dynamics of long-period and Halley-type comets, including comet Halley, and dark matter in the Solar System. The Kepler map applications in atomic physics of excited states are described at Microwave ionization of hydrogen atoms. Typical examples of the Poincaré sections given by the Kepler and generalized Kepler maps are shown in Fig.1 (see also Figs below), the orbit of comet Halley is presented in Fig.2, a schematic image of orbit is shown in Fig.3 with all main geometric notations.

# Map derivation and properties

The change of the particle energy, $$\Delta E \propto\bar{w}-w$$, is obtained by integration of energy variation, produced by the planet's orbital motion around the Sun, over an unperturbed Kepler orbit with eccentricity close to unity (in the limit over a parabolic orbit). Such integration is similar to the Melnikov integral appearing in the computation of the separatrix splitting (see Chirikov, 1979). This integration in the planar circular restricted three-body problem gives the following amplitude for the particle energy variation $\tag{2} J \approx \mu (1-\mu) (1-2\mu) 2^{5/4} \pi^{1/2} \omega^{5/2} q^{-1/4} \exp(-2^{3/2} \omega q^{3/2}/3) \; ,$ where $$\mu = m_p/(M+m_p)$$, $$q \approx \ell^2/2$$ is the perihelion distance ($$\ell$$ is the particle orbital momentum per unit of mass) and the variables are expressed in units with $$v_p=a_p=1$$ (see Petrosky, 1986,Petrosky and Broucke, 1988,Shevchenko, 2011,Lages et al., 2017). The frequency $$\omega$$ notes here the orbital frequency of the planet expressed in units of its Kepler orbital frequency $$\omega_0 = v_p/a_p$$. In the three-body problem, in whose framework the original Kepler map was derived, $$\omega=\omega_0=1$$, but for the problem of microwave ionization of hydrogen atoms (see also Benvenuto et al., 1994,Casati et al., 1987,Casati et al., 1988) or rotating rigid bodies (see below) these two frequencies can be different. The expression for $$J$$ is given for the orbit of the particle moving in a prograde direction (the same direction as the planet), for the case of retrograde particle's orbit the expression is slightly different (Petrosky and Broucke, 1988). In Eq.(2) it is also assumed that $$q \gg 1$$, thus there are no close encounters between the particle and the planet. In this regime $$J$$ is independent of $$w$$. Here, the planet produces a high frequency perturbation on the particle motion and due to analyticity $$J$$ decays exponentially with $$q$$.

The kick contribution to the energy change in Eq.(1) happens during the passage of the particle through perihelion with a time duration being much smaller than the particle period given by the third Kepler law (Kepler, 1619, see also WikiKepler, 2017) for the phase change in Eq.(1). The sinusoidal dependence of energy change on the planet phase corresponds to the main time dependent dipole contribution of gravitational interactions between the particle and the orbiting planet. The dipole term is dominant since $$q > a_p=1$$. In fact the problem of a central static charge (analogous to Sun) and rotating dipole (analogous to orbiting planet), interacting with an electron moving over a Kepler orbit in highly excited atomic states, can be exactly reduced to a problem of classical electron in a field of proton (hydrogen atom) and a circular microwave field (Benvenuto et al., 1994). In this formulation an electron with zero orbital momentum still have an energy change similar to that in Eq.(1) ($$J$$ remains constant at $$q \rightarrow 0$$) and its evolution is described by the Kepler map as it is shown in (Casati et al., 1987,Casati et al., 1988,Gontis and Kaulakis, 1987). We discuss mainly the map description for the circular restricted three-body problem but the analysis performed for the microwave ionization of hydrogen atoms shows that the Kepler map provides a good description also for the case of elliptic orbits or elliptic microwave polarization (Casati et al., 1988, Benvenuto et al., 1996).

It is important to note that for $$q \gg a_p$$, the first sine harmonic in the kick gives the main contribution while other harmonics have additional exponentially small factors. In this regime $$q$$ is approximately conserved for the planar case due to conservation of the Jacobi constant and the Tisserand relation (Shevchenko, 2015).

The Kepler map can be locally reduced to the Chirikov standard map and in this way the chaos border can be easily obtained. In order to do that in Eq.(1) the phase equation is linearized near the integer resonant values $$w_r = 1/r^{2/3}$$ determined by the condition of phase change by an integer $$\bar{x}-x= r$$ so that the phase equation becomes linear in rescaled energy $$\bar{\phi}=\phi - 3\pi (\bar{w}-w_r)/{w_r}^{5/2}$$ and the map is reduced to the Chirikov standard map

$\tag{3} \bar{y} = y + K \sin \phi ,\, \bar{\phi} = \phi + \bar{y} \;\; ,$ Figure 4: Contribution of the eight planets to the kick function $$F(x)$$ in Eq.(1) for the case of comet Halley. On each panel, $$F_i(x_i)$$ is obtained from the Melnikov integral calculation (black curve), the red dashed curve (the blue dotted-dashed curve) shows the Keplerian planet contribution (the Sun dipole contribution). Here $$x_i = t_i/T_i (\mod 1)$$ are the phases of the orbital motion a planet with a period $$T_i$$. On Jupiter and Saturn panels the kick functions extracted by Fourier analysis from observations and exact numerical computations are shown by red dots taken from Chirikov and Vecheslavov, 1989 (from Rollin et al., 2015a)

with the chaos parameter $$K = 3\pi J/{w_r}^{5/2}$$ (and $$y= -3\pi (w-w_r)/{w_r}^{5/2}$$). The integer resonances corresponding to $$r=1,2,3...$$ are well seen in the phase space shown in Fig.1. For the standard map the chaos becomes global at $$K > K_c \approx 0.97$$ that determines the chaos border for the Kepler map being $$w_c \approx 2.5 J^{2/5} < 1$$ (for parameters of Fig.1 this gives $$w_c \approx 0.34$$ in a good agreement with numerical data). Assuming that in the global chaotic domain the phases take random uncorrelated values one finds a diffusive growth $$w^2 \approx D \tau$$ with a diffusion constant $$D \approx J^2/2$$ (here $$\tau$$ is expressed in number of the particle orbital periods). This diffusion leads to the particle escape and determines a typical life time scale of the particle inside the system $$\tau_D \approx {w_c}^2/D \approx 12/J^{6/5}$$ expressed in number of map iterations. The transformation to real time can be estimated at $$t_D \sim T_p \tau_D/(2{w_c}^3)$$ but more exact computations are rather involved due to presence of stability islands, phase correlations and the particle period variation near the escape border $$w \rightarrow 0$$.

The Lyapunov exponent $$\Lambda$$ for the Kepler map can be computed numerically from the map iterations, like for the Chirikov standard map where $$\Lambda \approx \ln (K/2)$$ at $$K > 4$$ (see Chirikov, 1979). The theory and numerical results are described in Shevchenko, 2007. The computation of $$\Lambda$$ is done on a number of orbital periods being smaller than the diffusion escape time of light body from the binary. In the physical time one should take into account that the orbital period increases at $$w \rightarrow 0$$.

Numerically, the Lyapunov exponents are computed until the current Lyapunov exponent saturates at some plateau (emerging in the time dependence of the current exponent). For the Kepler map, the saturation is quick and takes place usually before the particle is ejected. Analytically, the Lyapunov exponents can be estimated either locally (at a given energy of a test particle, as approximated by the Chirikov standard map), or globally (averaged over the whole chaotic layer). For rough estimates, one may take the averaged Lyapunov exponent for the map (Chirikov's constant, see Shevchenko, 2007), divide it by the particle's current orbital period (which is just the time interval corresponding to one iteration of the map), and thus obtain an estimate of the local Lyapunov exponent in natural time units (yr$^{-1}$).

# Comet Halley dynamics

Detailed observational data exist for comet Halley due to various historical records dating back to year 240 B.C. (WikiHalley, 2017). Extensive numerical simulations performed by Yeomans and Kiang, 1981 determined 46 moments of time of perihelion passages of comet Halley up to year 1404 B.C. Using these 46 time moments and assuming that they are given by the Kepler map (1), with a certain kick function $$F(x)$$, Chirikov and Vecheslavov, 1989 determined this effective function expanding it in Fourier harmonics of Jupiter and Saturn periodic motion. Their empirical numerical points are show in Fig.4 (top two right panels). Even if such an approach is only approximate it demonstrates that in this Kepler map approximation the dynamics of comet is chaotic (see Fig.1 (top)) and its life time (in a bound orbit around the Sun) is relatively short being approximately only of the order of 10 million years. This time is obtained by averaging over many trajectories with initial orbital elements being close to the actual ones of comet Halley. This value is larger than the diffusive estimate for comet Halley life time $$t_d \sim t_H (w_c)^2/ D$$ where $$t_H \approx 75 \mathrm{yr}$$ is the current orbital period of comet Halley, $$w_c \approx 0.5$$ and $$D \approx (F_{max})^2/3 \approx 1.3 \times 10^{-5}$$ with the maximal kick amplitude $$F_{max} \approx 6 m_p/M \approx 0.006$$, thus giving $$t_H \approx 1.5 \times 10^6 \mathrm{yr}$$ (see Chirikov and Vecheslavov, 1989). As discussed above the difference is related to phase correlations and significant increase of comet period in a vicinity of escape border. The main contribution to the diffusion is given by Jupiter since for Saturn the maximal kick amplitude is by a factor 6 smaller while the diffusion rate is quadratic in amplitude.

The numerical computation of $$F(x)$$, taking into account only Jupiter in a circular orbit, showed visible deviations (see Fig.1a in Lages and Shepelyansky, 2013) from the kick function obtained in Chirikov and Vecheslavov, 1989. It turns out that here is present also a contribution from the Sun orbiting around the barycenter in a near circular orbit of radius $$a_S \approx m_p/M$$. Even if $$a_S$$ is small the Solar mass is large and this creates a dipole kick contribution comparable with the contribution given by a planet. These two contributions from the 8 planets and the Sun are shown in Fig.1 and described in detail in Rollin et al., 2015a. The kick from the Sun always has a sine shape since $$q \gg r_s$$ while a planet contribution has a more complicated shape of $$F(x)$$ since for Comet Halley its perihelion distance is comparable or even smaller than radii of planetary orbits (inclination of the comet orbital plane precludes close encounters). While the complete form of the kick function, taking into account the barycenter motion of the Sun and a planet, was analyzed by different groups (see e.g. Dvorak and Kribbel, 1990,Emelyanenko, 1992,Malyshkin and Tremaine, 1999,Zhou et al., 2000) the analysis performed in Rollin et al., 2015a allows one to have a clear physical insight on the two type of contributions in the kick function. Figure 5: Left panel: dependence of capture cross-section $$\sigma$$ on DMP energy $$w$$ in the continuum for the Jupiter case, the dashed lines show dependence $$\sigma \propto 1/w; 1/w^2$$; $$\sigma_p = \pi {a_p}^2$$. Right panel: dependence of rescaled captured number of DMP (per time unit) on $$w$$ for Jupiter, Saturn and a model planet with $$m_p/M=0.004$$ (solid, dashed and dot-dashed curves, respectively)(from Lages and Shepelyansky, 2013).

The Kepler map for Comet Halley provides understanding of the main physical properties of its chaotic dynamics. At the same time the map gives only an approximate description since the orbital momentum exhibits slow variation in time leading to change of $$q$$ as it is demonstrated in Dvorak and Kribbel, 1990 and Rollin et al., 2015a. The slow variation of the orbital momentum and kick amplitude was also analyzed for the microwave ionization of hydrogen atoms (see Casati et al., 1988). However, the comparison of the results of the Kepler map (in quantum regime) with the experimental data of Koch group shows that the map gives a good description of experimental data for real three-dimensional atoms (Casati et al., 1990). This comparison supports the expectation that the Kepler map gives a good average description of the long-term cometary dynamics in many cases. Another restriction for the map description of comet Halley evolution is related to the evaporation of its nucleus with time (see Chirikov and Vecheslavov, 1989,WikiHalley, 2017).

In the above description of comet Halley it is considered as a mathematical light body. For real physical comet, other effects, like evaporation during a passage through the Sun perihelion, play an important role and it is expected that the comet will evaporate or will split in the next few tens of thousands of years (WikiHalley, 2017). Figure 6: Left panel: radial density $$\rho \propto dN/dr$$ at present time $$t_S=4.5 \times 10^9 \mathrm{yr}$$ for the Solar System averaged over all captured DMP (normalization is fixed by $$\int_0^{6a_p}\rho dr =1$$). Right panel: volume density $$\rho/r^2$$ from the left panel data, the dashed line shows slope -2, the horizontal line shows average density for $$a_p/5 \leq r \leq a_p$$ (from Lages and Shepelyansky, 2013).

# Capture of dark matter by the Solar System

It is now well accepted that most of matter in the Universe is composed by dark matter which has only gravitational interaction (see e.g. Bertone et al., 2005). The understanding of the capture process of dark matter by the Solar System represents a significant interest. It is assumed that in a vicinity of the Solar System the velocity distribution of dark matter particles (DMP) has a Maxwell form $$f(v) dv = \sqrt{54/\pi} v^2/u^3 \exp(-3v^2/2u^2) dv$$ with the average module velocity $$u \approx 220 \mathrm{km/s}$$. The Galactic DMP mass density is expected to be $$\rho_g \sim 4 \times 10^{-25} \mathrm{g/cm}^3$$ (see Bertone et al., 2005). The first estimates for the DMP capture cross-section $$\sigma_p$$ and DMP dynamics, based on the Kepler map description, are presented in Khriplovich and Shepelyansky, 2009 with the further detailed analysis developed in Lages and Shepelyansky, 2013. The main effect is given by the interactions with the largest planet Jupiter. In the Kepler map picture only DMP with energies $$|w|= v^2/{v_p}^2 < w_{cap} = F_{max} \sim m_p/M$$ are captured under the condition that $$q < a_p$$ (see Fig.4, Jupiter case). The value of $$q$$ is determined by the DMP parameters at infinity, where its velocity is $$v$$ and its impact parameter is $$r_d$$, and hence $$q=(v r_d)^2/2kM$$ where $$k$$ is the gravitational constant. Since $$q \sim a_p$$ this leads to the DMP capture cross-section $\tag{4} \sigma \approx 3 \sigma_p M w_{cap}/(w m_p) \;, \; \sigma_p = \pi {a_p}^2 \; , \; w=v^2/{v_p}^2 \; , \; w_{cap} = m_p/M \; ,$

where the numerical coefficient is taken from the numerical simulations based on the Kepler map description of the DMP dynamics in Jupiter case (see details in Lages and Shepelyansky, 2013). The numerical results confirm the analytical estimate as it is shown in Fig.5. For typical capture velocities $$v_{cap} \sim v_p \sqrt{m_p/M} \sim 1 \mathrm{km/s}$$ the capture cross-section is by 1000 times larger than the area of Jupiter orbit. Thus the close encounters between DMP and a planet do not play significant role. This is also confirmed by a rapid drop of $$\sigma \approx 3 \sigma_p M {w_{cap}}^2/(m_p w^2)$$ for $$w > w_{cap}$$ as shown in Fig4. An additional confirmation is given by a numerical computation of the differential number of DMP captured per time unit $$d N = \sigma(w) n_g v_p^2 f(w) d |w| /2$$ where $$n_g=\rho_g/m_d$$ is the DMP density and $$f(w)$$ is the velocity distribution function given above. A number of DMP crossing the planet orbit area per time unit is $$N_p= \int_0^{1} n_g \sigma_p v_p^2 f(w) d|w|/2$$. The dependence $$dN/N_p dw$$, shown in Fig.5, drops quadratically for $$|w| > w_{cap}$$ showing that the contribution of close encounters is small. The results obtained with the Kepler map description are in agreement with the extensive numerical simulations of Newton dynamics of DMP in the Sun-Jupiter system (Peter, 2009), confirming that the close encounters do not give significant contribution for DMP capture and giving similar typical velocity for capture being $$1\mathrm{km/s}$$ as given above. Figure 7: Density of captured DMP in the Solar System at present time $$t_S$$. Top panel: DMP surface density $$\rho_s \propto dN/dzdr_\rho$$ shown on left in the cross plane $$(0,y,z)$$ perpendicular to the Jupiter orbit (data are averaged over $$r_\rho = (x^2+y^2)^{1/2}$$), on right in the Jupiter plane $$(x,y,0)$$; only the range $$|r|\leq 6a_p$$ around the Sun is shown. Bottom panels: corresponding DMP volume density $$dN/dxdydz$$ on left in the plane $$(0,y,z)$$, on right in the Jupiter orbit plane; only the range $$|r|\leq 2a_p$$ around the Sun is shown. Color grades signify the levels of density with yellow/black for maximum/zero density (from Lages and Shepelyansky, 2013). Figure 8: Logarithm of DMP density enhancement factor $$\log_{10}\zeta_g$$ shown by color and log-value-levels as a function of $$u/v_p$$ and $$J$$; two points are for $$J=0.005, u/v_p=17$$ (the Sun-Jupiter system) and $$u/v_p=0.035$$ (SMBH with $$v_p$$ of about $$2\%$$ of the speed of light) (from Rollin et al., 2015b).

# Chaotic dynamics of dark matter in the Solar System

The analysis of the long term DMP dynamics on the timescale of the Solar system age $$t_s \approx 4.5 \times 10^9 \mathrm{yr}$$ was performed in numerical simulations based on the generalized Kepler map (Lages and Shepelyansky, 2013). To determine the number of captured DMP $$N_{cap}(t)$$, in the Sun-Jupiter system, as a function of time a constant flow of scattered particles is modeled numerically with energy distribution $$d N_s \propto v f(v) dv$$ per time unit. The injection, capture, evolution and escape of DMP is described by the map (1) with multiple values of scattering parameters $$q, \theta, \varphi$$. The corresponding $F(x)$ function is computed numerically for each set of the scattering parameters with the scattering DMP distribution $$dN_s \propto dq d w$$. The scattering and evolution processes are followed for the whole time interval equal to the age of the Solar system with the total number of DMP, injected during time $$t_{S}$$ in the whole energy range $$0 \leq |w| \leq \infty$$ being $$N_{tot} \approx 1.5 \cdot 10^{14}$$ and with $$N_H =4 \cdot 10^9$$ scattered DMP in the Halley comet range $$0 < |w| \leq w_H \approx 5m_p/M$$. The efficiency of the map simulations allows one to consider the number of DMP five orders of magnitude larger than in the simulations of Newtonian dynamics in (Peter, 2009).

The results of this Kepler map modeling (Lages and Shepelyansky, 2013) show that at the initial stage the scattered DMP are captured in the Solar system and their number grows approximately linearly with time. However, after time scale $$t_H \approx 10^7 \mathrm{yr}$$ DMP start to escape from the Solar System and the steady state distribution of DMP is established. The dynamics of captured DMP is chaotic. The average DMP density distribution at time $$t_S$$ is shown in Fig. 6. It is interesting to note that at intermediate distances from the Sun there is a moderate growth of radial density $$\rho(r)$$ being similar to the velocity curve dependence $$v_m \propto r^{0.35}$$ observed in galaxies (see Fig.7 and Eqs.(1),(2) in Rubin et al., 1980). Indeed, in the virial theorem $${v^2}_m(r) \sim M(r)/r \sim \rho(r)$$, where $$M(r)$$ is the mass inside radius $$r$$. However, for the Sun-Jupiter system the total captured DMP mass is small (see below). The distributions of surface and volume DMP densities at time $$t_S$$ are shown in Fig.7; note visible specific regions of high DMP density.

The total mass of captured DMP inside the radius of Neptune orbit ($$r \leq 6a_p$$) is found to be rather small $$M_{r6} \approx 10^{-18} M$$. A typical average density inside the radius of Jupiter orbit is $$\rho_J \approx 10^{-4} \rho_g$$ with the corresponding mass $$M_{r1} \approx 5 \times 10^{-20} M$$. Thus, the density of captured DMP is much smaller than the galactic DMP density. However, it is by a factor $4 \times 10^3$ larger than the equilibrium DMP galactic density $$\rho_{gH} \approx 0.25 \rho_g/\kappa^{3/2} \approx 1.4 \cdot 10^{-32}\mathrm{g/cm^3}$$ taken in the energy range $$0<|w|<w_H = 5m_p/M$$ (with $$\kappa \approx 2u^2/(3v_p^2 w_H) \approx 4 \times 10^4$$) (Lages and Shepelyansky, 2013).

The main reason due to which there is no global density enhancement is related to the fact that the orbital velocity of Jupiter is small compared to the average velocity in the case of Galactic DMP Maxwell distribution $$u \approx 17 v_p$$. In such a case only a small fraction of DMP can be captured by one kick. In the case of a supermassive black hole (SMBH) with a star orbiting around it, it is possible to have situations with $$u \ll v_p$$ so that about a half of DMP flow is captured by one kick. The dependence of the density enhancement factor $$\zeta_g$$ on system parameters is analyzed in Rollin et al., 2015b using the generalized Kepler map simulation of $$10^{16}$$ DMP scattering in the energy range $$0 <|w|<\infty$$. The factor enhancement factor $$\zeta_g$$ is defined as the ratio of the DMP density in a binary system at distance $$r=a_p$$ to the Galactic density $$\rho_g$$. The results for $$\zeta_g$$ are well described by a semi-empirical formula, namely, $$\zeta_g = A \sqrt{J} (v_p/u)^3/[1+B J (v_p/u)^2]$$ with $$J=5 m_p/M$$, $$A =15.5$$ $$B=0.7$$ shown in Fig.8. This formula shows that the DMP density can be enhanced significantly in binary systems reaching the enhancement values $$\zeta_g \sim 10^4$$. Simple considerations explaining this formula are given in Rollin et al., 2015b.

Of course, the Kepler map description of DMP dynamics is an approximate one assuming that the orbital momentum, and hence the perihelion distance $$q$$, are conserved. But the obtained DMP density distributions (see Figs.5,6) are averaged over all scattering angles and momentum so that the averaged distributions are supposed to be close to those given by exact Newtonian dynamics. What is more, the Kepler map gives an estimate of typical energies of escaping DMP being determined by the last kick with $${v_d}^2 \sim F_{max} {v_p}^2 \sim m_p {v_p}^2/M$$ reaching high escape velocities $$v_d > 10^2 u$$ for SMBH with a star orbiting around with parameters $$m_p/M \sim 0.01$$, $$v_p \sim c/3$$, where $$c$$ is the speed of light (Yu and Tremaine, 2003,Rollin et al., 2015b). Such an ejection mechanism may play an important role in the evolution of binary black holes.

# Chaotic zones around rotating small bodies

Small irregular bodies such as asteroids and cometary nuclei, with consequently complex gravity fields (see Scheeres, 2012), quite often have bi-lobed shapes which can be equivalently modeled by either a rotating dumbbell (i.e., a massless rigid rod of length $$d$$ connecting two point masses $$m_1$$ and $$m_2$$) or by a rotating contact binary (see Fig.9). The case of a passively gravitating particle orbiting a gravitating small body is then analogous to the above depicted case of a particle orbiting the Solar System modeled by the Sun-Jupiter binary. The Kepler map can be applied taking into account two main differences: the masses of the lobes are comparable, so $$\mu=m_2/(m_1+m_2)\sim1/2$$, and the small body rotation frequency $$\omega$$ is arbitrary and can be very different from the Keplerian frequency $$\omega_0$$. For these small bodies, most often rubble pile objects, the rotation frequency $$\omega$$ should be less than the centrifugal disruption threshold $$\omega_0$$. The rotation of the small body modifies the test particle energy as the particle passes at the perihelion. The corresponding Kepler map can be written $\tag{5} \bar{E}=E+\Delta E(\mu,q,\omega,\phi),\quad \bar{\phi}=\phi+2\pi\omega/|2\bar{E}|^{3/2}$

Generalizing the work of Roy and Haddow, 2003 to the case of a binary with an arbitrary rotation frequency, the leading terms giving the exchange of energy between the test particle and the small body are (Lages et al., 2017) $\Delta E(\mu,q,\omega,\phi)=W_1\sin(\phi)+W_2\sin(2\phi).$

The first (second) term of the above sum is obtained from the octupole (quadrupole) term of the small body gravitational potential multipole expansion; analytical calculations give $W_1(\mu,q,\omega)=\mu (1-\mu) (1-2\mu) 2^{1/4} \pi^{1/2} \omega^{5/2} q^{-1/4} \exp(-2^{3/2} \omega q^{3/2}/3)$ and $W_2(\mu,q,\omega)\simeq-\mu(1-\mu)2^{15/4}\pi^{1/2}\omega^{5/2}q^{3/4}\exp\left(-2^{5/2}\omega q^{3/2}/3\right).$ Figure 11: Stability diagram for the orbital dynamics around Ida (left panel) and around Itokawa (right panel): the chaotic domain is shown by the reddish area. Chaos is determined by computing the Lyapunov exponent $$\Lambda$$ for a trajectory with initial orbital elements $$(q,e)$$. On the left panel, the black dot marks the dynamical locus of Dactyl, the moonlet of Ida (from Lages et al., 2017) Figure 12: Extent of the central chaotic zone around a rotating bilobed small body (with $$\mu=1/2$$) as a function of its rotation frequency $$\omega$$. The red (blue) area corresponds to chaotic (regular) dynamics of the particle orbiting the small body. Left panel: the central chaotic zone border has been obtained analytically. Right panel: obtained by iterating the Kepler map (from Lages et al., 2017).

For the study of a test particle dynamics in the vicinity of a small body (i.e., for $$\omega<\omega_0$$, $$q$$ about few $$d$$ (few dumb-bell sizes) and any $$\mu$$), a rapid comparison shows that the second harmonic coefficient $$W_2$$ dominates over the coefficient $$W_1$$ (see Lages et al., 2017). Fig.10 gives the Poincaré section $$(E,\phi)$$ computed using the Kepler map (5) illustrating the possible dynamics of Dactyl around 243 Ida. Using the Chirikov criterion (see the Map derivation and properties section), it is possible to estimate the chaos border location, i.e., the energy delimiting the chaotic component around the separatrix $$E=0$$ (see the last invariant KAM curve in red around $$E_{\mathrm{cr}}\simeq-0.096\omega_0^2d^2$$ in Fig.10 bottom right panel). The chaos border $$E_{\mathrm{cr}}$$ can then be used to locate the critical curve $$e_{\mathrm{cr}}(q)$$ separating, in the $$(e,q)$$ plane, initial orbital elements leading to chaotic or regular trajectories. Stability diagrams are shown in Fig.11 for the dynamics of a test particle around asteroids 243 Ida (Fig.11 left panel) and 25143 Itokawa (Fig.11 right panel). The stability diagrams are obtained by computing, for each initial condition $$(e,q)$$, the largest Lyapunov exponent $$\Lambda$$ of the Kepler map (5) as described above. In the diagrams, the initial conditions $$(e,q)$$ for orbits with high values of $$\Lambda$$ are marked by red colors and those with low $$\Lambda$$ values are marked by blue colors, according to the color bar scale given in Fig.11.

Let us define the central chaotic zone around the small body as the zone in $$q$$ such as at any initial eccentricity the particle’s dynamics is chaotic. This zone can be constructed analytically from the roots of the $$e_{\mathrm{cr}}(q)$$ function for any rotation rate $$\omega$$. The extent of the chaotic zone around small bodies is shown in Fig.12 (left panel: analytical determination using the Chirikov criterion; right panel: numerical determination by iterating the Kepler map (5)). The extent of the chaotic zone around a rotating small body is more than twice greater for rotation rate as slow as $$\omega\simeq0.1\omega_0$$ than for the Keplerian rotation rate $$\omega=\omega_0$$. Consequently, the relative size of the chaotic zone is greater around a rotating small body than around any non-bound gravitating binary system of the same size.

# Related topics Figure 13: The mass parameter - orbital period ratio relationship for the bi-planet and circumbinary exoplanet systems (dots). The vertical dashed (magenta) line shows the theoretical threshold in $$\mu$$, equal to $$0.05$$. The vertical dotted (cyan) line corresponds to $$\mu = 0.02$$ (From Shevchenko, 2015).
• Chaotic zones around gravitating binaries: Petrosky, 1986 used the Kepler map theory to show that the energy width of a one-sided chaotic band in the vicinity of the perturbed parabolic orbit scales as the power $2/5$ of the mass parameter $$\Delta E_\mathrm{cr} = \left|E_\mathrm{cr} \right| = - E_\mathrm{cr} \propto \mu^{2/5},$$ if $\mu \ll 1$. The particles with $E \in (-\Delta E_\mathrm{cr}, 0)$ move chaotically. This equation represents the $\mu^{2/5}$ law. Based on this concept of a chaotic layer near the parabolic motion (playing the role of a separatrix), and using analytical expressions for the Kepler map parameter, a strictly analytical criterion for the disintegration of a gravitating triple can be derived. Namely, based on the Kepler map theory, a theoretical criterion was proposed in Shevchenko, 2015 to describe a zone of chaotic orbits around a gravitationally-bound binary (double star, double black hole, double asteroid). The circumbinary continuous chaotic zone, where all circumbinary orbits are chaotic irrespective of their initial eccentricity, appears above a certain threshold in the mass parameter $\mu$ (the mass ratio of the central binary). It appears due to overlapping of the orbital resonances corresponding to integer ratios $p:1$ between the orbital periods of the particle and the central binary. Inside this zone, an unlimited chaotic diffusion of a test particle in the eccentricity $e$, up to unity, takes place until the particle is ejected from the system. The value of the mass parameter $\mu$, above which such a chaotic zone is present universally, was estimated to be $\sim 0.05$ (Shevchenko, 2015). The observed diversity of orbital configurations of bi-planet and circumbinary exoplanetary systems complies with the existence of this threshold in $\mu$: thus Fig.13 shows the observed dependence between mass parameter of the central binary $\mu$ versus the ratio of orbital periods $T_\mathrm{out}/T_\mathrm{in}$ of the outer planet and the central binary, the bi-planet systems are all on the left of the vertical line (the theoretical threshold $\mu \approx 0.05$ for the appearance of the central chaotic zone), whereas the circumbinary systems are on the right. The complete lack of exosystems with $T_\mathrm{out}/T_\mathrm{in} < 5$ at $\mu > 0.05$ is obvious. This is consistent with the theory: at $\mu > 0.05$, the central chaotic zone is formed, where the particles with any initial eccentricity are subject to the unlimited chaotic diffusion, until they are ejected from the system.
• Chaotic diffusion: The generalized Kepler map finds applications to description of chaotic diffusion in the dynamics of comets and meteor streams (see, e.g., Emelyanenko, 1992, Zhou et al., 2000)
• Survival probability: The probability that a particle stay in the system without escape (particle maintaining $$w>0$$) decays as $$P(t) \propto 1/t^{2/3}$$ at large physical times, since the orbital period between kicks is proportional to $$1/w^{3/2}$$ and the measure of such orbits is proportional to $$w$$ so that the orbits near $$w \rightarrow 0$$ play a dominant role (see Borgonovi et al., 1988). When the time is measured in the number of orbital periods $$\tau$$ then the decay is algebraic with $$P(\tau) \propto 1/{\tau}^{1.5}$$. This slow decay is related to the sticking of trajectories near the vicinity of the stability islands and the critical invariant curves. This survival probability is directly related to the statistics of Poincaré recurrences in symplectic maps with divided phase space (see more details at Chirikov standard map).
• Atomic physics: The applications of the Kepler map for the classical and quantum processes of microwave ionization of Rydberg and excited hydrodgen atoms are described in Microwave ionization of hydrogen atoms.
• Quantum chaos for dark matter in the Solar System: it is shown that the quantum effects for chaotic dark matter dynamics become significant for dark matter mass ratio to electron mass being smaller than $$2 \times 10^{-15}$$; this problem is shown to have similarities with a microwave ionization of Rydberg atoms studied previously experimentally and analytically; below the above mass border an energy diffusion over Rydberg states of dark matter atom becomes exponentially localized in analogy with the Anderson localization in disordered solids; the life time of dark matter in the Solar System (composed of Sun and Jupiter) is determined in dependence on mass ratio in the localized phase and a few graviton ionization regime (analogous to a few photon ionization in atomic physics) (see details in Shepelyansky, 2017).

# Historical notes

The second equation of the Kepler map is just a reformulation of the Kepler third law, published by Kepler, 1619, hence the name of the map. For the three-body problem, Sun, planet and comet, the Kepler map for the dynamics of comets in nearly parabolic orbits was derived by (Petrosky, 1986) and was applied to describe the long-term dynamics of the comet Halley using a generalized Kepler map by Chirikov and Vecheslavov at the end of 1986 (Chirikov and Vecheslavov, 1989). Independently, the Kepler map was derived for the classical and quantum evolutions of microwave electron excitation in Rydberg states in Casati et al., 1987 and for solely classical evolution in Gontis and Kaulakis, 1987. Probably Chirikov was a reviewer of Petrosky, 1986 paper (thus he pointed out the results of Petrosky to Shepelyansky in September 1986 when he showed to him the Kepler map description of microwave ionization of hydrogen atoms) that stimulated him to elaborate the generalized map description of chaotic, long-term dynamics of comet Halley. The studies of Chirikov and Vecheslavov were also stimulated by the Russian Vega mission to Venus and comet Halley in 1986 (WikiVega, 2017). The term Kepler map was coined in (Casati et al., 1987) for the problem of Microwave ionization of hydrogen atoms. For a microwave ionization process and comet dynamics with a larger perihelion distance ($$q \gg a_p$$) the kick function has a sine form and the map was named the Kepler map, for $$q \sim a_p$$ the kick function contains several harmonics that appeared for the Halley comet case in (Chirikov and Vecheslavov, 1989) corresponding to the generalized Kepler map (in some cases the term for this case is simplified to the Kepler map). The chaos border and diffusion rate for microwave ionization of hydrogen atom for different microwave frequencies were first obtained using the Chirikov criterion of overlapped resonances in 1983 (Delone et al., 1983). The capture cross-section by a binary was first obtained analytically and confirmed numerically by Heggie, 1975; and the map description of DMP capture was first developed in Khriplovich and Shepelyansky, 2009 and Lages and Shepelyansky, 2013. The relation between the problems of rotating dipole and quadrupole for a gravitating binary and autoionization of molecular Rydberg states was pointed out in Benvenuto et al., 1994.