Rotation of an oblate satellite: Chaos control

TL;DR

This paper explores controlling chaotic rotation of an oblate satellite by applying a small additive control term, transforming chaotic motion into periodic motion using Hamiltonian formalism.

Contribution

It introduces a novel control method that suppresses chaos in satellite rotation with minimal energy input, based on the Beletskii equation and Hamiltonian approach.

Findings

Control term effectively reduces phase space diffusion.

Chaotic motion is converted into strictly periodic motion.

Control method requires an order of magnitude smaller energy than the potential.

Abstract

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.

Equations

$$r=\frac{\left(1-e^2\right)}{1+e\cos f},$$

$$\dot{f}=\frac{{\left(1+e\cos f\right)}^2}{{\left(1-e^2\right)}^{3/2}},$$

$$\ddot{\theta }+\frac{{\omega }^2}{2r^3}\sin 2\left(\theta -f\right)=0,$$

$$(1+e\cos f)\frac{{\text{d}}^2\theta }{\text{d}{f}^2}-2e\sin f\frac{\text{d}\theta }{\text{d}f}-\frac{{\omega }^2}{2}\sin 2(f-\theta )=0,$$

$${\theta }^{\prime }$$

$${\eta }^{\prime }$$

$$\ddot{\alpha }+\frac{{\omega }^2}{2}\sin 2\alpha =0,$$

$$\dot{\psi }=\pm \sqrt{2\left(4\mathcal{E}+{\omega }^2\cos \psi \right)},$$

$$\psi (t)=4\arctan \left({\mathrm{e}}^{\pm \omega t}\right)-\pi ,$$

$$\left(1+e\cos f\right){\theta }^{\prime \prime }(f)-2e{\theta }^{\prime }(f)\sin f=0,$$

$${\theta }^{\prime }(f)=\frac{C_1}{{\left(1+e\cos f\right)}^2},$$

$$\begin{split}&\theta(f)=C_{1}\Bigg{[}\frac{e\sin f}{\left(e^{2}-1\right)\left(1+e\cos f\right)}+\\ &\hskip 42.67912pt\frac{2}{\left(\sqrt{1-e^{2}}\right)^{3/2}}\arctan\frac{\left(1-e\right)\tan\frac{f}{2}}{\sqrt{1-e^{2}}}\Bigg{]}+C_{2},\end{split}$$

$$\ddot{\theta }+\frac{{\omega }^2}{2}\sum _{k=-\infty }^{\infty }H\left(\frac{k}{2},e\right)\sin (2\theta -kt)=0,$$

$$\ddot{\gamma }+\frac{{\omega }^2}{2}H\left(\frac{k}{2},e\right)\sin 2\gamma =0,$$

$$H\left(\frac{k}{2},e\right)=\frac{1}{\pi (1-e^2{)}^{3/2}}\underset{0}{\overset{\pi }{\int }}(1+e\cos f)[k\tau (f)-2f]df,$$

$$\tau (f)=2\arctan \sqrt{\frac{1-e}{1+e}}\tan \frac{f}{2}-\frac{e\sqrt{1-e^2}}{1+e\cos f}\sin f.$$

$$\mathcal{L}=\frac{1}{2}G\left(\theta (f),f\right){\theta }^{\prime }(f{)}^2+F\left(\theta (f),f\right){\theta }^{\prime }(f)-\mathcal{V}\left(\theta (f),f\right),$$

$$\mathcal{L}(\theta ,{\theta }^{\prime })=\frac{1}{2}{\left(1+e\cos f\right)}^2{\theta }^{\prime 2}+\frac{{\omega }^2}{4}\left(1+e\cos f\right)\cos 2(f-\theta ).$$

$$\mathcal{H}(\theta ,p)=\frac{p^2}{2{\left(1+e\cos f\right)}^2}-\frac{{\omega }^2}{4}\left(1+e\cos f\right)\cos 2(f-\theta ).$$

$$\frac{\text{d}}{\text{d}t}\ln V(t)=\nabla \cdot \text{@vec}F.$$

