Issue 
A&A
Volume 606, October 2017



Article Number  A43  
Number of page(s)  9  
Section  Celestial mechanics and astrometry  
DOI  https://doi.org/10.1051/00046361/201731167  
Published online  05 October 2017 
Rotation of an oblate satellite: Chaos control
Astronomical Observatory, Jagiellonian University, Orla 171, 30244 Kraków, Poland
email: mariusz.tarnopolski@uj.edu.pl
Received: 13 May 2017
Accepted: 1 August 2017
Aims. This paper investigates the chaotic rotation of an oblate satellite in the context of chaos control.
Methods. A model of planar oscillations, described with the Beletskii equation, was investigated. The Hamiltonian formalism was utilized to employ a control method for suppressing chaos.
Results. An additive control term, which is an order of magnitude smaller than the potential, is constructed. This allows not only for significantly diminished diffusion of the trajectory in the phase space, but turns the purely chaotic motion into strictly periodic motion.
Key words: celestial mechanics / chaos / methods: numerical
© ESO, 2017
1. Introduction
After Voyager 2 (Smith et al. 1982) obtained highquality images of Hyperion (Bond 1848; Lassel 1848), the highly aspherical moon of Saturn, it became clear that it remained in an exotic rotational state (Klavetter 1989a,b; Black et al. 1995; Devyatkin et al. 2002; Hicks et al. 2008; Thomas 2010; Harbison et al. 2011). In a seminal paper, Wisdom et al. (1984) predicted that Hyperion rotates chaotically, analyzed the phase space of a model of planar rotation, showed that it has an unstable attitude, and computed the Lyapunov time to be about 10 times the orbital period (i.e. 10 × 21.28 d). Since then, the analyses of chaotic rotation of an oblate celestial body has become very common, regarding Hyperion, per se, as well as other solar system satellites.
For instance, Boyd et al. (1994) employed the method of close returns to a sparse data set of dynamical states of Hyperion simulated with Euler equations, and found that a time series spanning about 2.6 yr of data is sufficient to infer the temporary rotational state (chaotic/regular). These findings agree with the recent results of Tarnopolski (2015), who showed that to extract a maximal Lyapunov exponent from groundbased photometric observations of Hyperion’s light curve, at least one year of wellsampled data is required, but a threeyear data set would be desirable. Black et al. (1995) performed numerical simulations with the full set of Euler equations to model the longterm dynamical evolution, and found that the chaotic tumbling of Hyperion leads to transitions between temporarily regular and chaotic rotation with a period of hundreds of days up to thousands of years. It was also shown that the Voyager 2 images of Hyperion indicate that the motion was predictable at the time of the passage (Thomas et al. 1995).
Beletskii et al. (1996) examined a number of theoretical models, including the gravitational, magnetic, and tidal moments and analyzed the rotation in the gravitational field of two centers. The structure of the respective phase spaces was investigated with Poincaré surfaces of section. The stability of spinorbit resonances, with application to the solar system satellites, was inferred based on a series expansion of the equation of planar rotation (Celletti & Chierchia 1998, 2000; see also Celletti 1990). The Lyapunov spectra were exhaustively examined for several satellites (Shevchenko 2002; Shevchenko & Kouprianov 2002; Kouprianov & Shevchenko 2003, 2005), and the Lyapunov times spanned 1.5–7 orbital periods for Hyperion. An interesting possibility that Enceladus might be locked in a 1:3 librationalorbital secondary resonance was investigated using the model of planar rotation within the Hamiltonian formalism (Wisdom 2004). The dynamical stability was examined for a number of known satellites by Melnikov & Shevchenko (2010); in particular, the synchronous spinorbit resonance in case of Hyperion was confirmed to be unstable. The Hamiltonian formalism has also been employed in the research of secondary resonances (Gkolias et al. 2016), and has taken nonconservative forces into account (Gkolias et al. 2017). Finally, regarding the influence of a secondary body on the rotation of an oblate moon, Tarnopolski (2017) showed, using the correlation dimension and bifurcation diagrams, that the introduction of an additional satellite can change the rotation from regular to chaotic.
Welldefined methods to reduce chaos in physical systems have been known for a long time now (e.g. Ott et al. 1990; Pyragas 1992). In general, a chaos control scheme already demonstrated its usefulness in astrodynamical applications, for example in maneuvering the ISEE3/ICE3 satellite to reach the GiacobiniZinner comet in 1985 (Shinbrot et al. 1993). Regarding rotation of an oblate moon, the full set of modified Euler equations was investigated numerically with various methods in the context of chaos control (Tsui & Jones 2000) for a satellite with thrusters. Investigation of the MimasTethys system was successfully conducted by means of a Hamiltonian chaos control method (Khan & Shahzad 2008). While still rather futuristic, an ability to construct efficient control terms might prove to be important in future asteroid capture missions and mining attempts (Kargel 1994; Sonter 1997; Koon et al. 2000; LevasseurRegourd et al. 2006; Elvis et al. 2011).
The knowledge about the rotational state of a celestial body proved to be crucial in the 2011 comet Tempel 1 flyby of StardustNExT (Veverka et al. 2013). Tempel 1 was the target of the Deep Impact mission in 2005 (A’Hearn et al. 2005). The mission’s aim was to make the impactor collide with the comet to excavate a crater to allow investigation of its interior structure. The crater was not measured directly after the impact owing to a large cloud of dust that obscured the view of the orbiting spacecraft. Hence one of the objectives of StardustNExT was to image the site of the impactor’s collision with the object. The rotational state of the comet needed to be accurately predicted so that the site of interest would be visible to the spacecraft and well illuminated during the flyby. The rotational period was shown to be decreasing (Belton et al. 2011) and the time of arrival to the comet was adjusted by 8.5 h one year prior to the encounter (Veverka et al. 2013). While, unlike Hyperion, Tempel 1 is not rotating chaotically, it is a plain illustration of the importance of the rotational state of small solar system bodies, among which chaotic rotation is expected to be common (Jacobson & Scheeres 2011; Jafari Nadoushan & Assadian 2015).
Herein, a construction of a control term that reduces chaos substantially, down to a strictly periodic rotation, is presented. Using the Beletskii (1966) equation, the Hamiltonian setting of the problem is employed to investigate the prospects of chaos control within the framework of Ciraolo et al. (2004). Numerical examples are carried out with the parameters suitable for Hyperion, but the approach is general enough to be widely applicable and not restricted to only some ranges of parameter values; the latter is also illustrated with parameters typical for solar system asteroids.
This paper is organized in the following manner. In Sect. 2 the rotational model is introduced, and analytical solutions in some specific cases are presented for completeness. In Sect. 3 the Hamiltonian approach is discussed and the chaos control method is described. Section 4 shows that the method is able to suppress diffusion of a chaotic trajectory in the phase space effectively. Discussion and conclusions are gathered in Sect. 5. The computer algebra system Mathematica is applied throughout this paper.
2. Rotational model of an oblate moon
2.1. Equations of motion
The rotational equation of motion can be derived based on the following assumptions (Danby 1962; Goldreich & Peale 1966; Wisdom et al. 1984; Sussman & Wisdom 2001; Greiner 2010):

