Obserwatorium Astronomiczne
Oświadczenie
Ja, niżej podpisany: Mariusz Tarnopolski (nr indeksu: 1027879), doktorant Wydziału Fizyki, Astronomii i Informatyki Stosowanej Uniwersytetu Jagiellońskiego, oświadczam, że przedłożona przeze mnie rozprawa doktorska pt. „Chaotic rotational dynamics of Hype-rion” jest oryginalna i przedstawia wyniki badań wykonanych przeze mnie osobiście, pod kierunkiem dra hab. Zdzisława Goldy. Pracę napisałem samodzielnie.
Oświadczam, że moja rozprawa doktorska została opracowana zgodnie z Ustawą o prawie autorskim i prawach pokrewnych z dnia 4 lutego 1994 r. (Dziennik Ustaw 1994 nr 24 poz. 83 wraz z późniejszymi zmianami).
Jestem świadom, że niezgodność niniejszego oświadczenia z prawdą ujawniona w dowol-nym czasie, niezależnie od skutków prawnych wynikających z ww. ustawy, może spowodować unieważnienie stopnia nabytego na podstawie tej pracy.
Podpis doktoranta: . . . .
Astronomical Observatory
Mariusz Tarnopolski
Chaotic rotational dynamics of
Hyperion
A dissertation submitted to the Jagiellonian University
for the degree of Doctor of Philosophy
PhD Advisor:
dr hab. Zdzisław Golda
Pragnę wyrazić głęboką wdzięczność moim rodzicom, bez których żaden z moich sukcesów nie byłby możliwy. Chcę również podziękować mojemu promotorowi, drowi hab. Zdzisławowi Goldzie, za ciągłe wsparcie i cenne wskazówki.
I wish to express my deepest gratitude to my parents, without whom none of my successes would be possible. I want to thank my supervisor, dr hab. Zdzisław Golda, for his constant support and valuable tips.
Funding
During my PhD studies I received the following funding for my projects:
• 2015/2016 – the Minister of Science and Higher Education scholarship for outstanding achievements, grant nr D-0013/2015
• 2015/2016 – stipend from the Marian Smoluchowski CRACOW SCIENTIFIC CON-SORTIUM, in the framework of the KNOW scholarship, grant nr KNOW/48/
SS/PC/2015
• 2014/2015 – grant no. 7150/E-338/M/2014 from the Ministry of Science and Higher Education, for the project The study of nonlinear phenomena in Kepler catalogue of contact binaries, duration: 06/2014 – 04/2015 (position: Project Manager)
• 2013 – stipend from the Marian Smoluchowski CRACOW SCIENTIFIC CONSOR-TIUM, in the framework of the KNOW scholarship, grant nr KNOW/02/
SS/TM/2013
I want to thank also Fundacja Astronomii Polskiej and Fundacja Studentów i Absol-wentów UJ “Bratniak” for funding my participation in the conference Chaos 2015 in Paris, France.
This dissertation presents the results of the research conducted on the chaotic rotation of Hyperion. The main results are:
i) the construction—within a Hamiltonian formalism of the model of planar rotation—of a control term (an order of magnitude smaller than the potential) that turns chaotic trajectories into strictly periodic ones;
ii) analysis of the influence of the gravitational interaction with Titan on the planar rota-tional dynamics; it is found that Titan’s presence can change regular trajectories into chaotic ones;
iii) simulated lightcurves, in the framework of a three-dimensional rotation modeled with Euler equations, were analysed in order to explore conditions for Earth-based photo-metric observations to meet for estimating a maximal Lyapunov exponent (mLE). It was found that at least one year of dense data is required to detect a positive mLE—if present—in a reliable way.
The thesis is divided into three parts: Part I
Chapter 1 provides a general literature overview in the field of Hyperion’s rotation and relevant works about oblate satellites’ rotation in general.
Chapter 2 contains, for the sake of self-containedness, a brief description of the basic methods and tools of the theory of chaos that are exploited in the presented work. An overview of the implementions into software packages, present in the literature and employed in this work, is also provided.
Chapter 3 introduces the model of planar rotation that is the basis of the work presented herein, and shortly describes its basic properties. A three-dimensional rotation formalism, together with a setting for simulating Earth-based photometric lightcurves, is also presented. Part II
In Chapter 4, a Hamiltonian chaos control method is employed to the Beletskii (1966) equation. A control term an order of magnitude smaller than the potential is constructed. A spectacular suppression of chaos down to a strictly periodic rotation is achieved.
vi
Chapter 5 examines the influence of Titan’s gravitation on the rotational dynamics of Hyperion. The techniques of bifurcation diagrams and evolution of a cloud of initial con-ditions on the phase space are employed. It was found that perturbations (three orders of magnitude smaller than the interaction with the host planet) from an additional satellite cause a number of regular trajectories to become chaotic, hence enlarging the chaotic zone. Chapter 6 aims at examining photometric lightcurves of Hyperion, existent in the litera-ture as well as simulated using the full set of Euler equations, in order to list conditions that real-world astronomical observations should meet to extract a positive—if present—positive mLE in a reliable manner. It turned out that at least a one-year long ground-based obser-vational campaign from a single site is required for the chaoticity to be detected in the data gathered in this way.
Chapter 7 provides a summary of the discussed results. Part III
This part contains the publications that have been used as a basis of the dissertation. All of them have been included in the default journal format (in order of appearance in the List of Publications).
Keywords
Niniejsza rozprawa prezentuje wyniki badań nad chaotyczną dynamiką rotacyjną Hype-riona. Głównymi wynikami są:
i) konstrukcja, w ramach formalizmu hamiltonowskiego modelu rotacji w płaszczyźnie or-bity, członu sterującego (o rząd wielkości mniejszego niż potencjał) zmieniającego tra-jektorie chaotyczne w czysto periodyczne;
ii) analiza wpływu oddziaływania grawitacyjnego z Tytanem na dynamikę rotacji w pła-szczyźnie orbity; pokazano, że obecność Tytana może zmienić trajektorie okresowe w chaotyczne;
iii) krzywe zmian blasku, symulowane numerycznie w ramach modelu rotacji trójwymia-rowej opisanej równaniami Eulera, zostały przeanalizowane celem zbadania warunków jakie prowadzone z Ziemi obserwacje fotometryczne muszą spełnić, by móc wyznaczyć maksymalny wykładnik Lapunova (maximal Lyapunov exponent, mLE). Pokazano, że co najmniej jeden rok gęsto próbkowanych danych jest potrzebny by wyznaczyć dodatni mLE w wiarygodny sposób.
Niniejsza rozprawa została podzielona na trzy części: Część I
Rozdział 1 stanowi przegląd literatury dotyczącej rotacji Hyperiona oraz ważnych prac nt. rotacji asferycznych ciał niebieskich w ogólności.
Rozdział 2 zawiera krótkie wprowadzenie do metod i narzędzi teorii chaosu, które zostały użyte w niniejszej pracy, wraz z opisem implementacji zastosowanych pakietów obliczeniowych.
W rozdziale 3 zdefiniwany jest model rotacji w płaszczyźnie orbity oraz krótko omówione są jego własności. Opisany został również trójwymiarowy model rotacji, wraz z zastosowanym algorytmem symulacji krzywych zmian blasku.
Część II
Rozdział 4 przedstawia zastosowanie metody sterowania chaosem w ujęciu hamiltonow-skim do równania Bieleckiego (Beletskii 1966). Skonstruowano zaburzenie Hamiltonianu o rząd wielkości mniejsze niż potencjał. Zaobserwowano spektakularne stłumienie zachowań chaotycznych – w efekcie uzyskano trajektorie ściśle okresowe.
viii
Rozdział 5 przedstawia wyniki badań wpływu oddziaływania grawitacyjnego Tytana na rotację Hyperiona. Zastosowano diagramy bifurkacyjne oraz metodę śledzenie ewolucji zbiorów warunków początkowych w przestrzeni fazowej. Pokazano, że perturbacje pochodzące od drugiego satelity (trzy rzędy wielkości mniejsze niż oddziaływanie z Saturnem) sprawiają, że wiele trajektorii regularnych staje się chaotycznymi.
Rozdział 6 prezentuje wyniki analizy fotometrycznych krzywych zmian blasku Hyperiona obecnych w literaturze, zarówno opublikowanych rzeczywistych danych obserwacyjnych, jak i symulowanych przy użyciu pełnego układu równań Eulera. Celem tych analiz jest podanie warunków jakie obserwacje astronomiczne muszą spełnić, by móc wiarygodnie wyznaczyć mLE z uzyskanych szeregów czasowych. Okazuje się, że obserwacje naziemne powinny być prowadzone przez co najmniej jeden rok, by zachowania chaotyczne były widoczne w ze-branych danych.
Rozdział 7 zawiera podsumowanie otrzymanych wyników. Część III
W tej części zebrano publikacje autora, które stanowią podstawę przedstawionej rozprawy, w postaci w jakiej zostały opublikowane w poszczególnych czasopismach (w kolejności ich występowania w Liście Publikacji).
Słowa kluczowe
This dissertation has been written as a summary of the scientific activities reported in the following articles:
• Tarnopolski M., Rotation of an oblate satellite: Chaos control, submitted to A&A, arXiv:1704.02015 (2017)
• Tarnopolski M., Influence of a second satellite on the rotational dynamics of an oblate moon, Celestial Mechanics and Dynamical Astronomy, 127(2), 121–138 (2017); DOI: 10.1007/s10569-016-9719-7
• Tarnopolski M., Nonlinear time-series analysis of Hyperion’s lightcurves, Astrophysics and Space Science, 357:160 (2015); DOI: 10.1007/s10509-015-2379-3
Acknowledgements iii
Abstract v
Abstract in Polish (Streszczenie) vii
List of publications ix
Part I
Preliminaries
1
1 Introduction 3
1.1 Overview . . . 3
1.2 Aims and objectives of the thesis . . . 7
2 Methods 9 2.1 Poincaré surface of section . . . 10
2.2 Bifurcation diagrams . . . 11
2.3 Correlation dimension . . . 11
2.4 Lyapunov exponents . . . 12
2.4.1 Indicator of chaos . . . 12
2.4.2 Lyapunov spectrum . . . 13
2.5 Time series analysis . . . 14
2.5.1 Power spectrum . . . 14
2.5.2 Phase space reconstruction . . . 14
False nearest neighbours . . . 15
Mutual information . . . 16
2.5.3 Maximal Lyapunov exponents . . . 16
The Wolf et al. (1985) method . . . 16
The Kantz (1994) method . . . 17
2.5.4 Stationarity test . . . 17
3 Rotational model of an oblate moon 19 3.1 Planar rotation . . . 19
xii CONTENTS
3.2 Three-dimensional rotation . . . 22
3.2.1 Euler equations . . . 22
3.2.2 Lightcurve modeling . . . 23
Part II
Summary of the published articles
27
4 Chaos control 29 4.1 Chaos control . . . 294.2 Discussion and conclusions . . . 31
5 Titan’s influence 33 5.1 Rotational model . . . 33
5.2 Evolving clouds of initial conditions . . . 34
5.3 Bifurcation diagrams . . . 35
5.4 Discussion and conclusions . . . 38
6 Photometric lightcurves 41 6.1 Datasets . . . 41 6.1.1 Observational . . . 41 6.1.2 Simulated . . . 42 6.2 Results . . . 43 6.2.1 Observational data . . . 43 6.2.2 Simulated data . . . 44
6.3 Discussion and conclusions . . . 45
7 Summary, conclusions and outlook 49
References 51
Introduction
1.1
Overview
Saturn’s seventh moon, Hyperion (also called Saturn VII), orbiting around its host planet with a period of 21.28 d, was discovered in the XIX century by Bond (1848) and Lassel (1848), but it took more than a century to obtain its quality images due to Voyager 2 (Smith et al. 1982) and Cassini (Thomas 2010) missions. It then became apparent that it is the biggest known highly aspherical celestial body in the Solar System, with an elongated shape and dimensions 360 × 266 × 205 km (Thomas 2010). Since the rotational state of Hyperion was predicted, based on the spin–orbit coupling theory (Goldreich & Peale 1966), to remain in the chaotic zone (Wisdom, Peale & Mignard 1984) due to its high oblateness ω2 and relatively high eccentricity e (see Table 1.1), further analyses and observations, regarding Hyperion in particular as well as other oblate Solar System satellites in general, were conducted on a regular basis.
Table 1.1 Physical parameters of the Saturn–Titan–Hyperion system. (Based on Tarnopolski 2017a; see Table 1 therein.)
Parameter Symbol Value Reference
Saturn’s mass M 5.68 · 1026kg Jacobson et al. (2006)
Titan’s mass m1 1.35 · 1023kg Jacobson et al. (2006)
m1/M 2.4 · 10−4
Hyperion’s major semi-axis a 1 429 600 km Seidelmann et al. (2007); Thomas et al. (2007)
Titan’s major semi-axis a0 1 221 865 km http://ssd.jpl.nasa.gov/?sat_elem
a0/a 0.855
Hyperion’s oblateness ω2 0.79 Wisdom, Peale & Mignard (1984)
Hyperion’s eccentricity e 0.1 Wisdom, Peale & Mignard (1984)
Titan’s eccentricity e0 0.0288 http://ssd.jpl.nasa.gov/?sat_elem
Hyperion’s orbital period T 21.28 d Thomas et al. (2007)
Theoretical analyses and numerical treatment of the rotational dynamics of a generic oblate satellite (including Hyperion in particular) in various settings have been performed widely. After the seminal paper of Wisdom, Peale & Mignard (1984), Klavetter (1989a,b) performed ground-based photometric lightcurve observations, analysed them in the context of specifying Hyperion’s rotational state, and concluded it is not consistent with any regular
4 CHAPTER 1. INTRODUCTION motion. Boyd et al. (1994) applied the method of close returns to a sparse and short-term simulated observations of Hyperion to conclude that a 2.6 year-long data set suffices to infer the temporary dynamical state as regular or chaotic. Black et al. (1995) performed numerical experiments using the full set of Euler equations to model long-term dynamical evolution, and found that the chaotic tumbling leads to transitions between temporarily regular and chaotic rotational states 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, Pivovarov & Starostin (1996) considered a number of models, including the gravitational, magnetic and tidal moments as well as rotation in the gravitational field of two centers, and analysed the structure of the respective phase spaces with Poincaré surfaces of section. A model with a tidal torque was examined analytically by Khan, Sharma & Saha (1998) using Melnikov’s integrals to assess nonintegrability, and assymptotic methods to study the resonant cases. The stability of resonances with application to the Solar System satellites was inferred based on a series expansion of the terms in the equation of rotational motion (Celletti & Chierchia 1998, 2000; see also Celletti 1990). The Lyapunov exponents and spectra were exhaustively examined for a number of satellites (Shevchenko 2002; Shevchenko & Kouprianov 2002; Kouprianov & Shevchenko 2003, 2005), and in case of Hyperion yielded Lyapunov times spanning from 1.5 to 7 times the orbital period. An interesting possibility that Enceladus might be locked in a 1:3 librational–orbital secondary resonance was investigated using the model of planar rotation in the Hamiltonian formalism (Wisdom 2004). A model of an oblate satellite with a dissipative tidal torque was used to examine the basins of attraction in case of low eccentricities, especially with application to the Moon (Celletti & Chierchia 2008). The dynamical stability was examined for many of the known satellites by Melnikov & Shevchenko (2010), and the synchronous spin–orbit resonance in case of Hyperion was found to be unstable. The dynamical modeling using the full set of Euler equations was conducted again by Harbison, Thomas & Nicholson (2011), who also analysed the moments of inertia in light of the precessional period. Melnikov (2014) examined the conditions for emergence of strange attractors in the vicinity of the synchronous resonances when dissipation is introduced for a number of Solar System satellites. Tarnopolski (2015) argued that in order to extract a maximal Lyapunov exponent from the photometric lightcurves of Hyperion, at least one year of dense data is required, but a three year-long data set would be desirable. These findings are with agreement with the results of Boyd et al. (1994), who used the method of close returns. The Hamiltonian formalism was employed also in the research of secondary resonances (Gkolias et al. 2016), and to take non-conservative forces into account (Gkolias et al. 2017). Hyperion’s chaotic rotation, while unusual, is not exceptional in the Solar System: Pluto’s moons are also in a chaotic rotational state (Showalter & Hamilton 2015).
The orbital dynamics of Hyperion in the three-body Saturn–Titan–Hyperion system has been exhaustively examined due to the interesting 4:3 mean motion resonance between Hy-perion and Titan (Peale 1976; Taylor 1992; Stellmacher 1999; Rein et al. 2012). While the impact of Titan’s gravitation on Hyperion’s orbit has been established (Taylor 1987) and the stability of the 4:3 resonance has been considered in great detail (Colombo, Franklin & Shapiro 1974; Bevilacqua et al. 1980), introduction of the gravitational impact of a
sec-ondary body on the rotation of an oblate satellite has been done before for nearly spherical bodies such as Venus (Beletskii & Levin 1981) or low-eccentric orbits in the Pluto–Charon system (Correia et al. 2015). Numerical integrations were performed within the chaotic zone of the Saturn–Titan–Hyperion system (Tarnopolski 2017a and Chapter 5 herein) with pa-rameters ω2 and e such that the perturbation techniques are not valid (Maciejewski 1995).
The approach is general enough to be applicable to moons other than Hyperion.
Hyperion’s long-term observations were carried out twice in the post-Voyager 2 era. In 1987, Klavetter (1989a,b) performed photometric R band observations over a timespan of more than 50 days, resulting in 38 high-quality data points. In 1999 and 2000, Devyatkin et al. (2002) conducted C (integral), B, V and R band observations. The objective of both analyses was to determine whether Hyperion’s rotation is chaotic at the time and to fit a solution of the equation of motion (i.e. Euler equations) to the observations. To the best of the author’s knowledge (Melnikov, priv. comm.) there were no other long-term observations that resulted in a lightcurve allowing the delong-termination of Hyperion’s rotational state (see also Strugnell & Taylor 1990 and Dourneau 1993 for a list of earlier observations). The period determination in (Klavetter 1989a,b) was performed using the phase dispersion minimization (Stellingwerf 1978) and in (Devyatkin et al. 2002) the Deeming method (Deeming 1975) was used. Both authors noted that the periods found (6.6 d and 13.8 d in Klavetter 1989a,b, and 10.8 d in Devyatkin et al. 2002) were statistically insignificant, leading to the conclusion that Hyperion’s rotation is chaotic (or at least not periodic with sufficiently low periods). Although, shortly after the Cassini 2005 passage, a ground-based BVR photometry was conducted (Hicks, Buratti & Basilier 2008), resulting in 6 nights of measurements (and additional 3 nights of R photometry alone) over a month-long period. Unfortunately, this data was greatly undersampled and period fitting procedures yielded several plausible solutions.
In dynamical systems theory, a chaotic behaviour is typically recognised through a pos-itive maximal Lyapunov exponent (mLE), which describes the rate of divergence (or con-vergence in the negative case) of initially nearby phase space trajectories. The Lyapunov spectrum is relatively easy to calculate in the case when the differential equations are known (Benettin et al. 1980; Wolf et al. 1985; Baker & Gollub 1996; Ott 2002). On the other hand, there exist algorithms allowing to obtain an mLE from an experimental or observational time series (Wolf et al. 1985; Kantz 1994), although they are to be used with carefullness, using at least a few hundred data points (Rosenstein et al. 1993; Katsev & L’Hereux 2003) for nonlinear analysis. It is a hard task in astronomy to obtain long-term, well-sampled lightcurves. Although, despite this difficulty, it has been efficiently shown that pulsar spin– down rates exhibit chaotic dynamics (Seymour & Lorimer 2013): by re-sampling the original measurements, artificial time series were produced, equivalent to the original ones, contain-ing such a number of data points that the calculation of the correlation dimension of the attractor, reconstructed via Takens time delay embedding method (Takens 1981), and the mLE estimation, were possible. A similar approach was employed in (Tarnopolski 2015) in order to estimate the mLE based on existing observations and comparison of the numerical results present in the literature with the interpretation of their authors. Conditions that future observations should meet in order to reliably estimate the mLE were also discussed
6 CHAPTER 1. INTRODUCTION (see Chapter 6 herein). Such nonlinear analyses have become common in many fields, e.g. in case of the recently examined X-ray light curve of a microquasar GRS 1915+105 (Mannattil, Gupta & Chakraborty 2016). Contrary to the finding of Seymour & Lorimer (2013) regard-ing pulsars, Mannattil, Gupta & Chakraborty (2016) found no evidence of chaos for the examined microquasar. In particular, it was stressed in both (Tarnopolski 2015) and (Man-nattil, Gupta & Chakraborty 2016) that the stationarity of the time series (and obviously its length) is crucial for the nonlinear analysis of observational astronomical data.
Well-defined methods to reduce chaos in physical systems have been known for a long time now (e.g., Pyragas 1992; Ott 2002). In general, chaos control scheme already demonstrated its usefulness in an astrodynamical application, e.g. in maneuvering the ISEE-3/ICE-3 satel-lite 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 conducted by means of a Hamiltonian chaos control method (Khan & Shahzad 2008) whose performance turned out to be successful in suppressing chaotic rotation. While still rather futuristic, an ability to construct efficient control terms might prove to be important in future asteroid capture mis-sions and mining attempts (Kargel 1994; Sonter 1997; Koon et al. 2000; Levasseur-Regourd, Hadamcik & Lasue 2006; Elvis et al. 2011). It was already shown that a chaotic phase, lasting several thousands of years, is common in the fission process of near-Earth asteroids before they reach their final dynamical state (Jacobson & Scheeres 2011). Chaotic rotation might be prevalent in binary asteroid systems (Jafari Nadoushan & Assadian 2015).
Herein, a construction of a control term that reduces chaos substantially, down to a strictly periodic rotation, is presented. By transforming the equations of planar rotation into a form where the true anomaly is the independent variable [i.e., into the Beletskii (1966) equation], the Hamiltonian setting of the problem was exploited in investigating the prospects of chaos control within the framework of (Ciraolo et al. 2004) in (Tarnopolski 2017b) and Chapter 4 herein. 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 illustrated also with parameters typical for Solar System asteroids.
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 in order to excavate a crater to allow investigation of its interior structure. The crater was not measured directly after the impact due to a large cloud of dust that obscured the view of the orbitting spacecraft. Hence one of the objectives of Stardust-NExT was to image the site of the Deep Impact impactor’s collision with the object. The comet’s rotational state 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 Tempel 1 is not rotating chaotically like Hyperion, it is a plain illustration of the importance
of the rotational state of small Solar System bodies, among which chatic rotation is expected to be common (Jacobson & Scheeres 2011; Jafari Nadoushan & Assadian 2015).
1.2
Aims and objectives of the thesis
1.2.1
Aims
The purpose of the research presented in this thesis is to investigate the chaotic rotational dynamics of an oblate moon, Hyperion exemplary. In general, it aims at filling some gaps in knowledge about the planar rotation of an oblate moon that have not been adequately adressed in the literature so far. In particular, the purpose is to link potential astronomical observations with qualitative descriptors of chaos in a meaningful way, and to examine the impact of perturbations on the rotational dynamics as well as prospects of chaos control.
1.2.2
Objectives
1. To construct, within a Hamiltonian formalism and using the Beletskii equation, a con-trol term that will drastically diminish the diffusion of a chaotic trajectory in the phase space.
2. To establish the impact of the gravitational interaction with another satellite, i.e. whether the additional body suppresses or amplifies chaos.
3. To provide conditions that ground-based photometric lightcurves should meet in order to extract from them the mLE in a reliable way.
Methods
It is convenient to consider a chaotic physical system in the framework of dynamical systems. Consider hence an autonomous dynamical system in a d-dimensional phase space (Gucken-heimer & Holmes 1983; Tabor 1989; Strogatz 1994; Wainwright & Ellis 1997; Wiggins 2003):
˙
~x = ~F (~x), (2.1)
where ~x = ~x(t), and the overdot denotes differentiation with respect to time t. Note that a nonautonomous system can be turned into an autonomous one by simply setting time as an additional dynamical variable governed by an equation ˙t = 1 (Strogatz 1994).
It follows from the Poincaré–Bendixson theorem (Guckenheimer & Holmes 1983; Licht-enberg & Lieberman 1992; Strogatz 1994; Wainwright & Ellis 1997; Alligood et al. 2000; Lakshmanan & Rajasekar 2003; Wiggins 2003; Greiner 2010) that a two-dimensional phase space is too tight for chaos to emerge. Hence, an autonomous dynamical system from Eq. (2.1) needs to be at least three-dimensional for chaotic solutions to occur.1
If the dynamical system in Eq. (2.1) comes, via Hamilton equations ˙q = ∂H∂p, ˙p = −∂H∂q, from a Hamiltonian H(q, p, t), where p are momenta conjugate to the coordinates q, then the system is called Hamiltonian. Then the Liouville’s theorem (Tabor 1989; Lichtenberg & Lieberman 1992; Alligood et al. 2000; Morbidelli 2002; Ott 2002; Wiggins 2003; Greiner 2010) is true: The volume is conserved under the phase flow ⇔ ∇ · ~F = 0. It follows that a Hamiltonian system cannot have an attractor, as the volume of any cloud of initial conditions does not decay assymptotically to zero. A generalization of the Liouville’s theorem states that a subset of the phase space with initial volume V0 evolves according to the equation
(Lichtenberg & Lieberman 1992; Strogatz 1994; Baker & Gollub 1996; Alligood et al. 2000;
1For completeness, it should be noted that the corollary of the Poincaré–Bendixson theorem is valid for continuous-time systems, which will be investigated hereinafter. In particular, chaos can be present in a two-dimensional phase space if the governing map (a discrete-time dynamical system) is invertible (e.g., Chirikov standard map for K < 1, Chirikov 1979; Lichtenberg & Lieberman 1992; Meiss 1992; Zou et al. 2016), and a one-dimensional phase space is sufficient when the map is non-invertible, e.g. the logistic map (May 1976; Feigenbaum 1979). A delay differential system is infinite dimensional (Lakshmanan & Rajasekar 2003; Kantz 2004; Lakshmanan & Senthilkumar 2011), hence one dynamical variable suffices.
10 CHAPTER 2. METHODS
Figure 2.1 A Poincaré surface
of section is a mapping from
an n-dimensional phase space on an (n − 1)-dimensional hypersuface. The points (red) are recorded each time an orbit (black line) crosses such a surface (blue plane).
Ott 2002; Lakshmanan & Rajasekar 2003; Wiggins 2003): d
dtln V (t) = ∇ · ~F . (2.2)
In particular, a strange attractor, i.e. a set with Lebesgue measure zero, is encountered when V (t) → 0 for t → ∞.2
2.1
Poincaré surface of section
A chaotic solution of the dynamical system in Eq. (2.1) has an orbit in a d-dimensional phase space. Due to Poincaré–Bendixson theorem such an orbit must be embedded in a space that is at least three-dimensional. Since a chaotic orbit can be complicated, its visualisation might be difficult to interpret. On the other hand, a projection of a d-dimensional trajectory on a (d − 1)-dimensional hypersurface often leads to self-intersections, which obscure the insight into the underlying dynamics. An idea that manages to overcome these difficulties is the application of the so-called stroboscopic variables, i.e. to sample the trajectory when it crosses some arbitrary hypersurface. E.g. in the famous Hénon & Heiles (1964) 4-dimensional Hamiltonian model of motion of stars in an axial symmetric galactic potential, the energy E is a constant of motion; among the other three free variables, usually the points for such a projection are recorded each time when a trajectory crosses the x = 0 plane (Tabor 1989; Sussman & Wisdom 2001; Goldstein, Poole & Safko 2001; Lakshmanan & Rajasekar 2003; Bountis & Skokos 2012; Lowenstein 2012). Fig. 2.1 displays a schematic illustration of this idea. Such a mapping from d to (d − 1) dimensions is called the Poincaré surface of section (Tabor 1989; Lichtenberg & Lieberman 1992; Baker & Gollub 1996; Murray & Dermott 1999; Alligood et al. 2000; Sussman & Wisdom 2001; Goldstein, Poole & Safko 2001; Ott 2002; Lakshmanan & Rajasekar 2003; Wiggins 2003; Greiner 2010; Bountis & Skokos 2012; Lowenstein 2012). While the hypersurface may be arbitrary, it is natural to sample the points according to some constant period if the system is driven.
0.10 0.15 0.20 0.25 0.30 0.35 0.40 -0.5 0.0 0.5 1.0 γ x 0.150 0.152 0.154 0.156 0.158 0.160 0.20 0.25 0.30 0.35 0.40 0.45 γ x
Figure 2.2 Left: A bifurcation diagram for the Duffing pendulum, ¨x + 0.2 ˙x − x + x3 = γ cos t with (x0, ˙x0) = (1, 1), against the amplitude of the driving term γ ∈ [0.1, 0.4]. The integration was performed for t ∈ [0, 2000], and the first 1000 time units were skipped over to avoid transient behaviour. Right: A magnification of the red rectangle reveals a period-doubling scheme strikingly similar to that known from the logistic map (May 1976; Feigenbaum 1979).
2.2
Bifurcation diagrams
In dynamical systems theory, a bifurcation occurs when an infinitesimal change of a (nonlin-ear) parameter governing the system leads to a topological change in its behaviour. Generally speaking, at a bifurcation, the stability of equilibria, periodic orbits or other invariant sets is changed. The theory of bifurcation is a vast field (Guckenheimer & Holmes 1983; Tabor 1989; Crawford 1991; Lichtenberg & Lieberman 1992; Strogatz 1994; Baker & Gollub 1996; Alligood et al. 2000; Ott 2002; Lakshmanan & Rajasekar 2003; Wiggins 2003; Kuznetsov 2004; Peitgen, Jürgens & Saupe 2004). A bifurcation diagram is a diagram illustrating how does the limit set (fixed points, regular orbits, chaotic attractors etc.) of a dynamical system
˙x = F (x; a), where a is a nonlinear parameter, change with respect to a.
An exemplary bifurcation diagram for the Duffing pendulum is shown in Fig. 2.2 (see also Baker & Gollub 1996; Lakshmanan & Rajasekar 2003). Note that, in fact, any vertical cross section of a bifurcation diagram is a Poincaré surface of section projected on its vertical axis.
Usually, on the bifurcation diagrams there appear to be some curves running through the plot in the chaotic region. These are the so-called critical curves (Neidinger & Annen 1996; Peitgen, Jürgens & Saupe 2004), defined by x = Fn(x
0; a), which are polynomials of
order n.
2.3
Correlation dimension
A fractal dimension (or, more precisely, the Hausdorff dimension; Hausdorff 1919; Theiler 1990; Ott 2002) is often measured with the correlation dimension, dC (Grassberger &
Pro-caccia 1983; Grassberger 1986; Theiler 1990; Strogatz 1994; Alligood et al. 2000; Ott 2002; Lakshmanan & Rajasekar 2003), which takes into account the local densities of the points
12 CHAPTER 2. METHODS in the examined dataset. For usual 1D, 2D or 3D cases the dC is equal to 1, 2 and 3,
re-spectively. Typically, a fractional correlation dimension is obtained for fractals (Mandelbrot 1983; Strogatz 1994). In practical applications one often needs to estimate the fractal dimen-sion of a set consisting of a finite number of points, obtained from observations or numerical simulations, for which the strict Hausdorff dimension is equal to zero, which is obviously not a satisfying result. The correlation dimension gives an estimate of the fractal dimension in such cases.
The correlation dimension is defined as dC = lim
R→0
ln C(R)
ln R , (2.3)
with the estimate for the correlation function C(R) as
C(R) = 1 N2 N X i=1 N X j=i+1 Θ (R − ||~xi− ~xj||) , (2.4)
where the Heaviside step function Θ adds to C(R) only points ~xi in a distance smaller than R
from ~xj and vice versa. The total number of points is denoted by N , and the usual Euclidean
distance, ||.||, is employed. The correlation function C(R) is a measure of the probability of finding two points within a distance smaller than R. It hence gives an average ratio of points lying within a sphere of radius R. The limit in Eq. (2.3) is attained by using multiple values of R and fitting a straight line to the linear part of the obtained dependency. The correlation dimension is estimated as the slope of the linear regression. Eq. (2.3) comes from the fact that C(R) ∼ RdC. The last relation is justified by the fact that the probability of finding some points in a region of size comparable to the one of the examined set is close to 1, and the probability of finding some points in a ball of radius R and dimension dC is
proportional to the size of the ball.
The computations in this work were performed using the python code from (Tarnopolski 2014), with a slight modification so that ln R, instead of R, is uniformly sampled with a constant step ∆ ln R. The factor 1/N2 (Ott 2002) in Eq. (2.4) is often taken as 1/N (N − 1)
(Theiler 1990), but (i) any reasonable N makes this difference insignificant, and (ii) due to the logarithms in Eq. (2.3) such a multiplicative factor leads only to a vertical shift in the log C(R) vs. log R plot, leaving the slope (which is the estimate of dC) unaffected. This is
why in the code presented in (Tarnopolski 2014) this term is neglected altogether.
2.4
Lyapunov exponents
2.4.1
Indicator of chaos
The Lapunov exponent (LE) (Oseledec 1968; Benettin et al. 1980; Wolf et al. 1985; Tabor 1989; Lichtenberg & Lieberman 1992; Strogatz 1994; Abarbanel 1996; Baker & Gollub 1996; Ott 2002; Lakshmanan & Rajasekar 2003; Skokos 2010; Bountis & Skokos 2012; Lowenstein 2012) is a measure of the mean exponential divergence (or convergence) of two initially nearby
orbits of a dynamical system in its phase space in a time limit of infinity. In a d-dimensional system, there are d such LEs, denoted λi. Having an infinitesimal displacement between
two initially nearby orbits, ~w ∈ Rd, their exponential divergence in the i-th direction of the tangent space is characterised by λi > 0, related to chaotic behaviour (in fact it implies
local instability and sensitivity to initial conditions). Let us denote the maximal LE (mLE) by λ1, and order the LEs so that λ1 ≥ λ2 ≥ . . . ≥ λd. The condition λ1 > 0 is sufficient
to detect chaos. If the divergence is slower than exponential, or the trajectories converge, then λ1 ≤ 0, and the orbit is regular. For Hamiltonian systems,
Pd
i=1λi = 0, hence the LEs
come in pairs such that λj + λd−j+1 = 0 (Tabor 1989; Lakshmanan & Rajasekar 2003), and
at least one LE (corresponding to the direction of the flow; usually a time-like phase space direction) is equal to zero.
2.4.2
Lyapunov spectrum
Consider a dynamical system from Eq. (2.1), an initial condition ~x0, and two initially nearby
orbits (one characterised by ~x0) displaced by an infinitesimal deviation vector ~wt, evolving
in time according to the variational equation (Ott 2002; Lakshmanan & Rajasekar 2003; Bountis & Skokos 2012)
˙ ~
wt= J(t) · ~wt (2.5)
(coming from the linearization of the dynamical system, i.e. ~wtbelongs to the tangent space),
where the subscript denotes time dependence, and J(t) is the Jacobian matrix of ~F at time t. Then, the LEs are given by
λi = lim t→∞ 1 t ln || ~wti|| || ~wi 0|| . (2.6)
When the limit t → ∞ is replaced with t sufficiently large, it leads to the so-called finite time LEs (FTLEs).
In practical applications, the variational equation is considered along with the dynamical system itself (Alligood et al. 2000; Ott 2002; Lakshmanan & Rajasekar 2003):
˙
Φt(~x0) = D ~Ft(~x0) · Φt(~x0), (2.7)
where Φt(~x0), a d × d matrix, is the derivative of the flow with respect to the initial condition
~
x0, and D ~Ft(~x0) is the Jacobian evaluated at ~x0. The initial condition for Φ is Φt0(~x0) = Id (identity matrix). Because a chaotic trajectory is bounded in the phase space, the evolution of the flow cannot be performed for an arbitrarily long time; therefore, after a short interval ∆T the vectors forming the matrix Φ are being subject to the Gram–Schmidt procedure to ensure they remain orthogonal and to not allow their norms to become too big, which would lead to obscuring the mean exponential divergence (Benettin et al. 1980; Lakshmanan & Rajasekar 2003). Eventually, according to Eq. (2.6), the norms of such vectors, divided by k∆T (k being the current number of iterations, i.e. the total number of Gram–Schmidt orthogonalizations that took place so far) form the Lyapunov spectrum. Storing the values of LEs (actually, FTLEs) at each iteration allows to examine the convergence of the proce-dure, driven by the structure of the system under consideration. In general, a regular orbit
14 CHAPTER 2. METHODS characterised by a non-positive mLE exhibits a 1/k decay, i.e. a linear convergence to zero with slope −1 in a log–log plot, while a chaotic orbit pleateaus to a finite value. In general, the relation ∇ · ~F =P
i
λi holds (Baker & Gollub 1996; Lakshmanan & Rajasekar 2003).
2.5
Time series analysis
In a way, when the differential equations governing the system’s dynamics are known, virtu-ally all information about the system is accessible. Whether the access is easy and exhaustive is a matter of mathematical tools and computational capabilities. On the other hand, often one has access only to a scalar time series coming from experiments or observations. Such a quantity might come from a 3- or 100-dimensional system; it can be one of the explicit variables of a system, or be a nonlinear algebraic superposition of all 100 variables. In either case, such a time series still contains information about the underlying dynamics. Fortu-nately, the methods of time series analysis (Abarbanel 1996; Kantz 2004; Bradley & Kantz 2015) usually allow to quantify the chaotic dynamical properties of an underlying system (Seymour & Lorimer 2013). The fundament of the numerical analyses performed hereinafter was the TISEAN package3 (Hegger, Kantz & Schreiber 1999) supported with programmes
developed in (Kodba, Perc & Marhl 2005; Perc 2005a,b, 2006). In the following subsections a brief overview of some well establised core techniques is given.
2.5.1
Power spectrum
While the visual inspection of a purely periodic time series might suffice to characterise the dynamics, matters are not that abvious when noise is introduced. Moreover, often chaotic time series resemble noisy ones. For a qualitative discrimation between the three types of motion (i.e., periodic, chaotic and random) a power spectrum P (ν) is an informative tool (Swinney 1983; Tabor 1989; Abarbanel 1996; Baker & Gollub 1996; Lakshmanan & Rajasekar 2003); the differences are most pronounced on a semi-log plot. The P (ν) of a periodic data set consists of clearly distinguishable peaks corresponding to the frequencies ν of the harmonics. A chaotic P (ν) is a broad band spectrum, with all ν present, while lower frequencies have greater power than higher ones. Finally, in a random time series also all ν are present, but the spectrum is flat. Fig. 2.3 displays examples of regular and chaotic solutions for the Duffing pendulum, and a white noise for comparison. The presented spectra are Lomb–Scargle periodograms (Lomb 1976; Scargle 1982).
2.5.2
Phase space reconstruction
A basic technique to gain insight into the dynamics underlying an observed time series is to reconstruct the phase space trajectories via Takens time delay embedding method (Packard et al. 1980; Takens 1981; Eckmann & Ruelle 1985; Lichtenberg & Lieberman 1992; Strogatz 1994; Abarbanel 1996; Baker & Gollub 1996; Alligood et al. 2000; Lakshmanan & Rajasekar
3
0.0 0.5 1.0 1.5 2.0 2.5 3.0 10-6 10-4 10-2 100 102 ν P(ν) (a) 0.0 0.5 1.0 1.5 2.0 2.5 3.0 10-5 10-4 10-3 10-2 10-1 100 101 ν P(ν) (b) 0.0 0.5 1.0 1.5 2.0 2.5 3.0 10-3 10-2 10-1 100 101 ν P(ν) (c)
Figure 2.3 A Lomb–Scargle periodogram for the Duffing pendulum: (a) regular solution with prominent peaks, (b) a broad band spectrum for a chaotic solution. (c) A flat power spectrum of a white noise.
2003; Kantz 2004; Bradley & Kantz 2015). Having a series of scalar measurements x(t) uniformly distributed at times t one can form an m-dimensional location vector P (t) using only the values of x(t) at different times given by the formula
P (t) = x(t), x(t + τ ), x(t + 2τ ), . . . , x(t + (m − 1)τ ), (2.8) where m is the embedding dimension, and τ is the time delay chosen so that the components of P (t) are independent or uncorrelated. These parameters are estimated via the false nearest neighbour (FNN) algorithm (for the embedding dimension m) and the mutual information (MI) method (for the time delay τ ), both briefly described below. For a more detailed description of the FNN and MI algorithms and their implementations, see (Kodba, Perc & Marhl 2005; Perc 2005a,b, 2006; Hegger, Kantz & Schreiber 1999; Kantz 2004; Bradley & Kantz 2015). The programmes4 presented in (Kodba, Perc & Marhl 2005; Perc 2005a,b,
2006) are employed herein. The TISEAN software package (Hegger, Kantz & Schreiber 1999) is also widely used throughout this work.5
False nearest neighbours
The FNN fraction (Kennel, Brown & Abarbanel 1992; Abarbanel 1996; Hegger, Kantz & Schreiber 1999; Kodba, Perc & Marhl 2005; Perc 2005a) is calculated under the assumption that the trajectory folds and unfolds smoothly. Roughly speaking, when two initially nearby points diverge under forward iteration (typically not longer than for a time τ ; see next paragraph below), they do so not faster than Rth, where is the initial separation (not
greater than the standard deviation of the data), and Rth is a certain threshold. When this
behaviour is changed with the rise of m, the points are marked as FNNs; the FNN fraction
4
These programmes are available from the website http://www.matjazperc.com/ejp/time.html. 5The programmes from (Kodba, Perc & Marhl 2005; Perc 2005a,b, 2006) were used as they contain an implementation of the Wolf et al. (1985) algorithm absent in the TISEAN package, and the stationarity test; the program mutual gives the time delay τ estimated both via the MI and the autocorrelation function in one run, what simplifies comparisons.
16 CHAPTER 2. METHODS decreases with m and reaches a value significantly close to zero (in practical implementations to the threshold Rth) for the embedding dimension m considered to be a correct estimate of
the proper one. This is achieved by calculating the FNN fraction as Ri =
|xi+mτ − xt+mτ|
||P (i) − P (t)|| (2.9)
for nearby points P (i), P (t) such that ||P (i) − P (t)|| < .
Mutual information
The MI (Fraser & Swinney 1986; Kantz 1994; Abarbanel 1996; Hegger, Kantz & Schreiber 1999; Kantz 2004; Kodba, Perc & Marhl 2005; Perc 2005a) is an alternative for the also commonly used autocorrelation function (Grassberger & Procaccia 1983; Hegger, Kantz & Schreiber 1999; Lakshmanan & Rajasekar 2003; Perc 2006) using which the time delay τ used to be estimated as the delay at which autocorrelation drops to 1/e (to ensure e-folding), or, rarely, to zero, in order to form uncorrelated location vectors. The MI, on the other hand, gives the amount of information one can obtain about xt+τ given xt. In short, the absolute
difference between xmax and xmin is binned into j bins, j being large enough (the algorithm
is quite robust with respect to the size of the bins), and for each bin the MI as a function of τ , denoted by M I(τ ), is constructed from the probabilities that the variable lies in the h-th and k-th bins (ph and pkrespectively, and h, k = 1, . . . , j) and the joint probability ph,k that
xt and xt+τ are in bins h-th and k-th, respectively:
M I(τ ) = − j X h=1 j X k=1 ph,k(τ ) ln ph,k(τ ) phpk . (2.10)
The most proper τ is estimated as the argument for which the first minimum of M I(τ ) occurs.
2.5.3
Maximal Lyapunov exponents
The Wolf et al. (1985) method
The algorithm finds a nearest neighbour to an initial point and evolves them both until the separation becomes too big; next, the distance is being rescaled in order to stay in the small-scale structure, and this scheme repeats to the end of the time series. Then, the average of the logarithms of the displacement ratios is the estimate of the mLE (Kodba, Perc & Marhl 2005; Perc 2005a,b, 2006). The main drawbacks of this method are that (i) it fails to take advantage of all available data as it focuses on one fiducial trajectory (Rosenstein et al. 1993), and (ii) it does not test the presence of exponential divergence (a behaviour underlying chaotic dynamics) but assumes it explicitly ad hoc, what may lead to spurious results (Kantz 2004).
The Kantz (1994) method
This method differs from the previous one in that it takes several points in a neighbourhood of some particular point xi. Next, one computes the average distance of all obtained trajectories
to the reference, i-th one, as a dependence on the relative time n (incorporated to the k-th subscript as follows: xk+(m−1)τ +n). The average S(n) of logarithms of these distances (i.e.,
the stretching factors) is plotted as a function of n and the slope of the linear part is the mLE (see Kantz 1994; Hegger, Kantz & Schreiber 1999; Kantz 2004; Perc 2005a for further details). In the case of chaos, three regions should be distinct: a steep increase for small n, a linear part, and a plateau (Seymour & Lorimer 2013).
However, the methods described require the data to fulfill the stationarity assumption, i.e. that the statistical properties of a time series (e.g., mean and standard deviation) are constant in time. One therefore needs to perform a test to ensure these methods are applicable. For such a purpose the stationarity test is employed as follows.
2.5.4
Stationarity test
Let us consider a point P (t) as an event and find all similar events, i.e. those points P (i) that lie not further than from P (t). One averages all values of xi and call this a prediction
of a future observation based on the value of xt. The key now is to use cross-prediction, i.e.
to partition the whole dataset into non-overlapping segments and use the j-th segment to make a prediction of a k-th segment. One quantifies the correctness by calculating the error, δjk, via square root of mean square deviations from the mean in segment k and repeats this
for all j and k:
δjk = v u u t 1 N N X k=1 (˜xk− xk)2, (2.11)
where ˜xk is the prediction and xk is the true value in a k-th segment. If δjk is significantly
larger than the average, then either the dynamics is not conserved from one segment to another, or the data is undersampled. Both cases yield a conclusion that the data is non-stationary (Schreiber 1997; Kantz 2004; Perc 2005b, 2006).
Rotational model of an oblate moon
3.1
Planar rotation
The equation of rotational motion can be derived based on the following assumptions (Danby 1962; Wisdom, Peale & Mignard 1984; Murray & Dermott 1999; Sussman & Wisdom 2001; Greiner 2010):
1. the orbit of the satellite around the planet is Keplerian with eccentricity e and major semi-axis a:
r = a (1 − e
2)
1 + e cos f, (3.1)
where f is the true anomaly given by the differential equation ˙ f = √ GM [a (1 − e2)]3/2 (1 + e cos f ) 2 , (3.2)
with M being the mass of the planet, and G the gravitational constant;
2. in general, the physical model of the satellite is a triaxial ellipsoid; however, to simplify calculations, the satellite is simulated by a double dumbbell with four mass points 1 to 4 (see Fig. 3.1 for the geometrical setting) with equal mass m each, arranged in the orbital plane. The principal moments of inertia are A < B < C;
3. the satellite’s spin axis is fixed and perpendicular to the orbit plane; the spin axis is aligned with the shortest physical axis, i.e. the one corresponding to the largest principal moment of inertia.
The parameters from Table 1.1, corresponding of the Saturn–Titan–Hyperion system, are used unless stated otherwise.
Defining the oblateness as ω2 = 3(B−A)
C , and choosing the units so that the orbital period
T is equal to 2π and the major semi-axis a = 1 (which implies through Kepler’s third law GM = 1, and that the orbital mean motion n = 1), eventually the equation of motion takes the form ¨ θ + ω 2 2r3 sin 2 (θ − f ) = 0, (3.3) 19
20 CHAPTER 3. ROTATIONAL MODEL OF AN OBLATE MOON H f θ α O α= θ - f S r
Figure 3.1 Rotational model of an oblate moon. (Based on Tarnopolski 2017b; see Fig. 1 therein.)
where time is dimensionless and ˙θ is measured in units of n, with Eq. (3.1) for the orbit in the form
r = (1 − e
2)
1 + e cos f, (3.4)
and Eq. (3.2) for the true anomaly yields ˙
f = 1
(1 − e2)3/2 (1 + e cos f ) 2
. (3.5)
It will be useful to turn Eq. (3.3) and (3.5) into a dynamical system as in Eq. (2.1): ˙ θ = ϕ, (3.6a) ˙ ϕ = ω 2 2r3 sin 2(f − θ), (3.6b) ˙ f = (1 + e cos f ) 2 (1 − e2)3/2 , (3.6c)
where r is given by Eq. (3.4). Eq. (3.6a) and (3.6b) carry dynamical information, and Eq. (3.6c) is an auxiliary equation to resolve the time dependence.
Using the well-known general relations in the two-body motion, µ = GM , h = r2f ,˙ h2 = µl, l/r = 1 + e cos f [i.e., Eq. (3.1)], l = a(1 − e2) (Khan, Sharma & Saha 1998), one
can transform Eq. (3.3) so that f is the independent variable [the famous Beletskii (1966) equation1]: (1 + e cos f )d 2 θ df2 − 2e sin f dθ df − ω2 2 sin 2(f − θ) = 0, (3.7) which as a (nonautonomous) dynamical system reads
θ0 = η, (3.8a)
η0 = 2e sin f · η −
ω2
2 sin 2(θ − f )
1 + e cos f , (3.8b)
where the prime denotes differentiation with respect to f .
1which is an equation for δ = 2α = 2(θ − f ); herein the version in Eq. (3.7), i.e. an equation for θ, is used.
0 π/4 π/2 3π/4 π -0.5 0. 0.5 1. 1.5 2. 2.5 3. θ θ (a) 0 π/4 π/2 3π/4 π -0.5 0. 0.5 1. 1.5 2. 2.5 3. θ (b) 0 π/4 π/2 3π/4 π -0.5 0. 0.5 1. 1.5 2. 2.5 3. θ (c)
Figure 3.2 Phase space in stroboscopic variables for (a) e = 0.1, ω2= 0.79 (i.e., Hyperion’s parameters), (b) e = 0.005, ω2 = 0.79, (c) e = 0.1, ω2 = 0.04. Different colours correspond to different trajectories. Eq. (3.3) was integrated for 2,000 dimensionless time units, resulting in 319 points for each trajectory. The initial conditions were (0, ˙θ0) and (π/2, ˙θ0), where ˙θ0 ∈ {0.0, 0.1, . . . , 2.5, 2.6}. (Figure after Tarnopolski 2017a; see Fig. 2 therein.)
The phase space of Eq. (3.6) is 3-dimensional: Ω = {(θ, ˙θ, f ) : θ ∈ R mod 2π, ˙θ ∈ R, f ∈ R mod 2π}. But f is a regular, 2π-periodic function and does not carry much information (i.e., it is a time-like direction in the phase space). Moreover, in dimensionless units the orbital period of the oblate moon is also equal to 2π. Hence, a Poincaré surface of section, constructed by taking the values of (θ, ˙θ) with a time step of 2π provides insight into the rotational dynamics. Furthermore, the rotation of the satellite by 180◦ (i.e., θ → θ + π) gives an equivalent physical configuration, hence θ can be confined to the interval [0, π). Such exemplary surfaces of section are shown in Fig. 3.2.
As is common in Hamiltonian systems, the phase space is divided into regions occupied with chaotic trajectories, and regions of regular (periodic or quasiperiodic) motion. There are also motions called sticky orbits, when the trajectory initially behaves in a regular manner and diverges into the chaotic zone after some time. The phase space in Fig. 3.2(a) is dominated by a large chaotic sea, formed by merging the chaotic zones surrounding spin–orbit resonances from p = 1 : 2 to p = 2 : 1 when ω2 increases (Wisdom, Peale & Mignard 1984). Quasiperiodic motions are represented by smooth curves and by closed curves, e.g. the ones surrounding the synchronous p = 1 : 1 resonance (see Table 1 in Black et al. 1995 for locations of the surviving resonances). When the initial condition is located near the boundary between regular motion and the chaotic sea, sticky motion occurs. A narrow chaotic zone is present also in Fig. 3.2 (b) and (c) obtained for smaller values of the eccentricity and oblateness, respectively. In fact, every Solar System satellite has a chaotic period in its past (Spohn, Johnson & Breuer 2014). The critical value of ω for which wide-spread chaos is expected to occur, based on the Chirikov (1979) resonance overlap criterion, is ωRO = 1/(2 +√14e),
which is approximately 0.31 for Hyperion’s eccentricity (Wisdom, Peale & Mignard 1984), and the actual ω exceeds ωRO significantly (see Table 1.1).
22 CHAPTER 3. ROTATIONAL MODEL OF AN OBLATE MOON
3.2
Three-dimensional rotation
3.2.1
Euler equations
To realistically model the rotational dynamics of Hyperion, one needs to consider a three-dimensional rotation in terms of a full set of Euler equations:
A ˙ωa− ωbωc(B − C) = − 3 r3βγ(B − C), (3.9a) B ˙ωb− ωcωa(C − A) = − 3 r3γζ(C − A), (3.9b) C ˙ωc− ωaωb(A − B) = − 3 r3ζβ(A − B), (3.9c)
where ωi are the components of the angular velocity vector ~ω, and ζ, β, γ are the direction
cosines. The modeling, while straightforward on one hand, consists of several steps which make the procedure rather complex. A detailed description can be found in (Klavetter 1989b; Devyatkin et al. 2002); see also (Wisdom, Peale & Mignard 1984; Klavetter 1989a; Black et al. 1995; Shevchenko & Kouprianov 2002; Kouprianov & Shevchenko 2003, 2005; Harbison, Thomas & Nicholson 2011). Herein, only a brief overview is given.
The components of the angular velocity ~ω are
ωa= ˙θ sin ϕ sin ψ + ˙ϕ cos ψ, (3.10a)
ωb = ˙θ sin ϕ cos ψ − ˙ϕ sin ψ, (3.10b)
ωc= ˙θ cos ϕ + ˙ψ, (3.10c)
and the direction cosines are
ζ = cos(θ − f ) cos ψ − sin(θ − f ) cos ϕ sin ψ, (3.11a) β = − cos(θ − f ) sin ψ − sin(θ − f ) cos ϕ cos ψ, (3.11b)
γ = sin(θ − f ) sin ϕ. (3.11c)
After differentiating Eq. (3.10) and substituting into Eq. (3.9) one obtains ¨ θ = k1sin ψ + k2cos ψ sin ϕ , (3.12a) ¨ ϕ = k1cos ψ − k2sin ψ, (3.12b) ¨ ψ = k3 − ¨θ cos ϕ, (3.12c) where k1 = ωbωc− 3 r3βγ B − C
A − ˙θ ˙ϕ cos ϕ sin ψ − ˙θ ˙ψ sin ϕ cos ψ + ˙ϕ ˙ψ sin ψ, (3.13a) k2 = ωcωa− 3 r3γζ C − A
B − ˙θ ˙ϕ cos ϕ cos ψ + ˙θ ˙ψ sin ϕ cos ψ + ˙ϕ ˙ψ cos ψ, (3.13b) k3 = ωaωb− 3 r3ζβ A − B C + ˙θ ˙ϕ sin ϕ. (3.13c)
The Euler coordinates have a singularity at sin ϕ = 0 due to the denominator of Eq. (3.12a). For this reason, whenever the solution becomes nearly singular, one needs to change the coordinates to another set that is free of a singularity (for this particular state). These new coordinates are termed the Wisdom coordinates (Klavetter 1989b; Devyatkin et al. 2002), and the transformation from the Euler to Wisdom coordinates is as folows:
ψW = arctan − cos ϕ sin ψ
sin ϕ , (3.14a)
ϕW = arctan sin ϕ
cos ψWcos ψ cos ϕ, (3.14b)
θW = arctan cos θ sin ψ + sin θ sin ϕ cos ψ
cos θ sin ϕ cos ψ − sin θ sin ψ, (3.14c) with a similar transformation back to the Euler coordinates, and the corresponding direction cosines. The Wisdom coordinates have a singularity whenever cos ϕW = 0. During the
integration of the equations of motions, the transformation between the two coordinate systems is performed whenever | sin ϕ| < 0.1 or | cos ϕW| < 0.1.
3.2.2
Lightcurve modeling
The orientation of the satellite in three dimensions has to be converted to brightness values as would be seen on Earth. To do so, an instantaneous projection of the ellipsoid on the satellite–planet direction in planetocentric coordinate system (x, y, z) can be written as
Gy2+ 2J yz + Kz2 = 1, (3.15)
where G, J , K are coefficients dependent on θ, ϕ, ψ. Writting Eq. (3.15) in polar coordinates one can express the area A = A(ξ) of the projected ellipse in dependence on the angle of rotation ξ = 12arctanG−K2J describing the maximum distance from the center. This gives the projected area as seen from Saturn. To convert it to an Earth-based projection, rotations by an angle due to the difference between the longitudes of the Sun and Saturn, and by an angle between the orbital plane and the ecliptic need to be performed.
To relate the brightness and time in Julian date, the location of the satellite in its orbit is expressed with the mean longitude, the pericenter longitude, and orbit’s eccentricity. The longitudes give the mean anomaly M that relates with the eccentric anomaly E through the Kepler equation, M = E − e sin E , from which one can then calculate the true anomaly f . To take into account the phase angle κ for the Sun–Saturn–Hyperion system and the reflecting properties of the satellite’s surface, the model of (Lagerkvist & Williams 1987) is employed. The correction to the mean opposition magnitude is done via the formulae
m0 = m − 5 log r∆ rS∆S , (3.16a) m0(κ) = H − 2.5 log [(1 − G)Φ1(κ) + GΦ2(κ)] , (3.16b) Φi(κ) = exp −Ai tanκ 2 Bi , (3.16c)
24 CHAPTER 3. ROTATIONAL MODEL OF AN OBLATE MOON where m—observed magnitude, m0—mean opposition magnitude uncorrected for Solar phase, r—heliocentric distance, ∆—geocentric distance, rS—Saturn’s mean heliocentric distance,
∆S—Saturn’s mean geocentric distance, and H, G, Ai, Bi—suitable parameters (Klavetter
1989a,b; Devyatkin et al. 2002); in particular, H is the reduced magnitude corrected for mean opposition and zero Solar phase κ = 0.
Chaos control
Preface
A model of planar oscillations, described with the Beletskii equation, is investigated. The Hamiltonian formalism is utilised to employ a control method for suppressing chaos. An additive control term—an order of magnitude smaller than the potential—is constructed, and allows not only to diminish significantly the diffusion of the trajectory in the phase space, but turn the purely chaotic motion into strictly periodic.
Publications related to the chapter
• Tarnopolski M., Rotation of an oblate satellite: Chaos control, submitted to A&A, arXiv:1704.02015 (2017)
4.1
Chaos control
The approach of Ciraolo et al. (2004) was employed to construct a control term. The method aims at constructing a term F =
∞
P
s=2
Fsof the order O(ε2) that when added to a Hamiltonian
H = H0+ εV , creating thus a new Hamiltonian ˜H = H + αF = H0+ εV + αF , where α is a
(sufficiently small) scaling factor1, that the motion is much less diffused in the phase space, i.e. chaos is suppressed.
The Hamiltonian2 H = p 2 2 (1 + e cos f )2 + E | {z } H0 −ω 2 4 (1 + e cos f ) cos 2(θ − f ) | {z } V , (4.1)
1Note that herein this α is different than the angle α = θ − f defined in Section 3.1 and used also in Chapter 5.
2E ∈ const is a second action introduced in order to make the Hamiltonian autonomous.
30 CHAPTER 4. CHAOS CONTROL where ε ∼ ω2, gives Eq. (3.7) as the equation of motion via Hamilton equations. The series
forming F was truncated after the first term, and its zeroth-order expansion takes the final form:
F2(0)= −
ω4sin 2(f − θ) {2e [9 sin(f − 2θ) + sin(3f − 2θ)] + 9 sin 2(f − θ)}
288(1 + e cos f ) , (4.2)
which is indeed O(ω4), and the index in brackets emphasizes it is a zeroth-order expansion of F2. Inserting into the new Hamiltonian, ˜H = H0 + εV + αF2(0) gives the equation of
motion in the form
(1 + e cos f )θ00− 2e sin f θ0 = 1
72(1 + e cos f )2
n
9ω2sin 2(f − θ)2e2(1 + cos 2f ) + 8e cos f +
− αω2cos 2(f − θ) + 4 − αeω4[9 sin(3f − 4θ) + sin(5f − 4θ)]o.
(4.3)
Eq. (4.3) was integrated numerically for e = 0.1, ω = 0.89 (i.e., Hyperion’s values), f ∈ [0, 5000], and α ∈ [0, 9.9] with a step ∆α = 0.1. The case α = 0 corresponds to a lack of control term, i.e. the Hamiltonian in Eq. (4.1), hence will serve as a reference. Fig. 4.1 shows some of 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 did not influence the surface of section qualitatively, i.e. the diffusion was not diminished at all. Impressively, α = 2.2 − 2.4 suppressed chaos to a very small region of the phase space, forming closed curves depicting quasiperiodic motion. 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 behaviour is observed for α = 3, where the motion is turned to periodic. When α is increased further, the motion turns back to chaotic, with occasional signs of stickiness, e.g. α = 4.5, 6.4, 8.2 or 9.6.
The maximal value of the potential V is 0.218, while for the control term F2(0) (i.e.,
taking α = 1) it does not exceed 0.022—i.e. the control term is only about 10% of the potential. Note that α & 10 is useless since the control term is then of the same order of magnitude as the potential V . However, for α ∼ 2.2 it is still nearly an order of magnitude smaller than V , and sufficient to suppress the chaotic behaviour of the investigated system almost entirely, i.e. convert it to quasiperiodic rotation.
To examine a more astrodynamically realistic scenario, the above computations were repeated with ω2 = 0.25 and e = 0.005. The eccentricity is the mean value for the Solar
System satellites (Tarnopolski 2017a), and ω2 was decreased slightly (characterizing still a quite oblate object, and greater than the critical value ωRO = 1/(2 + √14e) from the
resonance overlap criterion), compared to Hyperion’s value, as in the method applied it 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 V is about 0.14, i.e. twice as small as when Hyperion’s parameters were employed.
Figure 4.1 Poincaré surfaces of section obtained by integrating Eq. (4.3) and recording the points (θ(f ), θ0(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. (Figure after Tarnopolski 2017b; see Fig. 3 therein.)
4.2
Discussion and conclusions
The range α ∈ [0, 2.2] in the control term in Sect. 4.1 was also sampled with a smaller step of 0.01, but the resulting surfaces of section did not reveal any signs of chaos suppression for α < 2.2. Also, the first-order series expansion of F2, denoted F2(1), was tested; it
allowed to suppress chaos down to a periodic motion with α = 0.7, but its overall maximal value attained 1/3 of the maximal value of the potential V , i.e. there was no improvement comparing to the performance of F2(0).
The results have a two-fold meaning: first, it is remarkable that a chaotic orbit can be turned into a strictly periodic one with such a simple control term, even though its magnitude, an order of magnitude smaller than V , might not be that impressive at first sight. Second, in the light of this result, higher-order series or other methods (e.g. Tsui & Jones 2000) are prospective and promising when ω2 is treated as a perturbing factor.
Titan’s influence
Preface
In this chapter the gravitational influence of Titan on the rotation of Hyperion is examined by evolving compact sets of initial conditions in the phase space and via bifurcation diagrams. The evolution of the former is quantified with the correlation dimension. It is found that the presence of another satellite causes some chaotic trajectories to become regular. Moreover, the highly structured picture revealed by bifurcation diagrams is destroyed when the influence is included.
Publications related to the chapter
• Tarnopolski M., Influence of a second satellite on the rotational dynamics of an oblate moon, Celestial Mechanics and Dynamical Astronomy, 127(2), 121–138 (2017)
5.1
Rotational model
A second satellite—Titan—is assumed to revolve around the planet on a circular orbit, with radius a0, in the same plane as the oblate moon. Based on Eq. (3.5), the true anomaly
depends linearly on time:
fT =
1
a3/20 t, (5.1)
where a0 is Titan’s semi-major axis (see Table 1.1). One obtains the distance between the
two satellites as1 (see Fig. 5.1)
rT H = r r 1 −2a0 r cos (fT − fH) + a0 r 2 . (5.2) 1In this chapter, f
H and fT denote the true anomalies of Hyperion and Titan, respectively.
34 CHAPTER 5. TITAN’S INFLUENCE
H
f
Hθ
α
O
S
r
a
0T
f
T rTHα
1x
1 2 3 4Figure 5.1 Geometry of the model including the orbital motion of a second satellite. S, T and H stand for Saturn, Titan and Hyperion (center of mass), respectively. (Figure after Tarnopolski 2017a; see Fig. 1 therein.)
The angle α1 is also required. One finds2
α1 = x + fT − θ − π (5.3) and x = arccos r 2 T H+ a20− r2 2rT Ha0 . (5.4)
Finally, one obtains the following equation of motion: ¨ θ = ω 2 2 ( sin 2 (fH − θ) r3 − m1/M r3 T H sin 2 arccos a0− r cos (fT − fH) rT H + fT − θ ) , (5.5) where m1/M is the mass ratio of Titan and Saturn. The initial conditions for the true
anomalies will be assumed throughout this chapter to be fH(0) = fT(0) = 0.
5.2
Evolving clouds of initial conditions
To investigate the differences between the absence and presence of the gravitational interac-tion with another satellite on the rotainterac-tion of an oblate moon, the evoluinterac-tion under the flow of two arbitrarily chosen sets of initial conditions is analysed:
1. IC1—a total of 101 × 101 = 10, 201 points distributed uniformly on a 0.1 × 0.1 square centered on (π/2, 0.55);
2. IC2—similar to IC1, but centered on (π/2, 1.5).
IC1 was chosen so that it is located on the edge of the chaotic sea and the domain of quasiperiodic motion. For each point in both sets IC1 and IC2, the equations of motion in Eq. (3.6) and (5.5) were solved numerically for 20 orbital periods and (θ, ˙θ) are recorded at every periapse. The sets to which IC1 and IC2 evolved after k orbital periods are dis-played in Fig. 5.2. The correlation dimension dC after each k revolutions was computed. 2Here, as in (Tarnopolski 2017a), the convention of (Murray & Dermott 1999; Greiner 2010) is employed, i.e. the angles are undirected, and hence, from the triangle OSH, α = fH− θ, and Eq. (5.3) comes from the triangle T SH. See (Murray & Dermott 1999; Greiner 2010; Tarnopolski 2017a) for details.