$$\frac{\text{d}}{\text{d}f}\ln V(f)=\nabla \cdot \text{@vec}F=\frac{2e\sin f}{1+e\cos f}.$$

$$V(f)=\frac{V_0(1+e{)}^2}{(1+e\cos f{)}^2},$$

$$\mathcal{V}=\sum _{\text{@vec}k\in {\mathbb{Z}}^n}{\mathcal{V}}_{\text{@vec}k}\exp (\text{i}\text{@vec}k\cdot \mathbf{\phi }).$$

$$\Gamma \mathcal{V}=\sum _{\begin{array}{c}\text{@vec}k\in {\mathbb{Z}}^n\\ \mathbf{\xi }\cdot \text{@vec}k\ne 0\end{array}}\frac{{\mathcal{V}}_{\text{@vec}k}}{\text{i}\mathbf{\xi }\cdot \text{@vec}k}\exp (\text{i}\text{@vec}k\cdot \mathbf{\phi }),$$

$$\mathcal{R}\mathcal{V}=\sum _{\begin{array}{c}\text{@vec}k\in {\mathbb{Z}}^n\\ \mathbf{\xi }\cdot \text{@vec}k=0\end{array}}{\mathcal{V}}_{\text{@vec}k}\exp (\text{i}\text{@vec}k\cdot \mathbf{\phi }),$$

$$\mathcal{N}\mathcal{V}=\sum _{\begin{array}{c}\text{@vec}k\in {\mathbb{Z}}^n\\ \mathbf{\xi }\cdot \text{@vec}k\ne 0\end{array}}{\mathcal{V}}_{\text{@vec}k}\exp (\text{i}\text{@vec}k\cdot \mathbf{\phi }),$$

$$\mathcal{F}=\sum _{s=2}^{\infty }{\mathcal{F}}_s,$$

$${\mathcal{F}}_s=-\frac{1}{s}\left\{\Gamma \mathcal{V},\mathcal{V}\right\}$$

$${\mathcal{F}}_2=-\frac{1}{2}\left\{\Gamma \mathcal{V},\mathcal{V}\right\}.$$

$$\mathcal{H}=\underset{{\mathcal{H}}_0}{\underbrace{\frac{p^2}{2{\left(1+e\cos f\right)}^2}+E}}\underset{\mathcal{V}}{\underbrace{-\frac{{\omega }^2}{4}\left(1+e\cos f\right)\cos 2(\theta -f)}}.$$

$$\text{@vec}k=(k_1,k_2)\in \{(-2,1),(-2,2),(-2,3),(2,-1),(2,-2),(2,-3)\},$$

$$\begin{split}&\mathcal{F}_{2(0)}=-\frac{1}{288(1+e\cos f)}\Bigg{\{}\omega^{4}\sin 2(f-\theta)\Big{\{}2e\big{[}9\sin(f-2\theta)+\\ &\hskip 51.21504pt+\sin(3f-2\theta)\big{]}+9\sin 2(f-\theta)\Big{\}}\Bigg{\}},\end{split}$$

Full text

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

1 Introduction

