EDP Sciences
Free Access
Volume 606, October 2017
Article Number A43
Number of page(s) 9
Section Celestial mechanics and astrometry
DOI https://doi.org/10.1051/0004-6361/201731167
Published online 05 October 2017

© ESO, 2017

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 × 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 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 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 semi-axis 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).

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

Open with DEXTER

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, μ = GM, h = r2, h2 = μl, l/r = 1 + ecosf (i.e. Eq. (1)), l = a(1−e2) (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 non-autonomous1.

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π-periodic functions f = f(t) and r = r(t).

thumbnail 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 C1 = 1, C2 = 0.

Open with DEXTER

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 + f0, 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 C1 and C2 are constants of integration.

2.3. Fourier series expansion

Equation (3) can be expanded into a Fourier-like 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 spin-orbit 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 Euler-Lagrange 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 C1 = 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 V0 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 V0 = 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 two-dimensional phase space (q,p): V0 = q0p0 = qtpt = VtV(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 2n-dimensional 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 = (k1,k2,...,kn) 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 non-resonant 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.

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

Open with DEXTER

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 = (k1,k2) ∈ { (−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 pc(1 + ecosf)2θ, , i.e. , which are major spin-orbit 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 zeroth-order series takes the final form, (28)which is indeed O(ω4), and the index in brackets emphasizes it is a zeroth-order expansion of 2.

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

Open with DEXTER

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 Lomb-Scargle 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 spin-orbit 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 first-order 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 zeroth-order first term of the series from Eq. (24), i.e. the zeroth-order 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 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 ω2, which is a perturbative parameter, the suppression should be expected to be even more efficient.


Writing Eqs. (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.


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

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

Open with DEXTER
In the text
thumbnail 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 C1 = 1, C2 = 0.

Open with DEXTER
In the text
thumbnail 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.

Open with DEXTER
In the text
thumbnail 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.

Open with DEXTER
In the text

Current usage metrics show cumulative count of Article Views (full-text 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 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.