1.
The orbit of the satellite around its host planet is Keplerian witheccentricity e, i.e. the distance r between the satellite and its host planet is (with major semiaxis a = 1) (1)where f is the true anomaly given by (2)where the overdot denotes differentiation with respect to time (see Fig. 1 for the geometrical setting). The units are chosen so that the orbital period T = 2π, where a = 1 leads to the orbital mean motion n = 1.

2.
The satellite is modeled as a triaxial ellipsoid, with principal moments of inertia A<B<C.

3.
The spin axis is fixed and perpendicular to the orbit plane; it is aligned with the shortest physical axis (i.e. the axis related to the largest principal moment of inertia C).
With the oblateness defined as , the equation of motion takes the form (3)where θ is the angular orientation of the satellite relative to some arbitrary line in space (see Fig. 1).
Fig. 1 Rotational model of an oblate moon (H) orbiting around a central body (denoted with S); α is the angle between the long axis of the moon and the direction to the planet, and θ is the dynamical angle from Eqs. (3) and (4), where f is the true anomaly. 
It should be emphasized that while the last assumption is valid for most solar system satellites, as the angular momentum is assumed to be constant, it was shown that the chaotic state is attitude unstable in case of Hyperion (Wisdom et al. 1984). Therefore, the model investigated herein should be thought of as an illustrative first approximation. However, a more astrodynamically realistic setting is also examined further on.
Using the wellknown general relations in the twobody motion, where M is the mass of the planet, μ = GM, h = r^{2}ḟ, h^{2} = μl, l/r = 1 + ecosf (i.e. Eq. (1)), l = a(1−e^{2}) (Khan et al. 1998), one can transform Eq. (3) so that f is the independent variable (the famous Beletskii (1966) equation; see Appendix A), (4)which as a dynamical system ẋ = F(x) reads where the prime denotes differentiation with respect to f. This system is nonautonomous^{1}.
It should be stressed that the equation of motion in the temporal domain, i.e. Eq. (3), can also be obtained from a Hamiltonian (Sussman & Wisdom 2001; Flynn & Saha 2005; Celletti 2010) describing the rotational motion, with an auxiliary differential equation governing the evolution of the true anomaly, Eq. (2), as a constraint coming from the description of the orbit. Such a Hamiltonian describes a 1degree of freedom, nonautonomous system with implicit time dependence through 2πperiodic functions f = f(t) and r = r(t).
Fig. 2 Solution of Eq. (9). Dependence of panel a: the angular velocity θ′(f) and panel b: the orientation θ(f) on the true anomaly f. Constants of integration are C_{1} = 1, C_{2} = 0. 
2.2. Analytical instances
For e = 0 and ω^{2} = 0, i.e. for a circular orbit and spherical symmetry of the satellite, Eq. (3) leads to a uniform rotation, .
When e = 0 and ω^{2} ≠ 0, Eq. (2) is solved trivially by f(t) = t + f_{0}, and Eq. (3) becomes the pendulum equation, (6)where α = θ−f (see Fig. 1), which is solved by the Jacobi elliptic functions (Lowenstein 2012; Lara & Ferrer 2015). With ψ = 2α, (7)where is a constant of motion. In Eq. (7), 4ℰ is the total energy of the pendulum, while ω^{2} is the maximal value of the potential energy. For , the motion is a libration in the orbital plane, and it is a rotation for (Murray & Dermott 1999). For the border case, , the motion takes place on a separatrix, with an initial condition α(0) = 0, i.e. (8)which is an unstable trajectory. The separatrix consists of two branches that form a cross in the phase space; this point is an unstable stationary point . Asymptotically, every initial condition ends in this point.
The last simple case is ω^{2} = 0 and e ≠ 0. This again leads, via Eq. (3), to . But now, Eq. (4) takes the form (9)which allows us to investigate how the rotation of a satellite depends on its orbital location. Equation (9) is integrable and yields the angular velocity (10)which in turn can be also directly integrated (see Fig. 2), (11)where C_{1} and C_{2} are constants of integration.
2.3. Fourier series expansion
Equation (3) can be expanded into a Fourierlike series, formally valid for all ω^{2}, (12)which can be naturally averaged over one orbital period to give (13)i.e. the pendulum equation from Eq. (6), with , to which the discussion from Sect. 2.2 applies. Here, is a resonant variable.
The coefficients H(k/ 2,e) were calculated by Cayley (1859), and the first few of these coefficients can be found tabulated in (Goldreich & Peale 1966; Wisdom et al. 1984; Celletti & Chierchia 1998; Murray & Dermott 1999). For a set k, H(k/ 2,e) is a power series in e, with its dominant term ∝e^{ k−2 }. Thence, higher order terms are negligible in most cases (see also Celletti & Chierchia 2000). The coefficients are given by an integral formula (Murray & Dermott 1999; Beletskii 2001) (14)where (15)Truncations of the series allow us to infer the stability of the spinorbit resonances in the solar system (Celletti & Chierchia 2000). These resonances occur when , i.e. . Locations in the phase space of the main resonances can be obtained by an analysis of Eqs. (3) or (4) and are given in Wisdom et al. (1984) and Black et al. (1995).
Based on the Chirikov (1979) criterion (see also Lichtenberg & Lieberman 1992), the onset of chaos is predicted to occur when the 1:1 and 3:2 resonances overlap, which happens when the critical value of ω is . For Hyperion’s eccentricity, e = 0.1, ω^{R} ≈ 0.31. This was confirmed with numerical simulations (Wisdom et al. 1984; Tarnopolski 2017).
3. Hamiltonian formalism and chaos control
3.1. Hamiltonian
Using the EulerLagrange equations on a general Lagrangian (16)one finds that Eq. (4) is obtained for (17)With , the Hamiltonian can be obtained through the Legendre transformation ℋ = pθ′−ℒ as (18)The Hamiltonian ℋ is explicitly dependent on f, which has the meaning of time in this setting, hence it is not a constant of motion. In particular, when ω^{2} = 0, the Hamilton equations lead to p ∈ const. and to Eq. (10) with C_{1} = p, as should be expected.
3.2. Phase space volume conservation
Consider a generic Hamiltonian ℋ(q,p). The Liouville theorem states that the flow is conservative, ∇·F = 0, i.e. the volume is preserved as the system evolves with time. Now consider a general dynamical system ẋ = F(x). A subset of the phase space with initial volume V_{0} evolves according to the equation (19)If the system under consideration is Hamiltonian, then this relation reverts to the Liouville theorem via the Hamilton equations , . Thence, the flow in the (θ,p) phase space described by the Hamiltonian in Eq. (18) is obviously conservative.
Consider, however, the dynamical system in Eq. (5). Its phase space is (θ,θ′), and (20)This equation can be solved for V(f) as follows: (21)where V_{0} = V(0). This function is 2π periodic.
Owing to Eq. (20) the flow is not conservative for the system in Eq. (5). Indeed, the flow is conservative in canonical coordinates (θ,p) with p = (1 + ecosf)^{2}θ′ (see Sect. 3.1). The volume does not need to be conserved in other coordinates. In general, consider a rectangular area in a twodimensional phase space (q,p): V_{0} = q_{0}p_{0} = q_{t}p_{t} = V_{t} ≡ V(t). The relation between p and is , where ζ(t) is a scalar function. Then, the corresponding volume in the space is , and . , but , which is the form of Eq. (21).
3.3. Chaos control method
The method of Ciraolo et al. (2004) is employed and briefly described as follows. Consider an autonomous Hamiltonian with a 2ndimensional phase space, where ℋ_{0} is the Hamiltonian of an integrable system and ε is a multiplicative parameter (relatively small). The aim is to find a control term ℱ of order O(ε^{2}) such that the motion described by a new Hamiltonian with the control term is much less diffused in the phase space, i.e. it is much more ordered. This is accomplished in the next steps.
First, with A being the actions and ϕ the angles, is decomposed as (22)For this purpose, integers k = (k_{1},k_{2},...,k_{n}) need to be found. Next, the frequency vector is constructed, . With and ξ, one defines the following operators via their action on : where is the resonant and is the nonresonant part of .
For the purpose of this work, it is sufficient to describe the construction of ℱ in the case when . The control term is given by a series (24)where (25)where and {·,·} is the Poisson bracket. The first term of the series in Eq. (24) takes the form (26)The whole infinite series does not need to be calculated. In fact, truncating the series after the first term turns out to effectively suppress chaos. To additionally improve the performance, a scaling factor η is introduced multiplicatively, i.e. ℱ_{2} → ηℱ_{2}, and its optimal value is found numerically.
Fig. 3 Poincaré surfaces of section obtained by integrating Eq. (30) and recording the points (θ(f),θ′(f)) with a step of Δf = 2π. The highlighted plot corresponds to η = 0, i.e. no control term. The most prominent suppression of chaos is obtained with η = 3. 
4. Results
Consider the Hamiltonian from Eq. (18). It has 1.5 degrees of freedom, so one needs to introduce a second action, E, to make it autonomous. It takes then the form (27)One denotes A = (p,E), ϕ = (θ,f). The frequency vector is . The potential can be decomposed into the form in Eq. (22) with k = (k_{1},k_{2}) ∈ { (−2,1),(−2,2),(−2,3),(2,−1),(2,−2),(2,−3) } , and are simple monomials in e and ω. Thence, for example . Overall, ξ·k ≠ 0 for p ≠ c(1 + ecosf)^{2}θ′, , i.e. , which are major spinorbit resonances. It is assumed hereinafter that the trajectory is far from these resonances, which means that , . Equation (23a) gives a lengthy output, as does ℱ_{2} in terms of p (see Appendix B for explicit formulae). With the substitution p = (1 + ecosf)^{2}θ′ one gets a much more compact expression, which expanded into a zerothorder series takes the final form, (28)which is indeed O(ω^{4}), and the index in brackets emphasizes it is a zerothorder expansion of ℱ_{2}.
Fig. 4 Panel a: solution of Eq. (30) for the suppressed case η = 3. Panel b: its power spectrum in a linear scale and panel c in log scale. Panel d: a broadband spectrum for the chaotic case η = 0. 
The Hamiltonian with the control term is (29)where the scaling factor η was already introduced (compare with Sect. 3.3). Eventually the Hamilton equations yield (30)Equation (30) was integrated numerically for e = 0.1, ω = 0.89, f ∈ [ 0,5000 ], and η ∈ [ 0,9.9 ] with a step of 0.1. The case η = 0 corresponds to the lack of a control term, i.e. the Hamiltonian in Eq. (18) yielding Eq. (4) as the equation of motion, hence serves as a reference. Figure 3 shows the obtained Poincaré surfaces of section. The plot corresponding to the lack of control term is highlighted with the thick frame. All plots up to η = 2.1 do not qualitatively influence the surface of section, i.e. the diffusion was not diminished at all. Impressively, η = 2.2−2.4 suppresses chaos to a very small region of the phase space, forming closed curves. At a slightly higher η, chaotic motion bursts out, to be suppressed again for a wide range of the scaling factor, η = 2.7−3.7. The most prominent reduction of the chaotic behavior is observed for η = 3, where the motion is turned to periodic. To emphasize the extent of this reduction, Fig. 4 shows a power spectrum (a LombScargle periodogram; Lomb 1976; Scargle 1982) of the θ′ time series; the case η = 0 (chaotic) is displayed for comparison. When η is increased further, the motion turns back to chaotic with occasional signs of stickiness, for example η = 4.5 or 6.4. The spectacular suppression occurs again for η = 7.5,7.8,7.9,8.3,8.4; on the other hand, stickiness is manifested for η = 8.2 and 9.6. The stickiness occurring for these higher values of η is different than for such lower values: for η = 4.5 the points gather around the 1:1 spinorbit resonance (compare with Black et al. 1995), while for η = 8.2 the gathering happens near the upper separatrix. The maximal value of the potential is 0.218, while for the control term ℱ_{2(0)} (i.e. taking η = 1) it does not exceed 0.022; i.e. the control term is only about 10% of the potential. The value η ≳ 10 is useless since the control term is then of the same order as the potential . However, for η ~ 2.2 (compare with Fig. 3) it is still nearly an order of magnitude smaller and is sufficient to suppress the chaotic behavior of the investigated system almost entirely. The range η ∈ [ 0,2.2 ] was also sampled with a smaller step of 0.01, but the resulting surfaces of section did not reveal any signs of chaos suppression. Also the firstorder series expansion of ℱ_{2}, denoted ℱ_{2(1)} (not shown), was tested; it allowed chaos to be suppressed down to a periodic motion with η = 0.7, but its overall maximal value attained 1/3 of the maximal value of the potential , i.e. there was no improvement compared to the performance of ℱ_{2(0)}.
Finally, to examine a more astrodynamically realistic scenario, with ω^{2} = 0.25 and e = 0.005 the above computations were repeated. The eccentricity is the mean value for the solar system satellites (Tarnopolski 2017), and ω^{2} was decreased slightly (still characterizing a very oblate object) when compared to the value of Hyperion because in the employed method this parameter is treated as a small perturbation. The chaotic zone in case η = 0, while much more narrow, is diffused in a large portion of the phase space (not shown). For values η = 1.7−2.2, chaos is suppressed as spectacularly as previously demonstrated, in which case the ratio of the control term (for η = 1.7) to the potential is about 0.14, i.e. twice as small as when the parameters of Hyperion were employed.
5. Discussion and conclusions
The Hamiltonian ℋ in Eq. (27), yielding Eq. (4), was employed in Sect. 4 to construct a control term ℱ that is able to suppress the chaotic behavior, which forms a new Hamiltonian . The zerothorder first term of the series from Eq. (24), i.e. the zerothorder series expansion of Eq. (25), ℱ_{2(0)}, was constructed as described in Sect. 3.3, resulting in Eq. (28). An additional multiplicative scaling term η was introduced, and the resulting equation of motion in Eq. (30) was integrated with η ∈ [ 0,9.9 ]; η = 0 corresponds to the lack of control term and served as a reference. The results represented in Fig. 3 show that the diffusion in the phase space cannot only be diminished, but when η = 3 (meaning that the control term is an order of magnitde smaller than the potential ) the motion becomes exactly periodic, which is confirmed with the power spectrum in Fig. 4. With η = 2.2, the motion becomes quasiperiodic with a very small diffusion in the phase space. Therefore, chaos is successfully suppressed.
While the model of planar rotation is in fact not applicable to Hyperion (used herein as a demonstration because, for its parameters, the phenomenon investigated herein is prominently visible), it should more adequately describe other oblate solar system satellites (Melnikov 2001). This is indeed the case, as the diffusion of the chaotic trajectory, with parameters typical for a solar system asteroid, is suppressed to nearly periodic rotation. The control term is then only about 14% of the potential. For even smaller ω^{2}, which is a perturbative parameter, the suppression should be expected to be even more efficient.
References
 A’Hearn, M. F., Belton, M. J. S., Delamere, W. A., et al. 2005, Science, 310, 258 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Beletskii, V. V. 1966, Motion of an Artificial Satellite about Its Center of Mass (Jerusalem: Israel Program For Scientific Translations) [Google Scholar]
 Beletskii, V. V. 2001, Essays on the motion of celestial bodies (Basel: Birkhäuser) [Google Scholar]
 Beletskii, V. V., Pivovarov, M. L., & Starostin, E. L. 1996, Chaos, 6, 155 [NASA ADS] [CrossRef] [Google Scholar]
 Belton, M. J. S., Meech, K. J., Chesley, S., et al. 2011, Icarus, 213, 345 [NASA ADS] [CrossRef] [Google Scholar]
 Black, G. J., Nicholson, P. D., & Thomas, P. C. 1995, Icarus, 117, 149 [NASA ADS] [CrossRef] [Google Scholar]
 Bond, W. C. 1848, MNRAS, 9, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Boyd, P. T., Mindlin, G. B., Gilmore, R., & Solari, H. G. 1994, ApJ, 431, 425 [NASA ADS] [CrossRef] [Google Scholar]
 Cayley, A. 1859, MmRAS, 29, 191 [NASA ADS] [Google Scholar]
 Celletti, A. 1990, Z. Angew. Math. Phys., 41, 174 [CrossRef] [Google Scholar]
 Celletti, A. 2010, Stability and Chaos in Celestial Mechanics (New York: Springer) [Google Scholar]
 Celletti, A., & Chierchia, L. 1998, Planet. Space Sci., 46, 1433 [NASA ADS] [CrossRef] [Google Scholar]
 Celletti, A., & Chierchia, L. 2000, Celest. Mech. Dyn. Astr. 76, 229 [NASA ADS] [CrossRef] [Google Scholar]
 Chirikov, B. V. 1979, Phys. Rep., 52, 263 [NASA ADS] [CrossRef] [Google Scholar]
 Ciraolo, G., Chandre, C., Lima, R., Vittot, M., & Pettini, M. 2004, Celest. Mech. Dyn. Astr., 90, 3 [NASA ADS] [CrossRef] [Google Scholar]
 Danby, J. M. A. 1962, Fundamentals of Celestial Mechanics (New York: The MacMillan Company) [Google Scholar]
 Devyatkin, A. V., Gorshanov, D. L., Gritsuk, A. N., et al. 2002, Sol. Syst. Res., 36, 248 [NASA ADS] [CrossRef] [Google Scholar]
 Elvis, M., McDowell, J., Hoffman, J. A., & Binzel, R. P. 2011, Planet. Space Sci., 59, 1408 [NASA ADS] [CrossRef] [Google Scholar]
 Flynn, A. E., & Saha, P. 2005, AJ, 130, 295 [NASA ADS] [CrossRef] [Google Scholar]
 Gkolias, I., Celletti, A., Efthymiopoulos, C., & Pucacco, G. 2016, MNRAS, 459, 1327 [NASA ADS] [CrossRef] [Google Scholar]
 Gkolias, I., Efthymiopoulos, C., Pucacco, G., & Celletti, A. 2017, Commun. Nonlinear Sci. Numer. Simulat., 51, 23 [NASA ADS] [CrossRef] [Google Scholar]
 Goldreich, P., & Peale, S. 1966, AJ, 71, 425 [NASA ADS] [CrossRef] [Google Scholar]
 Greiner, W. 2010, Classical Mechanics, Systems of Particles and Hamiltonian Dynamics (Berlin Heidelberg: Springer) [Google Scholar]
 Harbison, R. A., Thomas, P. C., & Nicholson, P. C. 2011, Celest. Mech. Dyn. Astron., 110, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Hicks, M. D., Buratti, B. J., & Basilier, E. N. 2008, Icarus, 193, 352 [NASA ADS] [CrossRef] [Google Scholar]
 Jacobson, S. A., & Scheeres, D. J. 2011, Icarus, 214, 161 [NASA ADS] [CrossRef] [Google Scholar]
 Jafari Nadoushan, M., & Assadian, N. 2015, Nonlinear Dyn., 81, 2031 [CrossRef] [Google Scholar]
 Kargel, J. S. 1994, J. Geophys. Res., 99, 129 [CrossRef] [Google Scholar]
 Khan, A., & Shahzad, M. 2008, AJ, 136, 2201 [NASA ADS] [CrossRef] [Google Scholar]
 Khan, A., Sharma, R., & Saha, L. M. 1998, ApJ, 116, 2058 [CrossRef] [Google Scholar]
 Klavetter, J. J. 1989a, AJ, 97, 570 [NASA ADS] [CrossRef] [Google Scholar]
 Klavetter, J. J. 1989b, AJ, 98, 1855 [NASA ADS] [CrossRef] [Google Scholar]
 Koon, W. S., Lo, M. W., Marsden, J. E., & Ross, S. D. 2000, Chaos, 10, 427 [NASA ADS] [CrossRef] [MathSciNet] [PubMed] [Google Scholar]
 Kouprianov, V. V., & Shevchenko, I. I. 2003, A&A, 410, 749 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kouprianov, V. V., & Shevchenko, I. I. 2005, Icarus, 176, 224 [NASA ADS] [CrossRef] [Google Scholar]
 Lara, M., & Ferrer, S. 2015, Eur. J. Phys., 36, 055040 [CrossRef] [Google Scholar]
 Lassel, W. 1848, MNRAS, 8, 195 [NASA ADS] [CrossRef] [Google Scholar]
 LevasseurRegourd, A. C., Hadamcik, E., & Lasue, J. 2006, Adv. Space Res., 37, 161 [NASA ADS] [CrossRef] [Google Scholar]
 Lichtenberg, A. J., & Lieberman, M. A. 1992, Regular and Chaotic Dynamics (New York: Springer) [Google Scholar]
 Lomb, N. R. 1976, Ap&SS, 39, 447 [NASA ADS] [CrossRef] [Google Scholar]
 Lowenstein, J. H. 2012, Essentials of Hamiltonian Dynamics (Cambridge: Cambridge University Press) [Google Scholar]
 Melnikov, A. V. 2001, Cosm. Res., 39, 68 [NASA ADS] [CrossRef] [Google Scholar]
 Melnikov, A. V., & Shevchenko, V. V. 2010, Icarus, 209, 786 [NASA ADS] [CrossRef] [Google Scholar]
 Murray, C. D., & Dermott, S. F. 1999, Solar System Dynamics (Cambridge: Cambridge University Press) [Google Scholar]
 Ott, E., Grebogi, C., & Yorke, J. A. 1990, Phys. Rev. Lett., 64, 1196 [NASA ADS] [CrossRef] [MathSciNet] [PubMed] [Google Scholar]
 Pyragas, K. 1992, Phys. Lett. A, 17, 421 [NASA ADS] [CrossRef] [Google Scholar]
 Scargle, J. D. 1982, ApJ, 263, 835 [NASA ADS] [CrossRef] [Google Scholar]
 Shevchenko, I. I. 2002, Cosm. Res., 40, 296 [NASA ADS] [CrossRef] [Google Scholar]
 Shevchenko, I. I., & Kouprianov, V. V. 2002, A&A, 394, 663 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Shinbrot, T., Grebogi, C., Ott, E., & Yorke, J. A. 1993, Nature, 363, 411 [NASA ADS] [CrossRef] [Google Scholar]
 Smith, B. A., Soderblom, L., Batson, R. M., et al. 1982, Science, 215, 504 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Sonter, M. J. 1997, Acta Astronaut., 41, 637 [NASA ADS] [CrossRef] [Google Scholar]
 Sussman, G. J., & Wisdom, J. 2001, Structure and Interpretation of Classical Mechanics (Cambridge: MIT Press) [Google Scholar]
 Thomas, P. C. 2010, Icarus, 208, 395 [NASA ADS] [CrossRef] [Google Scholar]
 Thomas, P. C., Black, G. J., & Nicholson, P. D. 1995, Icarus, 117, 128 [NASA ADS] [CrossRef] [Google Scholar]
 Tarnopolski, M. 2015, Ap&SS, 357, 160 [NASA ADS] [CrossRef] [Google Scholar]
 Tarnopolski, M. 2017, Celest. Mech. Dyn. Astr., 127, 121 [NASA ADS] [CrossRef] [Google Scholar]
 Tsui, A. P. M., & Jones, A. J. 2000, Phys. D, 135, 41 [CrossRef] [Google Scholar]
 Veverka, J., Klaasen, K., A’Hearn, M., et al. 2013, Icarus, 222, 424 [NASA ADS] [CrossRef] [Google Scholar]
 Wisdom, J. 2004, AJ, 128, 484 [NASA ADS] [CrossRef] [Google Scholar]
 Wisdom, J., Peale, S. J., & Mignard, S. 1984, Icarus, 58, 137 [NASA ADS] [CrossRef] [Google Scholar]