After Voyager 2 (Smith et al., 1982) obtained high-quality 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\times 21.28\,{\rm 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 years 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 ground-based photometric observations of Hyperion’s light curve, at least one year of well-sampled data is required, but a three-year data set would be desirable. Black et al. (1995) performed numerical simulations with the full set of Euler equations to model the long-term 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 spin-orbit 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 librational-orbital 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 spin-orbit 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 non-conservative 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.

Well-defined 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 ISEE-3/ICE-3 satellite to reach the Giacobini-Zinner 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 Mimas-Tethys 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; Levasseur-Regourd 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 Stardust-NExT (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 Stardust-NExT 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\,{\rm 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. Sect. 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):

The orbit of the satellite around its host planet is Keplerian with eccentricity $e$, i.e. the distance $r$ between the satellite and its host planet is (with major semi-axis $a=1$)

[TABLE]

where $f$ is the true anomaly given by

[TABLE]

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\pi$, where $a=1$ leads to the orbital mean motion $n=1$. 2. 2.

The satellite is modeled as a triaxial ellipsoid, with principal moments of inertia $A<B<C$. 3. 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 $\omega^{2}=\frac{3(B-A)}{C}$, the equation of motion takes the form

[TABLE]

where $\theta$ is the angular orientation of the satellite relative to some arbitrary line in space (see Fig. 1).

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 well-known general relations in the two-body motion, where $M$ is the mass of the planet, $\mu=GM$, $h=r^{2}\dot{f}$, $h^{2}=\mu l$, $l/r=1+e\cos f$ [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],

[TABLE]

which as a dynamical system $\dot{\@vec{x}}=\@vec{F}(\@vec{x})$ reads

[TABLE]

where the prime denotes differentiation with respect to $f$. This system is non-autonomous.111Writing Eq. (2) and (3) together as a system of three autonomous first-order differential equations yields a so-called 1.5 degree of freedom system; however, only two equations—those from a second-order Eq. (3)—carry dynamical information.

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 1-degree of freedom, non-autonomous system with implicit time dependence through $2\pi$-periodic functions $f=f(t)$ and $r=r(t)$.

2.2 Analytical instances

For $e=0$ and $\omega^{2}=0$, i.e. for a circular orbit and spherical symmetry of the satellite, Eq. (3) leads to a uniform rotation, $\theta(t)=\dot{\theta}_{0}t+\theta_{0}$.

When $e=0$ and $\omega^{2}\neq 0$, Eq. (2) is solved trivially by $f(t)=t+f_{0}$, and Eq. (3) becomes the pendulum equation,

[TABLE]

where $\alpha=\theta-f$ (see Fig. 1), which is solved by the Jacobi elliptic functions (Lowenstein, 2012; Lara & Ferrer, 2015). With $\psi=2\alpha$,

[TABLE]

where $\mathcal{E}=\frac{\ddot{\alpha}}{2}-\frac{\omega^{2}}{4}\cos 2\alpha$ is a constant of motion. In Eq. (7), $4\mathcal{E}$ is the total energy of the pendulum, while $\omega^{2}$ is the maximal value of the potential energy. For $\mathcal{E}<\frac{\omega^{2}}{4}$, the motion is a libration in the orbital plane, and it is a rotation for $\mathcal{E}>\frac{\omega^{2}}{4}$ (Murray & Dermott, 1999). For the border case, $\mathcal{E}=\frac{\omega^{2}}{4}$, the motion takes place on a separatrix, with an initial condition $\alpha(0)=0$, i.e.

[TABLE]

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 $\left(\pm\pi,0\right)$. Asymptotically, every initial condition ends in this point.

The last simple case is $\omega^{2}=0$ and $e\neq 0$. This again leads, via Eq. (3), to $\theta(t)=\dot{\theta}_{0}t+\theta_{0}$. But now, Eq. (4) takes the form

[TABLE]

which allows us to investigate how the rotation of a satellite depends on its orbital location. Eq. (9) is integrable and yields the angular velocity

[TABLE]

which in turn can be also directly integrated (see Fig. 2),

[TABLE]

where $C_{1}$ and $C_{2}$ are constants of integration.

2.3 Fourier series expansion

Eq. (3) can be expanded into a Fourier-like series, formally valid for all $\omega^{2}$,

[TABLE]

which can be naturally averaged over one orbital period to give

[TABLE]

i.e. the pendulum equation from Eq. (6), with $\omega\rightarrow\omega\sqrt{|H(k/2,e)|}$, to which the discussion from Sect. 2.2 applies. Here, $\gamma=\theta-\frac{k}{2}t$ 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 $\propto 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)

[TABLE]

where

[TABLE]

Truncations of the series allow us to infer the stability of the spin-orbit resonances in the solar system (Celletti & Chierchia, 2000). These resonances occur when $\dot{\gamma}\approx 0$, i.e. $|\dot{\theta}-\frac{k}{2}|\ll\frac{1}{2}$. Locations in the phase space of the main resonances can be obtained by an analysis of Eq. (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 $\omega$ is $\omega^{R}=1/\left(2+\sqrt{14e}\right)$. For Hyperion’s eccentricity, $e=0.1$, $\omega^{R}\approx 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 Euler-Lagrange equations $\frac{\textrm{d}}{\textrm{d}f}\frac{\partial\mathcal{L}}{\partial\theta^{\prime}}=\frac{\partial\mathcal{L}}{\partial\theta}$ on a general Lagrangian

[TABLE]

one finds that Eq. (4) is obtained for

[TABLE]

With $p=\frac{\partial\mathcal{L}}{\partial\theta^{\prime}}=\left(1+e\cos f\right)^{2}\theta^{\prime}$, the Hamiltonian can be obtained through the Legendre transformation $\mathcal{H}=p\theta^{\prime}-\mathcal{L}$ as

[TABLE]

The Hamiltonian $\mathcal{H}$ 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 $\omega^{2}=0$, the Hamilton equations lead to $p\in{\rm const.}$ and to Eq. (10) with $C_{1}=p$, as should be expected.

3.2 Phase space volume conservation

Consider a generic Hamiltonian $\mathcal{H}(q,p)$. The Liouville theorem states that the flow is conservative, $\nabla\cdot\@vec{F}=0$, i.e. the volume is preserved as the system evolves with time. Now consider a general dynamical system $\dot{\@vec{x}}=\@vec{F}(\@vec{x})$. A subset of the phase space with initial volume $V_{0}$ evolves according to the equation

[TABLE]

If the system under consideration is Hamiltonian, then this relation reverts to the Liouville theorem via the Hamilton equations $\dot{q}=\frac{\partial\mathcal{H}}{\partial p}$, $\dot{p}=-\frac{\partial\mathcal{H}}{\partial q}$. Thence, the flow in the $(\theta,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 $(\theta,\theta^{\prime})$, and

[TABLE]

This equation can be solved for $V(f)$ as follows:

[TABLE]

where $V_{0}=V(0)$. This function is $2\pi$ periodic.

Owing to Eq. (20) the flow is not conservative for the system in Eq. (5). Indeed, the flow is conservative in canonical coordinates $(\theta,p)$ with $p=\left(1+e\cos f\right)^{2}\theta^{\prime}$ (see Sect. 3.1). The volume does not need to be conserved in other coordinates. In general, consider a rectangular area in a two-dimensional phase space $(q,p)$: $V_{0}=q_{0}p_{0}=q_{t}p_{t}=V_{t}\equiv V(t)$. The relation between $p$ and $\dot{q}$ is $\dot{q}=\frac{p}{\zeta(t)^{2}}$, where $\zeta(t)$ is a scalar function. Then, the corresponding volume in the $(q,\dot{q})$ space is $\tilde{V}_{0}=\frac{q_{0}p_{0}}{\zeta(0)^{2}}$, and $\tilde{V}_{t}=\frac{q_{t}p_{t}}{\zeta(t)^{2}}$. $\tilde{V}_{t}\neq\tilde{V}_{0}$, but $\tilde{V}_{t}\zeta(t)^{2}=\tilde{V}_{0}\zeta(0)^{2}\,\Rightarrow\,\tilde{V}_{t}=\frac{\tilde{V}_{0}\zeta(0)^{2}}{\zeta(t)^{2}}$, 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 $\mathcal{H}=\mathcal{H}_{0}+\varepsilon\mathcal{V}$ with a $2n$-dimensional phase space, where $\mathcal{H}_{0}$ is the Hamiltonian of an integrable system and $\varepsilon$ is a multiplicative parameter (relatively small). The aim is to find a control term $\mathcal{F}$ of order $O(\varepsilon^{2})$ such that the motion described by a new Hamiltonian with the control term $\tilde{\mathcal{H}}=\mathcal{H}+\mathcal{F}=\mathcal{H}_{0}+\varepsilon\mathcal{V}+\mathcal{F}$ is much less diffused in the phase space, i.e. it is much more ordered. This is accomplished in the next steps.

First, with $\@vec{A}$ being the actions and $\bm{\varphi}$ the angles, $\mathcal{V}$ is decomposed as

[TABLE]

For this purpose, integers $\@vec{k}=(k_{1},k_{2},\ldots,k_{n})$ need to be found. Next, the frequency vector is constructed, $\bm{\xi}=\frac{\partial\mathcal{H}_{0}}{\partial\@vec{A}}$. With $\mathcal{V}$ and $\bm{\xi}$, one defines the following operators via their action on $\mathcal{V}$:

[TABLE]

where $\mathcal{R}\mathcal{V}$ is the resonant and $\mathcal{N}\mathcal{V}$ is the non-resonant part of $\mathcal{V}$.

For the purpose of this work, it is sufficient to describe the construction of $\mathcal{F}$ in the case when $\mathcal{R}\mathcal{V}=0$. The control term is given by a series

[TABLE]

where

[TABLE]

where $\mathcal{F}_{1}=\mathcal{V}$ and $\left\{\cdot,\cdot\right\}$ is the Poisson bracket. The first term of the series in Eq. (24) takes the form

[TABLE]

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 $\eta$ is introduced multiplicatively, i.e. $\mathcal{F}_{2}\rightarrow\eta\mathcal{F}_{2}$, and its optimal value is found numerically.

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

[TABLE]

One denotes $\@vec{A}=(p,E)$, $\bm{\varphi}=(\theta,f)$. The frequency vector is $\bm{\xi}=\left(\frac{\partial\mathcal{H}_{0}}{\partial p},\frac{\partial\mathcal{H}_{0}}{\partial E}\right)=\left(\frac{p}{(1+e\cos f)^{2}},1\right)$. The potential $\mathcal{V}$ can be decomposed into the form in Eq. (22) with

[TABLE]

and $\mathcal{V}_{\@vec{k}}$ are simple monomials in $e$ and $\omega$. Thence, for example $\bm{\xi}\cdot(-2,1)=1-\frac{2p}{(1+e\cos f)^{2}}$. Overall, $\bm{\xi}\cdot\@vec{k}\neq 0$ for $p\neq c(1+e\cos f)^{2}\theta^{\prime}$, $c=\frac{1}{2},1,\frac{3}{2}$, i.e. $\theta^{\prime}\neq\frac{1}{2},1,\frac{3}{2}$, which are major spin-orbit resonances. It is assumed hereinafter that the trajectory is far from these resonances, which means that $\mathcal{R}\mathcal{V}=0$, $\mathcal{N}\mathcal{V}=\mathcal{V}$. Eq. (23a) gives a lengthy output, as does $\mathcal{F}_{2}$ in terms of $p$ (see Appendix B for explicit formulae). With the substitution $p=(1+e\cos f)^{2}\theta^{\prime}$ one gets a much more compact expression, which expanded into a zeroth-order series takes the final form,

[TABLE]

which is indeed $O(\omega^{4})$, and the index in brackets emphasizes it is a zeroth-order expansion of $\mathcal{F}_{2}$.

The Hamiltonian with the control term is

[TABLE]

where the scaling factor $\eta$ was already introduced (compare with Sect. 3.3). Eventually the Hamilton equations yield

[TABLE]

Eq. (30) was integrated numerically for $e=0.1$, $\omega=0.89$, $f\in[0,5000]$, and $\eta\in[0,9.9]$ with a step of 0.1. The case $\eta=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. Fig. 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 $\eta=2.1$ do not qualitatively influence the surface of section, i.e. the diffusion was not diminished at all. Impressively, $\eta=2.2\!-\!2.4$ suppresses chaos to a very small region of the phase space, forming closed curves. At a slightly higher $\eta$, chaotic motion bursts out, to be suppressed again for a wide range of the scaling factor, $\eta=2.7\!-\!3.7$. The most prominent reduction of the chaotic behavior is observed for $\eta=3$, where the motion is turned to periodic. To emphasize the extent of this reduction, Fig. 4 shows a power spectrum (a Lomb-Scargle periodogram; Lomb 1976; Scargle 1982) of the $\theta^{\prime}$ time series; the case $\eta=0$ (chaotic) is displayed for comparison. When $\eta$ is increased further, the motion turns back to chaotic with occasional signs of stickiness, for example $\eta=4.5$ or 6.4. The spectacular suppression occurs again for $\eta=7.5,7.8,7.9,8.3,8.4$; on the other hand, stickiness is manifested for $\eta=8.2$ and 9.6. The stickiness occurring for these higher values of $\eta$ is different than for such lower values: for $\eta=4.5$ the points gather around the 1:1 spin-orbit resonance (compare with Black et al. 1995), while for $\eta=8.2$ the gathering happens near the upper separatrix. The maximal value of the potential $\mathcal{V}$ is 0.218, while for the control term $\mathcal{F}_{2(0)}$ (i.e. taking $\eta=1$) it does not exceed 0.022; i.e. the control term is only about 10% of the potential. The value $\eta\gtrsim 10$ is useless since the control term is then of the same order as the potential $\mathcal{V}$. However, for $\eta\sim 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 $\eta\in[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 first-order series expansion of $\mathcal{F}_{2}$, denoted $\mathcal{F}_{2(1)}$ (not shown), was tested; it allowed chaos to be suppressed down to a periodic motion with $\eta=0.7$, but its overall maximal value attained 1/3 of the maximal value of the potential $\mathcal{V}$, i.e. there was no improvement compared to the performance of $\mathcal{F}_{2(0)}$.

Finally, to examine a more astrodynamically realistic scenario, with $\omega^{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 $\omega^{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 $\eta=0$, while much more narrow, is diffused in a large portion of the phase space (not shown). For values $\eta=1.7\!-\!2.2$, chaos is suppressed as spectacularly as previously demonstrated, in which case the ratio of the control term (for $\eta=1.7$) to the potential $\mathcal{V}$ is about 0.14, i.e. twice as small as when the parameters of Hyperion were employed.

5 Discussion and conclusions

The Hamiltonian $\mathcal{H}$ in Eq. (27), yielding Eq. (4), was employed in Sect. 4 to construct a control term $\mathcal{F}$ that is able to suppress the chaotic behavior, which forms a new Hamiltonian $\tilde{\mathcal{H}}=\mathcal{H}+\eta\mathcal{F}$. The zeroth-order first term of the series from Eq. (24), i.e. the zeroth-order series expansion of Eq. (25), $\mathcal{F}_{2(0)}$, was constructed as described in Sect. 3.3, resulting in Eq. (28). An additional multiplicative scaling term $\eta$ was introduced, and the resulting equation of motion in Eq. (30) was integrated with $\eta\in[0,9.9]$; $\eta=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 $\eta=3$ (meaning that the control term is an order of magnitde smaller than the potential $\mathcal{V}$) the motion becomes exactly periodic, which is confirmed with the power spectrum in Fig. 4. With $\eta=2.2$, the motion becomes quasi-periodic 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 $\omega^{2}$, which is a perturbative parameter, the suppression should be expected to be even more efficient.

Appendix A Beletskii equation

From Fig. 1 it follows that $\theta=\alpha+f$; hence $\theta^{\prime}=\alpha^{\prime}+1$ and $\theta^{\prime\prime}=\alpha^{\prime\prime}$, where the prime denotes differentiation with respect to $f$. Substituting this into Eq. (4),

[TABLE]

which after defining $\delta=2\alpha$ and reordering becomes the Beletskii (1966) equation, i.e.

[TABLE]

Appendix B Explicit formulae

The potential from Eq. (27) in the form of Eq. (22) is as follows:

[TABLE]

The term $\Gamma\mathcal{V}$, i.e.

[TABLE]

leads to the control term

[TABLE]

which with the substitution $p=\left(1+e\cos f\right)^{2}\theta^{\prime}$ turns into

[TABLE]

After a power series expansion in $\theta^{\prime}$ to a zeroth order, Eq. (28) is obtained.