Appendix A: Beletskii equation
From Fig. 1 it follows that θ = α + f; hence θ′ = α′ + 1 and θ′′ = α′′, where the prime denotes differentiation with respect to f. Substituting this into Eq. (4), (A.1)which after defining δ = 2α and reordering becomes the Beletskii (1966) equation, i.e. (A.2)
Appendix B: Explicit formulae
The potential from Eq. (27) in the form of Eq. (22) is as follows: (B.1)The term , i.e. (B.2)leads to the control term (B.3)which with the substitution p = (1 + ecosf)^{2}θ′ turns into (B.4)After a power series expansion in θ′ to a zeroth order, Eq. (28) is obtained.
All Figures
Fig. 1 Rotational model of an oblate moon (H) orbiting around a central body (denoted with S); α is the angle between the long axis of the moon and the direction to the planet, and θ is the dynamical angle from Eqs. (3) and (4), where f is the true anomaly. 

In the text 
Fig. 2 Solution of Eq. (9). Dependence of panel a: the angular velocity θ′(f) and panel b: the orientation θ(f) on the true anomaly f. Constants of integration are C_{1} = 1, C_{2} = 0. 

In the text 
Fig. 3 Poincaré surfaces of section obtained by integrating Eq. (30) and recording the points (θ(f),θ′(f)) with a step of Δf = 2π. The highlighted plot corresponds to η = 0, i.e. no control term. The most prominent suppression of chaos is obtained with η = 3. 

In the text 
Fig. 4 Panel a: solution of Eq. (30) for the suppressed case η = 3. Panel b: its power spectrum in a linear scale and panel c in log scale. Panel d: a broadband spectrum for the chaotic case η = 0. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.