Jean-Yves Chemin

Michael Joyce

Michel Capderou

Kevin Baillie

Nicolas Rambaud

Alexis Petit

PhD Auriane Siegfried Egal

Eggl

Anatoliy Maria

Melaine Ivantsov

Irina Saillenfest

Aaron Kovalenko

Jérôme Rosengren

Vishnu Daquin

Raphaël Viswanathan

Giacomo Errani

Lari

Gaël Noirot

Matthieu Chretien

Brice Francois

Evan Gaël

Siwei Eams

Anelise Audibert

Paola Dimauro

Caterina Umilta

Daniele Palaferri

Elena Bellomi

Letizia Cap- Itanio

Salvatore Colombo

Carolina Parroni

Fiorella Polles

Isadora Bicalho

Pedro Hasselmann

Fernando Caro

Adarsh Ranjan

Bogdan Ciambur

Clément Brisa Gianluca

Julien Mancillas

Dassa-Terrier

Rose Corentin

Meriem E L Coogan

Miguel Yajouri

Valeria Marques

Vladan Markov To

Elodie Belleuvre

Ayoub Ladaci

Benny Mandhub

Hakim Bennabi

Ada Kruczkowska

Ansaf Hassani

Bakari Nouhou

Madjid Namane

Etienne Marti

George Farkouh

Manuela Bratti

Daniela Brucalossi

Daria Strelnikova

Samar Dridi

Sofia Kris- Tensson

John Jayne

Lotfi Kouadria

Pierre Moon

Paula Chabrol

Roberta Casil- Las

Shane Lada

Léopold Carron

Maria Del Mar

Charaf Lachachi

Pamela Gutierra

Sofia Bernal

Iliya Kurtev

Sonia Amy

Hortense Amir

Zuzanna Decool

Beatriz Nowicka

Lauren Vigiato

Isabella Kacher

Olesya Kildyushevskaya

Perrin Brouhot

Amélie Thaler

Louisa Baumgartel

Paolo Massimeti

Ibrahim Nouali

Katarina Giglio

Henny Ryan Lucas

Messin

Dédicasse Au

Than-Cha RAZE Jules Rambaud

Audrey Amoussou

Morgane Portal

Wendy Herry

Aliénor Bernardie

Cyril Sol- Biac

Amélie Duplessis

Sarah Bensadoun

Alice Limoli

Giulia Gritsch

Roberta Scattolini

Chiara Palazzi

Pascal Laugier

Sarah Man- Gou

Claire Pascale

Kearvina Sohun

Hortense Niklas

Ar- Ianna Comin

Alexandre Garnier

Tanguy Rodelar To Germinal

Arnaud Gourion

Gerardo To

Lucia Aviles

Sophie Ameijeiras

Anissa Philippe

Marion Kalfa

Carla Bousquet

Tugba Abaci Palominos

Félix Polycarpe

Dorothée Lano

Rachel Polycarpe

William Polycarpe

Next I want to thank many students and comrades that I acknowledge the support

The work I present here is the last and absolutely not the least effort I have done during my studying years. After nine years in academia, studying courses, passing exams, looking for advisors and writing theses, here I stand now as a PhD writing these lines of acknowledgements. Many people asked me if it was worth it. I never asked myself that question and I guess this is why the answer is simply yes. All the work, the knowledge, the self-teaching ... I don't know if I could have done it somewhere else or later in my life. And a satisfying thought comes to my mind when I think of all this : I was doing what I think was worth doing, I was doing what I wanted to do, I was doing what I like. The difficult moments and the pain are just part of a

Introduction

The system of Saturn Saturn is the sixth planet of the Solar System. As it is observable with a naked eye to humans, it has been subject to many considerations throughout History. Saturn takes his name from the ancient Roman god of sowing or seed, also celebrated by the ancient Greeks as Cronus. Saturn, like other planets and celestial bodies, was studied by many civilisations that existed and was known as the most distant planet until the discovery of Uranus in 1789 by William Herschel, thus being a central character in various mythologies. During the early ages, Saturn was subject to many observations and measurements of its motion in the sky by the Babylonians but also later on by Islamic astronomers during the last millennium, when the geocentric model of Ptolemy was still in use.

It is at the beginning of the modern era that ideas were changing. First by the work of several Renaissance scientists : Copernicus, Johannes Kepler, Tycho Brahe and Galileo Galilei, from whom the heliocentric model of the solar system was brought to reality. But also by the development of optical instruments capable of precise observations, such as the refracting telescope. In 1610, using one of his own telescopes Galileo managed to distinguish rings around Saturn, without really being able to characterise them as such. But a few decades later Christian Huygens, using a more powerful telescope was the first to claim that Saturn had rings. This was in 1656 and their composition was still under investigation. Huygens is also known as being the discoverer of Titan, the first body discovered orbiting around Saturn. Later, two other satellites were discovered : In 1671 and 1672, Giovanni Domenico Cassini, director of the new Paris Observatory, observed Iapetus and Rhea. He also managed to spot a gap in the rings, which is known today as the Cassini division. Then he added two more bodies on the list of Saturn's satellite, Tethys and Dione. After a century with no significant discoveries around Saturn, Herschel was able to observe two close moons around Saturn, namely Mimas and Enceladus in 1789. Then in the 19 th century, Hyperion was discovered between the orbits of Titan and Iapetus.

Today there are a confirmed 61 moons orbiting Saturn1 , but the satellites enumerated above are known to be the major ones. Their masses and orbits make them play an important role in the dynamics of the system and it is one of the general aims of this thesis to understand their motion as well as their history.

Space missions and planetology

The 20 th century saw the beginning of space exploration. Many spacecraft were sent all around the solar system and Saturn was the centre of much interest. In 1979, Pioneer 11 made the first flyby of the planet, giving details about the planetary rings. A year later, Voyager 1 confirmed that Titan had a thick atmosphere, completely close to 1:2 meaning that while Mimas orbits twice around Saturn, Tethys accomplishes only one orbit. Therefore, one can understand that after one resonant cycle, both satellites will find themselves in similar positions as before the cycle, therefore the gravitational pull between the two satellites is subjected to cyclic behaviours. Today one can observe two other resonances around Saturn : Enceladus and Dione have their period ratio also close to 1:2 and Titan and Hyperion with a 4:3 commensurability.

More than just a relation between two mean motions, a resonance also involves the variation of the pericentres and the ascending node. For instance, in the case of Titan and Hyperion one would have at first that 3n 6 -4n 7 ≈ 0

(1)

with n 6 and n 7 the mean motion of Titan and Hyperion. But a closer analysis of the dynamics shows that the angle

φ 4:3 = 3λ 6 -4λ 7 + ̟ 7 (2)
oscillates around 180 degrees during their evolution.

FIGURE 2: Evolution of the resonant angle in the case of the Titan-Hyperion 4:3 mean motion. This figure produced from the outputs of numerical simulation involving all the satellites, Jupiter and the Sun as perturbers. The initial conditions were taken from SAT389.

One calls it an eccentricity-type resonance as the element involved the longitude of the pericentre of Hyperion. This resonant behaviour makes the closest approach of those two bodies happening at the apocentre of Hyperion, and Hyperion having an average eccentricity around 0.1, no real close encounters can happen between Titan and Hyperion. This geometry makes this resonance stable. In the case of Mimas and Tethys, there is a special inclination-type resonance, the resonant angle being φ 4:2 = 2λ 1 -4λ 3 + Ω 1 + Ω 3

(3)

Contents

and it involves both ascending nodes. Mimas 6.60 × 10 -8 1.24 × 10 -3 2.18 × 10 -2 1.57 Enceladus 1.90 × 10 -7 1.59 × 10 -3 6.35 × 10 -3 9.89 × 10 -3 Tethys 1.09 × 10 -6 1.97 × 10 -3 9.68 × 10 -4 1.09 Dione 1.93 × 10 -6 2.52 × 10 -3 2.93 × 10 -3 2.90 × 10 -2 Rhea 4.06 × 10 -6 3.52 × 10 -3 8.00 × 10 -4 3.19 × 10 -1 Titan 2.37 × 10 -4 8.17 × 10 -3 2.86 × 10 -2 3.60 × 10 -1 Hyperion 9.79 × 10 -9 9.93 × 10 -3 1.27 × 10 -1 1.02 Iapetus 3.18 × 10 -6 2.38 × 10 -2 2.79 × 10 -2 15.5

Fast tidal migration and new approach to the system

There exists a large number of mean motion resonances in the Solar System. It is understood that two bodies are unlikely to have been created with their orbital period verifying a commensurable relation [START_REF] Roy | On the occurrence of commensurable mean motions in the solar system[END_REF] but instead, mean motion resonances are likely the outcome of changes in semi-major axes due to tidal interaction (Goldreich, 1965a) with the planet.

da dt = 3k 2 R 5 p nm m S a 4 arctan 1 Q (5)
where a the semi-major axis, t the time, k 2 the Love Number, R p the planetary radius, n the mean motion, m the mass of the satellite, m S the mass of Saturn and Q the quality factor. We will elaborate the tidal theory in chapter 1. In the case of Saturn, a first estimate of the tidal migration was given by [START_REF] Goldreich | Q in the Solar System[END_REF] when considering the orbit of satellites. Assuming that these satellites should have been created at the beginning of the Solar System (around 4.5 billion years ago), Goldreich estimates Q to be around 60 × 10 3 for Saturn. Assuming a similar value [START_REF] Allan | Evolution of Mimas-Tethys Commensurability[END_REF] managed to estimate an age for the Mimas-Tethys commensurability of around 240 million years [START_REF] Sinclair | On the origin of the commensurabilities amongst the satellites of Saturn[END_REF] gave numerical probability of around 4% for this resonance to occur. However, these early studies also pointed out that this pair should have encountered the 2λ 1 -4λ 3 +2Ω 1 first with a higher probability of 7% of capture, which could have made the tidal migration scenario collapse. We also mention that the probability of capture is highly dependant on the model considered. Indeed, [START_REF] Champenois | The Role of Secondary Resonances in the Evolution of the Mimas-Tethys System[END_REF] invokes the importance of secondary resonances and the fact that the eccentricity of Tethys is an important parameter for the capture into the ii ′ sub-resonance to happen. Also, the high value of Q could not really account for the apparent advanced resonant state between Enceladus and Dione (low amplitude) and there were some concerns about the possibility of Titan driving Hyperion eccentricity up regarding the low migration rate induced by this low dissipation in Saturn [START_REF] Sinclair | On the origin of the commensurabilities amongst the satellites of Saturn[END_REF].

In contrast to the first statements made on resonant captures, [START_REF] Colombo | On the formation of the orbit-orbit resonance of Titan and Hyperion[END_REF], as well as [START_REF] Greenberg | Evolution of satellite resonances by tidal dissipation[END_REF] argues that the value Q given in [START_REF] Goldreich | Q in the Solar System[END_REF] does not need to apply to Titan. Indeed, such a parameter may depend on the satellite's mean motion, raising tidal bulges with different amplitudes and behaviour. Q and k 2 can be frequency dependant and Colombo, Franklin, Contents and Shapiro, 1974 concludes that Q = 2000 can permit an automatic capture of Titan and Hyperion into their 4:3 mean motion resonance, Hyperion's eccentricity being under 0.02 prior to their capture.

The system of Saturn was the aim of many astrometrical campaigns during the last century. The effect of tides on the orbit of the moons were thought to be beyond the precision of instruments, but from a wide collection of observations spanning 123 years of data and assuming frequency-dependant tidal responses for the moons [START_REF] Lainey | Strong Tidal Dissipation in Saturn and Constraints on Enceladus' Thermal State from Astrometry[END_REF] found a best fit for the moons' positions for the tidal parameters

k 2 Q = (2.3 ± 0.7) × 10 -4 (6)
which yields Q ≈ 1700 for a given k 2 = 0.390 [START_REF] Gavrilov | Love numbers of the giant planets[END_REF]. Those values were found at the tidal frequency of Enceladus, Tethys and Dione. Later, using Cassini data [START_REF] Cooper | The Caviar software package for the astrometric reduction of Cassini ISS images: description and examples[END_REF], more calculations done independently by JPL confirm those low values of Q for the inner satellites but moreover, [START_REF] Lainey | New constraints on Saturn's interior from Cassini astrometric data[END_REF] announced a value of around 300 for Rhea, revealing the frequency-dependent nature of the tides in Saturn. Such revelations about the tidal migration of the icy-moons brought more insight about the dissipative processes acting in the planet [START_REF] Remus | The equilibrium tide in stars and giant planets. I. The coplanar case[END_REF], and furthermore opened the door to new considerations about the past history of saturnian moons. Saturnian satellites are today thought to be much younger than the Solar System and from this idea, [START_REF] Charnoz | Accretion of Saturns mid-sized moons during the viscous spreading of young massive rings: Solving the paradox of silicate-poor rings versus silicate-rich moons[END_REF] proposed that the inner midsized moons were created from a primordial massive ring [START_REF] Salmon | Long-term and large-scale viscous evolution of dense planetary rings[END_REF] one by one, assuming that the quality factor for Saturn was below 2000.

FIGURE 5: Satellite mass as a function of their position. This distribution of matter around the planet would be a consequence of satellites forming one after another from a massive ring and migrating fast. Rhea would be created first from the initial ring and migrate outward, then Dione formed from a less important reservoir which explains its smaller mass compared to Rhea. After migrating and eventually colliding with previously formed bodies, satellites would arrange their orbits as we see them today, leaving the initial massive ring with a final mass of the order of Mimas' mass.

Also strong tidal dissipation inside Saturn was brought forward in order to explain the eccentricity equilibrium of Enceladus Lainey et al., 2012. FIGURE 6: Saturn and Enceladus quality factors for preserving the eccentricity equilibrium state of Enceladus. Tidal heating of Enceladus was revealed by Cassini in 2006 [START_REF] Porco | Cassini Observes the Active South Pole of Enceladus[END_REF] and the amount of energy dissipated inside the satellite could not explain its eccentric orbit (Table 1). A high momentum exchange from the planet due to tidal energy dissipation inside the planet explains this equilibrium.

The figure was taken from [START_REF] Lainey | Strong Tidal Dissipation in Saturn and Constraints on Enceladus' Thermal State from Astrometry[END_REF].

Titan and Iapetus

In the linage of those last results, this thesis is dedicated to reformulate a past history of the System of Saturn. Even though, results stand mainly for the icy satellites, we have assumed all along this thesis a low value of Q for Titan with the aim of studying a past 5:1 mean motion resonance between Titan and Iapetus. Today their period ratio is close to 4.97, meaning that the resonance could have been crossed not long ago. This resonance crossing is possible if Q is under 14 000 at the frequency of Titan.

A first argument for the importance of this study is the fact that Titan is the most massive satellite in this system and it should have played a major role in the early evolution of the system. A history of Saturn can not be done without this important body. Secondly, Iapetus' orbit has several unexplained features, such as its eccentricity (around 0.03) and a 8 degrees orbital tilt with respect to its local Laplace plane (Nesvorný et al., 2014, chapter 4 is dedicated to the Laplace plane). Some authors have proposed some plausible scenarios to explain inclination. [START_REF] Ward | Orbital inclination of Iapetus and the rotation of the Laplacian plane[END_REF] invoked a fast dispersal of the circumplanetary disk from which satellites were formed. Iapetus would have evolved to the orbit we see today if the disk was dispersed in around 10 2 or 10 3 years, which is a fast time scale for this process to happen. Also, in the context of the early solar system instability [START_REF] Tsiganis | Origin of the orbital architecture of the giant planets of the Solar System[END_REF], [START_REF] Nesvorný | Excitation of the Orbital Inclination of Iapetus during Planetary Encounters[END_REF] argue that Iapetus could have been excited by close encounters between an ice giant and the system of Saturn. We are looking here for an alternative scenario where strong tidal dissipation and the 5:1 resonance crossing would be responsible for the orbital behaviour of Iapetus.

The first part of this thesis is dedicated to theoretical modelling. The first chapter is dedicated to tides, inside energy dissipation of the planet and the consequential tidal acceleration of moons. The second chapter deals with basic orbital dynamics, such as the N-Body problem and the planetary gravitational field. Analytical development of the equations of motion is outlined afterwards and constitutes a third chapter. The last chapter of this first part presents the equations of the local Laplace plane. Results of the 5:1 mean motion resonance crossing are presents in the second part of the thesis.

Part I Theoretical framework

Chapter 1

Tides and tidal acceleration in gravitational systems

This first chapter is dedicated to tidal interactions. We will elaborate the main ideas and equations behind tides and draw a general overview. The first approach concerning dynamical astronomy is to consider a given number of point-like masses which undergo mutual gravitational interaction between one another. In an inertial frame, a mass m would create a gravitational potential per unit mass

U = - Gm r (1.1)
where G stands for the gravitational constant and r the radial distance from the center of the body. From there, bodies surrounding the first one will undergo an acceleration r = -Gm r 3 r (1.2) and all bodies will evolve in time, attracting each other according to the inverse square law1 . This is the basis of orbital dynamics involving several masses and the study of such systems lies in what is called the N-Body problem.

Extended levels of complexity can be introduced from there, for instance by adding non-gravitational forces, such as radiation pressure or atmospheric drag. Those forces are known to be important in studying the motion of an artificial satellite in low Earth orbit. Furthermore one can consider taking into account the more complex gravitational field of bodies. The inverse square law applies to perfectly spherical bodies or point-like masses, but departure from it will induce a different force. If the shape of the body undergoing a force is considered then one must take into account that the gravitational potential will not be identical throughout its volume. This will lead to tidal interactions between bodies and the aim of this chapter is to make a general overview of modern tidal theories and then to outline the basic equations governing the dynamical states of bodies, in terms of their orbital evolutions and rotational motions

General overview 1.Differential acceleration, the source of tides

We first consider two bodies with different masses in the inertial frame of their centre of mass. Let there be for instance a satellite of mass m orbiting a more massive planet with mass M . The gravitational pull of the satellite acting on the closest side of the planet is greater than the pull acting on its centre, which is greater than the pull acting on its farthest side. This is the simple consequence of the inverse square law. m M FIGURE 1.1: Simple schema of the field applied on the different regions of the primary. Arrows correspond to instantaneous forces applied on the body.

The resulting consequences can be separated into two effects. The first and more classic effect will see the primary undergoing an acceleration (in an inertial frame), just like any point-like masses. The second and more interesting will be a distortion of its shape. Indeed, seen in the frame of the primary, the force field at the surface of the body will look like FIGURE 1.2: Seen in the frame of the primary, the surface will stretch under the force field produced by the satellite's gravity.

Central parts will also undergo similar forces.

Tidal bulge

As bodies are not perfectly rigid, this force field will create a tidal bulge which is a priori aligned with the satellite's direction. But many features have been neglected so far. First, the primary is usually spinning around its rotational axis and the satellite travels on its orbit around the primary. Therefore, in the frame of the rotating planet, the tidal bulge is moving around the primary with a frequency

ω tides = Ω p -n (1.3)
where Ω p = 2π P is the rotational spin frequency of the planet, with P its period and n is the mean motion of the satellite. Secondly, and most importantly, the tidal bulge does not instantaneously align itself with the satellite's direction but undergoes a lag, due to the non-elasticity of the primary body. As the tidal force field is constantly changing in time for an element of mass in the planet, the body never fully reaches a state of equilibrium but constantly changes its shape, chasing the tidal equilibrium. For a system, where the satellite orbits beyond the synchronous orbit2 , the tidal bulge is in advance with respect to the satellite's direction. On the contrary, if a satellite is below it, the bulge tends to have a delay. This is the case for the martian satellite Phobos. The tidal bulge, driven by the faster rotation of the primary is ahead of the direction of the satellite. The angle β represents the geometrical lag between the satellite's direction and the axis of the bulge.

Energy dissipation and angular momentum exchange

Depending on the structure of the primary, the strain on the primary and its reaction to it will produce mechanical energy dissipation which is restored in the form of heat. This energy loss is due to the non-elasticity of the solid matter in the case of telluric bodies, but one can also find viscous dissipation in the case of fluid layers. The result of the misalignment of the bulge with the satellite's direction will produce an orbital torque on the satellite. Indeed, the part of the bulge closer to the satellite will create an attraction stronger than the farthest part of the bulge and this imbalance will then create a resulting force with one of the components parallel to the satellite's motion. This torque adds momentum to the satellite's motion and increases its semi-major axis and eccentricity.

Formalism m M

Or bi ta l m ot ion F t FIGURE 1.4: Scheme showing how the tidal bulge adds a perturbing force to the motion of the satellite. The resulting force exerted by the two ends of the bulge has a component transverse to the satellite's motion. This is due to the shorter distance separating the satellite and the closest part of the bulge.

To balance the gain of momentum of the satellite, the planet's spin undergoes an opposite reaction by losing rotational momentum. In other words, the planet spins down due to tides raised by its satellite

Formalism

Using the notation in Kaula, 1964 3 , the potential generated by a moon (Equation 1.1) can be replaced in the reference frame of the planet by its disturbing function

U 0 = Gm * |r -r * | (1.4)
which can be expanded as a sum of zonal harmonics

U 0 = Gm * r * ∞ l=2 r r * l P l (cos(S)) (1.5)
where the superscripts * denotes the variables and parameters of the moon. r is here the radial distance from the planet's centre where the potential is calculated, S is the angle between the moon's direction and r and P l is the l th order Legendre polynomial. If the size of the planet is negligible compared to the distance of the moon, then one can truncate the series, keeping only the term associated with l = 2, as done in [START_REF] Mignard | The evolution of the lunar orbit revisited. I[END_REF]. In Kaula's paper [START_REF] Kaula | Tidal Dissipation by Solid Friction and the Resulting Orbital Evolution[END_REF], this potential is expressed in terms of the moon's orbital elements, and in his second article (Mignard, 1980), Mignard keeps a more compact form which makes use of the position of the moon. The planet deforms itself as a response to this potential and this adds an extra tidal potential to the planet's gravitational field, which makes use of the Love numbers k l

U (r) = ∞ l=2 k l Gm * R 2l+1 p r l+1 r * l+1 P l (cos(S)) , (1.6)
where cos(S) = r.r * rr * This potential would be synchronous with the moon's position. But, as we outlined in the previous section, the deformation response of the planet is delayed with respect to the moon's potential. The delayed response can be illustrated by introducing a fictitious satellite. This body would raise instantaneous tides on the planet and would be lagging behind the real tide raising satellite with a position

r f = r * + f (1.7)
where the position lag

f is f = ∆t (Ω p × r * -v) . (1.8) m* M Or bi ta l m ot ion f r f r * FIGURE 1.6:
The lag f is usually not constant in the case of an eccentric orbit for instance but some approximation can be made

f = ∆tr|Ω p -n|
Love Numbers are usually difficult to obtain either from theory or observation. Only the k 2 is usually considered in tidal theories, so we neglect here higher orders. Also we substitute r with r * because we are looking for the potential at the real moon's position and r * is substituted with r f . Then following [START_REF] Efroimsky | Physics of bodily tides in terrestrial planets and the appropriate scales of dynamical evolution[END_REF] the force, F = ∇U , has the expression

F = - 3k 2 Gm 2 R 5 p r 10 rr 2 -f r 2 -2r (r.f ) (1.9)
from which, after using Equation 1.8, the acceleration becomes

a = - 3k 2 G (m S + m * ) m * R 5 p r 10 m S ∆t 2r (r.v) + r 2 (r × Ω + v) (1.10)
where ∆t is the time lag defined by [START_REF] Lainey | New constraints on Saturn's interior from Cassini astrometric data[END_REF])

∆t = arctan(1/Q) 2π T (1.11)
and T is the period of the tidal excitation.

T = 2π 2|Ω p -n| .
(1.12)

Energy dissipation in the primary

Several parameters were introduced in the tidal acceleration. First, the Love Number k 2 appears as a constant in the equation. But the assumption that the tidal potential raised by a satellite is independent from the satellite's orbit or its mass can not always hold. Instead one can assume a frequency dependant tidal response, for which k 2 is a function of the main tidal frequency χ = 2|Ω p -n|.

(1.13) Also, one important feature of the tidal theory is the quality factor Q introduced in the time lag (Equation 1.11). It is a measure of energy dissipated during a tidal cycle [START_REF] Goldreich | Q in the Solar System[END_REF])

1 Q = 1 2πE 0 - dE dt dt (1.14)
where E 0 is the maximum energy stored in the tidal distortion. In the case of Saturn, tidal dissipation can have several sources. • Non-elasticity of the core. As explained in the previous section, the solid body will have a delayed response to the tidal potential of the moon and the reason for this delay is the inelastic nature of the body. Energy is lost as the body reshapes itself continuously.

• Viscous friction in gaseous envelope. For giant planets, the existence of fluid envelopes adds a component to tidal response. The fluid response to the tidal potential of the satellite will follow the Navier-Stokes equations for fluid dynamics [START_REF] Remus | The equilibrium tide in stars and giant planets. I. The coplanar case[END_REF]. It will create an additional flow to the existing convective motion of the fluid envelope and their interaction produces viscosity. This will generate a toroidal velocity field and change the gravity field of the planet.

It is not yet known wether tidal dissipation primarily comes from the solid core or the convective envelope.

Constant lag versus frequency dependent tidal response

Equation 1.10 stands for the tidal acceleration undergone by the satellite. As shown in [START_REF] Efroimsky | Physics of bodily tides in terrestrial planets and the appropriate scales of dynamical evolution[END_REF], the secular effect on the satellite's semi-major axis, assuming a low inclination and a low eccentricity, is

da dt = - 6k 2 R 5 p nm∆t m S a 4 (n -Ω p ) (1.15)
which can be rewritten

da dt = G(m S + m) 3k 2 R 5 p m m S a 11/2 arctan 1 Q (1.16)
after substitution of the time lag (Equation 1.11) into the previous equation4 . Now, several assumptions can be made. Like previously stated, the tidal response of the planet can consist in a constant time lag with respect to the satellite's direction. Or, in a more complex fashion, the response can be dependant on several parameters such as the orbital period. Using this approach, tidal parameters k 2 and Q will then be frequency dependent and [START_REF] Lainey | New constraints on Saturn's interior from Cassini astrometric data[END_REF] shows evidence for such a behaviour. However, if we assume constant parameters in Equation 1.16, then the differential equation for the semi-major axis is separable in time and semi-major axis and therefore solvable.

a(t) = 39 2 G(m S + m) k 2 R 5 p m m S arctan 1 Q t + a 13/2 0 2/13
(1.17) This approach can give us a first insight into the possible values that Q can take. Table 1.1 indicates the value of Q needed for each satellite to migrate from the synchronous orbit to its semi-major axis in the age of the solar system, using Equation 1.17. [START_REF] Goldreich | Q in the Solar System[END_REF] derived Q ≈ 60 000 for the same consideration for Mimas, but used a different value for k 2 . Also in the case of Titan, values between 2000 and 4000 were brought forward in order to explain the advanced state of the resonance between Titan and Hyperion [START_REF] Colombo | On the formation of the orbit-orbit resonance of Titan and Hyperion[END_REF]. Invoking atmospheric tidal response in Saturn, Greenberg, 1973 also argues for values around 4000. It seems that it was difficult to have an agreement on Q for all satellites. Constant tidal lag is a convenient model to make calculation but the reality may be more complex. Studies of frequency dependant tidal response have been brought forward in the last decade [START_REF] Efroimsky | Physics of bodily tides in terrestrial planets and the appropriate scales of dynamical evolution[END_REF] and specially for stars and fluid envelops of giant planets [START_REF] Ogilvie | Tidal Dissipation in Rotating Giant Planets[END_REF][START_REF] Auclair-Desrotour | Impact of the frequency dependence of tidal Q on the evolution of planetary systems[END_REF][START_REF] Remus | The equilibrium tide in stars and giant planets. I. The coplanar case[END_REF]. Very recently, Fuller, Luan, and Quataert, 2016 made an explicit model of tides acting in the fluid envelope of the planet. Their model predicts frequency dependant quality factors for Saturn's moons, which are comparable to those found in [START_REF] Lainey | Strong Tidal Dissipation in Saturn and Constraints on Enceladus' Thermal State from Astrometry[END_REF]Lainey et al., 2017 (Figure 1.10). The evolution of a moon's orbit is dependant on a frequency dependant quality factor

Q ≡ 3k 2 m M R p a 5 nt tides (1.18)
where k 2 is the Love number, R p is the radius of the planet and t tides the time scale of the orbit expansion. The mechanism involved is resonance locking which can be reached due to the structure evolution of the planet and the consequential evolution of the oscillation mode frequency in time. During a lock, an oscillation mode stays resonant with the tidal forcing produced by a moon and increases greatly the tidal dissipation, which in turn, produces fast outward migration. In the frequency domain, the satellite would "surf" on dissipation peaks. The orbital evolution consists of a rapid evolution of a satellite's semi-major axis in such a way that nearby orbits diverge. This could go against the formation of a resonance for a pair of satellites, as it is understood that convergent orbits are needed to get a resonance capture in mean motion.

Summary and Conclusion

We have explained the tidal theory for a moon raising tides on a planet. After explaining the basic concepts, we have discussed the equations governing the motion of a moon undergoing a tidal acceleration. Most of the mathematics exposed in this chapter comes from the literature. Although different theories predict different evolutions, we have chosen to implement Equation 1.10

a = - 3k 2 G (m S + m * ) m * R 5 p r 10 m S ∆t 2r (r.v) + r 2 (r × Ω + v)
in our software for the tidal interaction, because we only need a constant increase of the semi-major axis of Titan to simulate the mean motion resonance crossing. Fuller's model predicts a quality factor at Titan's frequency of around 20 assuming a tidal time scale of around 2 billion years (Equation 1.18). For Iapetus, Fuller's model also predicts a fast migration but we find the rate unrealistic (Q ≈ 2 × 10 -4 ). However, we assume for this thesis that Iapetus does not migrate.

Lately many theories have argued in favour of a frequency dependant nature of the tidal response of a planet [START_REF] Fuller | Resonance locking as the source of rapid tidal migration in the Jupiter and Saturn moon systems[END_REF] and astrometric results [START_REF] Lainey | New constraints on Saturn's interior from Cassini astrometric data[END_REF] confirm those findings. It is obvious that a general formulation for tides is still to come, but the recent studies brought new insights to the complexity of the problem.

Chapter 2

Equations of motion for the system of Saturn

This chapter is dedicated to some basic equations for dynamical astronomy. The reader will find the general equations governing the motion of a satellite around Saturn. We will derive the equations for the relative acceleration of a satellite with respect the planet as well as approximations for the icy satellites. The gravity field of the planet will also be considered.

Saturnian system as a N-Body system

As for the previous chapter, we start with the well known Newtonian gravitational interaction between bodies. Just for the very beginning, we assume being in a inertial frame, in which we consider Saturn and two other bodies which can be satellites or the Sun. If we denote by r S , r i and r j the positional vectors of Saturn and bodies i and j, in the Galilean frame of the barycentre of mass verify M r S + m i r i + m j r j = 0

(2.1) and be differentiating with respect to time, we see that the same relation will hold for velocities and accelerations. Here, M denotes Saturn's mass and m i and m j are both inertial masses of the bodies. Then the Newtonian acceleration for Saturn is

rS = - Gm i |r S -r i | 3 (r S -r i ) - Gm j |r S -r j | 3 (r S -r j ) (2.2)
and for bodies i and j, the acceleration will be similar, we just need to swap indexes

ri = - GM |r i -r S | 3 (r i -r S ) - Gm j |r i -r j | 3 (r i -r j ) (2.3) rj = - GM |r j -r S | 3 (r j -r S ) - Gm i |r j -r i | 3 (r j -r i ) .
(2.4)

Those differential equations are representative of what one calls a three-body problem. The reference frame considered enables us to express those laws of motion in their simplest form. However, in our case, the Sun plays a major role in satellite dynamics and if considered in a certain model, the barycentre of mass will then be very close to the center of the Sun, due to its superior mass. This choice is not convenient if we study the orbits of natural satellites around a planet. For the rest of this work, we choose to express the dynamics of satellites in the frame of Saturn. In such a frame, the position, velocity and acceleration of Saturn are, by definition, null. This is a simple consequence of setting the center of the reference frame to Saturn itself. The position of the body i, r ′ i , with respect to Saturn is then

r ′ i = r i -r S .
(2.5) Differentiating twice with respect to time and using expressions for accelerations (2.2) and (2.4), the acceleration of the body i in the saturnian frame is

r′ i = - G(M + m i ) r ′ 3 i r ′ i + Gm j r j -r i |r j -r i | 3 - r j -r S |r j -r S | 3 .
(2.6)

Here, we use the translation law for j, as for the body i

r ′ j = r j -r S (2.7)
and also the fact that

r j -r i = r ′ j -r ′ i .
(2.8)

Finally, we denote r ij = |r jr i | and we redefine the unprimed vectors as being the one set in the frame of Saturn. We obtain the relative form of the equations of motion for the body i around Saturn perturbed by the body

j ri + G(M + m i ) r 3 i r i = Gm j r j -r i r 3 ij - r j r 3 j .
(2.9)

The first term in the right-hand side is known as being the direct acceleration due to the gravitational attraction of the perturber j. The second term is called the indirect part and represents the acceleration of Saturn due to the attraction of the body j. It's a manifestation of ficticious forces acting on a body due the choice of a non-inertial frame. In the presence of a fourth body and more, the right hand side will simply be a summation taken over all the masses

ri + G(M + m i ) r 3 i r i = N -1 j=1,j =i Gm j r j -r i r 3 ij - r j r 3 j .
(2.10)

In the absence of a third mass perturbing i, its acceleration will be the equation (2.3) without the second term involving the mass of j. Here, the only difference between the motion described in the inertial and the planeto-centric frame is the gravitational parameter, µ, which is GM in the inertial frame, and G(M +m i ) for the frame centred on the planet. This is known as the two-body problem1 and the solution of this equation is the well-known conic

r = a 1 -e 2 1 + e cos(f ) (2.11)
In the general equation (2.9), the perturbing part on the right hand side is expressible in terms of a potential. We multiply equation (2.9) by the infinitesimal vector

dr i ri + G(M + m i ) r 3 i r i .dr i = Gm j r j -r i r 3 ij - r j r 3 j .dr i .
(2.12)

For the left hand part of the previous equation, we change dr i = ṙi dt, and we have that ri .dr

i = ri .ṙ i dt = d 1 2 ṙ2 (2.13)
and

G(M + m i ) r i dr i r 3 i = d -G(M + m i ) 1 r i (2.14)
and for the right hand side, recalling that r ij = (r jr i ) . (r jr i )

r j -r i [(r j -r i ) . (r j -r i )] 3/2 - r j r 3 j = ∇ r i 1 (r j -r i ) . (r j -r i ) - r i .r j r 3 j (2.15)
At this point, we introduce the disturbing function for the body j acting on i

R j = Gm j 1 r ij - r i .r j r 3 j (2.16)
so that the acceleration of the body i orbiting Saturn and perturbed by the body j can be written

ri + G(M + m i ) r 3 i r i = ∇ r i R j .
(2.17)

The same analysis applies to the body j, perturbed this time by the body i, so that we have similar equation for j

rj + G(M + m j ) r 3 j r j = ∇ r j R i (2.18)

Saturn's gravity field

The gravitational field of Saturn and its moons were greatly improved with the Cassini mission. The Radio Science Subsystem on-board the spacecraft permitted precise position and velocity measurements during the whole time of the mission.

The gravitational potential of an object can be be expressed in terms of spherical harmonics

V = - Gm S r 1 - ∞ i=2 J i R p r i P i (sin(α)) (2.19)
where α = z r with z the third component of the position of the satellite and r the distance from the body's centre. Chapter 2. Equations of motion for the system of Saturn Gravitational spherical harmonics were determined up to J 10 for the even coefficients and up to J 3 for the odd coefficient. However, in this thesis we only kept the major perturbing part of the gravitational field which is the force associated with J 2 . Due to the great distance of Titan and Iapetus from Saturn, all other harmonics are neglected.2 . The disturbing function associated with the J 2 effect is

Parameter

R J 2 = - G(m S + m i )J 2 R 2 p 2 1 r 5 3z 2 -r 2 (2.20)
and the acceleration r = ∇ r R J 2 , written in the form of its components is

ẍ = 3 2 G(m S + m i )J 2 R 2 p 1 r 7 5z 2 -r 2 x, ÿ = 3 2 G(m S + m i )J 2 R 2 p 1 r 7 5z 2 -r 2 y, z = 3 2 G(m S + m i )J 2 R 2 p 1 r 7 5z 2 -3r 2 z.
(2.21)

Inner Moons

The gravitational interaction between the inner moons and the outer moons was simplified. Instead of including the gravitational acceleration of inner moons in the system, we have used an approximation which consists in adding the gravitational contribution of the inner moons to the second order harmonic of the planet. This leads to defining an upgraded J 2 coefficients [START_REF] Tremaine | Satellite Dynamics on the Laplace Surface[END_REF])

J ′ 2 = J 2 + 1 2 i a i R p 2 m i m S (2.22)
By using this approximation, we get rid of low frequencies which would appear in the system. Indeed, If Mimas orbit was added, the lowest time scale of the system would be around 1 day, which corresponds to Mimas' orbital period. Instead, the lowest time scale is the one associated with Titan's orbital period which is around 16 days. Timesteps in numerical integrations are then greatly increased and therefore numerical integrations get faster. Also, the secular effects of all inner moons are preserved Satellite

Summary of forces on a moon

To summarize our numerical model for the study of the 5:1 mean motion resonance between Titan and Iapetus, we are integrating the positions and velocities of Titan, Iapetus and the Sun. The gravitational pull of all satellites from Mimas to Rhea were approximated by adding their contribution to the J 2 coefficient of Saturn (Equation 2.22). However Hyperion is sometimes integrated for a few simulations in order to study the effect of the 5:1 mean motion resonance crossing on its orbit (we recall that it is actually in a resonance with Titan). Jupiter is also added sometimes, but most of our work is done using 4 bodies. It is computed using cos(J) = cos(i) cos(i LP ) + sin(i) sin(i LP ) cos(Ω -Ω ⊙ ) where the LP stands for the elements of the Laplace plane (chapter 4), Ω ⊙ is the ascending node of the Sun and J is the tilt. Therefore a direct N-Body integration of the four bodies considered shows a good agreement with an integration made with all the bodies. Figure 2.1 shows the behaviour of the orbital plane of Iapetus. We can see that its inclination oscillates around a forced value (see chapter 4 for the Laplace plane) and has a 8 degree tilt to it (Figure 2.2).

Summary of the numerical model

Using a Gauss-Radau integrator (Appendix A), direct numerical integrations of the 4-Body system are performed (Saturn, Titan, Iapetus and the Sun). N-Body interactions was implemented using Equation 2.10 and the J 2 effect of Saturn is added via Equation 2.21 using Equation 2.23 for the value of J 2 . Added to these basic forces is the tidal acceleration (Equation 1.10) defined in the previous chapter.

Chapter 3

Semi-analytic modelling

Integrating numerically the equations of motion is the most direct method to solve the motion of the moons. However, this method is usually machine-intensive and time consuming. To avoid making long numerical simulations, we have decided to use a semi-analytic approach, along with the N-Body code, by first expanding the disturbing function in terms of the elliptical elements of the bodies involved, then averaging the expansion over the elements varying rapidly, and finally integrating numerically the averaged equations. By doing so, we are getting rid of the short period variations of the orbits, which makes the dynamical time scale much longer and therefore the integration time much smaller. Several kinds of perturbations are taken into account. The first one, is the averaging of the flattening effect of Saturn and the inner satellites. Then we will handle the averaged effect of the Sun acting on the satellites. Finally, we will deal with the mutual gravitational interaction between Titan and Iapetus.

General statements

The disturbing function as a sum of cosines

Starting from the main perturbing acceleration (Equation 2.17), our main interest in defining the disturbing function is to express it in terms of the orbital elements of the body and the perturber. This task was subjected to centuries of research and enhancement after the basis of orbital mechanics was laid by Isaac Newton in his Principia Mathematica [START_REF] Newton | Philosophiae naturalis principia mathematica[END_REF]. A few notable works in which the authors have been working on expanding the disturbing function are [START_REF] Laplace | Théorie de Jupiter et de Saturne[END_REF][START_REF] Delaunay | Théorie du mouvement de la lune[END_REF]. In general, one has Equation 2.17

ri + G(M + m i ) r 3 i r i = ∇ r i R j
for a perturbing body j acting on a body i. Now, as we have shown in chapter 2 or in Appendix C, the disturbing function is slightly different whether one studies a body perturbing an inner or an outer body. However, one can show that, using several expansion technics, such as Fourier or Taylor series, a general disturbing function R can be expanded as a trigonometric series

R = S(a, a ′ , e, e ′ , i, i ′ ) cos j 1 λ ′ + j 2 λ + j 3 ̟ ′ + j 4 ̟ + j 5 Ω ′ + j 6 Ω (3.1)
where functions S, playing the role of amplitudes, are functions of the metric elements of both bodies, namely the semi-major axis, the eccentricity and the inclination, whereas the arguments of the cosine terms are linear combinations of the angular elements : the mean longitude, the longitude of the pericenter and the longitude of the ascending node. Different types of expansions can be made but there are always three different constrains on the amplitudes and arguments.

• An argument is associated with a set of six j coefficients which have to verify what one calls the d'Alembert relation [START_REF] Murray | Solar System Dynamics[END_REF] 6

i=1 j i = 0 (3.2)
• Furthermore, there is a constrain on the coefficients associated with longitudes of ascending nodes, namely

|j 5 | + |j 6 | is even or zero (3.3)
• And finally we have an important property relating an amplitude to its argument named the d'Alembert characteristics [START_REF] Brouwer | Methods of celestial mechanics[END_REF].

The lowest power of an eccentricity/inclination function 1 in the amplitude equals to the absolute value of the j coefficient associated with its pericenter/ascending node.

After an expansion, one can distinguish terms acting on a short time scale, or with high frequencies and others whose effects govern the long term behaviour of the system; one names those secular terms. The idea is to average out the high frequencies of the system. Satellites mean anomalies usually have the shortest period, therefore they are always averaged out except for the combination 5λ 8λ 62 which has a long period and is the central part of our dynamics. The motion of the Sun around Saturn has a period of almost 30 years (Table 3.1). It acts like an intermediate frequency between mean anomalies and the long varying angles such as the longitudes of pericentre and the motion of orbital planes around the Laplace equilibrium (see chapter 4). Initially the sun's mean longitude was averaged out but we have kept a few terms in subsequent simulations Table 3.3.

Averaging the flattening effect

We have seen in the previous chapter that the potential of the planet Saturn can be expanded as (Equation 2.19)

V = - G(M + m) r + G(M + m) r 3 J 2 R 2 p P 2 (sin(α)) + G(M + m) r 5 J 4 R 4 p P 4 (sin(α)) + ... (3.4) 1 namely s = sin i
where the first term stands for the point-like interaction of Saturn with its satellites. The next term is the most important one describing the flattening of the planet. Any other terms are neglected for the motion of Titan and Iapetus due to their distance to the planet3 . Recalling that r sin(α) = z, we write the disturbing function associated with the flattening

R J 2 = - G(M + m) r 3 J 2 R 2 p 1 2 3 z r 2 -1 ; (3.5)
then replacing z by its expression in terms of the orbital elements (Equation C.8) and r by the expression of a conic Equation 2.11, we rewrite Equation 3.5

R J 2 = - 1 2 J 2 R 2 p G(M + m) a 3 (1 -e 2 ) 3 (1 + e cos(f )) 3 3 (sin(i) sin(ω + f )) 2 -1 . (3.6)
This function is 2π-periodic in the variable f , the true anomaly. However, we want the average value over one period in time, not geometrically. Therefore the average is made over the mean anomaly

R J 2 = 1 2π 2π 0 - 1 2 J 2 R 2 p G(M + m) a 3 (1 -e 2 ) 3 (1 + e cos(f )) 3 3 (sin(i) sin(ω + f )) 2 -1 dM.
(3.7) Using useful relations between anomalies [START_REF] Murray | Solar System Dynamics[END_REF] we have

dM = 1 √ 1 -e 2 r a 2 df. (3.8) Therefore R J 2 = - 1 2 J 2 R 2 p G(M + m) a 3 (1 -e 2 ) 3/2 1 2π 2π 0
(1 + e cos(f )) 3 (sin(i) sin(ω + f )) 2 -1 df (3.9) which integrates to

R J 2 = G(M + m)J 2 R 2 p 8a 3 (1 -e 2 ) 3/2 (1 + 3 cos(2i)) (3.10)

Solar terms

Secular

The acceleration of a satellite due to the gravity of the Sun is expressed as

rs + G(M + m s ) r 3 s r s = ∇ rs R ⊙ (3.11)
where the disturbing function of the Sun is

R ⊙ = Gm ⊙ 1 |r ⊙ -r s | - r s .r ⊙ r 3 ⊙ . (3.12)
This function can be expanded as a sum of Legendre polynomial as shown in Appendix C and in Murray and Dermott, 2000

R ⊙ = Gm ⊙ r ⊙ ∞ l=2 r s r ⊙ l P l (cos(Ψ)). (3.13)
Now, since the Sun's distance from Saturn is big with respect to any saturnian satellite's position around the planet, the ratio rs r ⊙ is small during the evolution, one can neglect the terms in the expansion for l greater than 3. Therefore the expression for the Sun's disturbing function becomes

R ⊙ = n 2 ⊙ a 2 s r s a s 2 a ⊙ r ⊙ 3 1 2 3 cos 2 (Ψ) -1 (3.14)
where we have made use of Kepler's third law to introduce the mean motion of the Sun

G (m ⊙ + m S ) ≈ Gm ⊙ = n 2 ⊙ a 3 ⊙ . (3.15)
rs as 2 and a ⊙ r ⊙ 3 are replaced with power series of eccentricities (Equation C.13 and Equation C.14) and cos(Ψ) is replaced with a sum of cosines of angular elements (Equation C.9). After these last replacement and collecting all cosine terms we average the series over the mean longitude of the satellite and the Sun.

R ⊙ = 1 2π 2π 0 1 2π 2π 0 R ⊙ dλ s dλ ⊙ (3.16)
This yields eight different cosine terms exposed in Table 3.2. We write

R ⊙ = n 2 ⊙ a 2 s 8 i=1 S i cos (φ i ) (3.17) i Argument (φ i ) Amplitude (S i ) ∅ 1 32 8 + 12e 2 ⊙ + 15e 4 ⊙ + 6e 2 s 2 + 3e 2 ⊙ 1 -6s 2 s + 6s 4 s 1 -6s 2 ⊙ + 6s 4 ⊙ 2Ω s -2Ω ⊙ 3 8 8 + 12e 2 ⊙ + 15e 4 ⊙ + 6e 2 s 2 + 3e 2 ⊙ s 2 s 1 -s 2 s s 2 ⊙ 1 -s 2 ⊙ Ω s -Ω ⊙ 3 8 8 + 12e 2 ⊙ + 15e 4 ⊙ + 6e 2 s 2 + 3e 2 ⊙ s s 1 -s 2 s 1 -2s 2 s s ⊙ 1 -s 2 ⊙ 1 -2s 2 ⊙ 2ω s -2Ω ⊙ 15 8 e 2 s 2 + 3e 2 ⊙ 1 -s 2 s 2 s 2 ⊙ 1 -s 2 ⊙ 4Ω s -2Ω ⊙ -2ω s 15 8 e 2 s 2 + 3e 2 ⊙ s 4 s s 2 ⊙ 1 -s 2 ⊙ 3Ω s -Ω ⊙ -2ω s 15 4 e 2 s 2 + 3e 2 ⊙ s 3 s 1 -s 2 s s ⊙ 1 -s 2 ⊙ 1 -2s 2 ⊙ Ω s + Ω ⊙ -2ω s -15 4 e 2 s 2 + 3e 2 ⊙ s s 1 -s 2 s 3/2 s ⊙ 1 -s 2 ⊙ 1 -2s 2 ⊙ 2ω s -2Ω s 15 8 e 2 s 2 + 3e 2 ⊙ s 2 s 1 -s 2 s 1 -6s 2 ⊙ + 6s 4 ⊙ TABLE 3
.2: Secular terms for the perturbation of the satellite by the Sun. Note that no approximation was made concerning both inclinations. The subscript s stands for either Titan or Iapetus.

Dependence on the Sun's mean longitude

Equation 3.16 averages the disturbing function over the longitude of the satellite and the Sun. But one can also average on the satellite's mean longitude only, avoiding further approximation for the Sun. The reason is that the period of the Sun is quite big compared to the satellite's period and acts on a time scale between the fast orbital period of the satellites and the slow variations of their nodes and pericentres. After averaging only over the satellite's mean longitude, we are left with the eight secular terms Table 3.2 plus some terms involving the mean longitude of the Sun in the form k × λ ⊙ . Because of the large number of terms and the fact that k × λ ⊙ has a short period when k is big, we only choose terms in 1 × λ ⊙ . Table 3.3 gathers all of theses terms. We have an integration with the semi-analytic code comprising terms in Table 3.2 (in gold) and a simulation integration using in addition terms in Table 3.3 (in green). These evolutions do not represent any real dynamics but are test-simulations to assess the reliability of the solar expansion.

R ⊙,λ ⊙ = n 2 ⊙ a 2 s 35 i=1 S i cos(φ i ). (3.18) i Argument (φ i ) Amplitude (S i ) 1 λ ⊙ + 2Ω s -3̟ ⊙ 1 32 e 3 ⊙ s 2 s 1 -6s 2 s 1 -s 2 ⊙ 2 2 λ ⊙ + Ω s + Ω ⊙ -3̟ ⊙ -1 16 e 3 ⊙ s s 1 -s 2 s 1 -2s 2 s s ⊙ 1 -s 2 ⊙ 3/2 λ ⊙ + 4Ω s -4Ω ⊙ -2̟ s + ̟ ⊙ -15 16 e 2 s e ⊙ s 4 s s 4 ⊙ λ ⊙ + 3Ω s -3Ω ⊙ -2̟ s + ̟ ⊙ -15 4 e 2 s e ⊙ s 3 s 1 -s 2 s s 3 ⊙ 1 -s 2 ⊙ λ ⊙ + 2Ω s -2Ω ⊙ -2̟ s + ̟ ⊙ -45 8 e 2 s e ⊙ s 2 s 1 -s 2 s s 2 ⊙ 1 -s 2 ⊙ λ ⊙ + Ω s -Ω ⊙ -2̟ s + ̟ ⊙ -15 4 e 2 s e ⊙ s s 1 -s 2 s 3/2 s ⊙ 1 -s 2 ⊙ 3/2 λ ⊙ -4Ω s + 2̟ s + ̟ ⊙ -15 16 e 2 s e ⊙ s 4 s 1 -s 2 ⊙ 2 λ ⊙ -4Ω ⊙ + 2̟ s + ̟ ⊙ -15 16 e 2 s e ⊙ s 4 s 1 -s 2 s 2 λ ⊙ -Ω s -3Ω ⊙ + 2̟ s + ̟ ⊙ 15 4 e 2 s e ⊙ s s 1 -s 2 s 3/2 s 3 ⊙ 1 -s 2 ⊙ λ ⊙ -2Ω s -2Ω ⊙ + 2̟ s + ̟ ⊙ -45 8 e 2 s e ⊙ s 2 s 1 -s 2 s s 2 ⊙ 1 -s 2 ⊙ λ ⊙ -3Ω s -Ω ⊙ + 2̟ s + ̟ ⊙ 15 4 e 2 s e ⊙ s 3 s 1 -s 2 s s ⊙ 1 -s 2 ⊙ 3/2
Figure 3.1 shows that the approximation of the solar gravitational pull is well represented by terms inTable 3.2 and Table 3.3.

Titan-Iapetus interaction

The largest effort was put towards the disturbing function concerning the gravitational interaction between Titan and Iapetus. The perturbation of Iapetus on Titan is written

r6 + G(M + m 6 ) r 3 6 r 6 = ∇ r 6 R 8 (3.19)
whereas similarly, the equation for Titan perturbing Iapetus is

r8 + G(M + m 8 ) r 3 8 r 8 = ∇ r 8 R 6 (3.20)
where the subscript 6 denotes the variables for Titan and 8 for Iapetus. The disturbing functions are

R 8 = Gm 8 1 |r 8 -r 6 | - r 8 .r 6 r 3 8 (3.21) R 6 = Gm 6 1 |r 6 -r 8 | - r 6 .r 8 r 3 6 . (3.22)
Both expressions are similar, except for the indirect parts which differ and the mass factors. During the evolution of the satellites, Iapetus will act like an external perturber for Titan, whereas Titan acts like an internal one

r 8 > r 6 (3.23)
Following [START_REF] Murray | Solar System Dynamics[END_REF], we rewrite both disturbing functions

R 8 = Gm 8 a 8 (R D + αR E ) (3.24) R 6 = Gm 6 a 8 R D + 1 α 2 R I (3.25)
where the direct part, R D is common to both functions

R D = a 8 |r 8 -r 6 | (3.26)
and α is the semi-major axis ratio α = a 6 a 8 .

(3.27)

Indirect parts differ. For Iapetus we have

R E = - r 6 a 6 a 8 r 8 2 cos(Ψ) (3.28)
while for Titan

R I = - r 8 a 8 a 6 r 6 2 cos(Ψ) (3.29)
where we recall that r 6 r 8 cos(Ψ) = r 6 .r 8 . The disturbing functions written in Equation 3.21 and Equation 3.22 can also be expanded in a series involving Legendre polynomials, just as the disturbing function for the Sun Equation 3.13. Both functions would read

R 8 = Gm 8 r 8 ∞ l=2 r 6 r 8 l P l (cos(Ψ)) (3.30)
and

R 6 = Gm 6 r 8 ∞ l=2 r 6 r 8 l P l (cos(Ψ)) + Gm 6 r 6 r 2 8 cos(Ψ) -Gm 6 r 8 r 2 6 cos(Ψ). (3.31)
But these expansions were only used to check our main expansion (see section 3.6) which we will derive here. Instead we have made an extensive use of explicit expressions of amplitudes associated with arguments. For the direct part of the disturbing function, one has (see Murray and Dermott, 2000 Chapter 6, p.247 and Equation C.23)

R D = imax i=0 2i! i! (-1) i 2 2i+1 α i × i s=s min nmax n=0 (2s -4n + 1)(s -n)! 2 2n n!(2s -2n + 1)! s-2n m=0 κ m (s -2n -m)! (s -2n + m)! × (-1) s-2n-m F s-2n,m,p (I 6 )F s-2n,m,p ′ (I 8 ) i-s l=0 (-1) s 2 2s (i -s -l)!l! × ℓmax ℓ=0 (-1) ℓ ℓ! ℓ k=0 ℓ k (-1) k α ℓ d ℓ dα ℓ b (j) i+ 1 2 (α) × X i+k,-j 2 -j 4 -j 2 (e 6 )X -(i+k+1),j 1 +j 3 j 1 (e 8 ) × cos [j 1 λ 8 + j 2 λ 6 + j 3 ̟ 8 + j 4 ̟ 6 + j 5 Ω 8 + j 6 Ω 6 ]
which gives an explicit expression for the amplitude associated to an argument φ = j 1 λ ′ + j 2 λ + j 3 ̟ ′ + j 4 ̟ + j 5 Ω ′ + j 6 Ω in terms of the semi-major axis ratio, both eccentricities and inclinations. Details for Equation C.23 can be found in Appendix C. For the indirect parts, terms of Equation 3.28 can be written

R E = -κ m (1 -m)! (1 + m)! F 1,m,p (I 6 )F 1,m,p ′ (I 8 )X 1,-j 2 -j 4 -j 2 (e 6 )X -2,j 1 +j 3 j 1 (e 8 ) × cos [j 1 λ 8 + j 2 λ 6 + j 3 ̟ 8 + j 4 ̟ 6 + j 5 Ω 8 + j 6 Ω 6 ]
and for Equation 3.29

R I = -κ m (1 -m)! (1 + m)! F 1,m,p (I 6 )F 1,m,p ′ (I 8 )X -2,-j 2 -j 4 -j 2 (e 6 )X 1,j 1 +j 3 j 1 (e 8 ) × cos [j 1 λ 8 + j 2 λ 6 + j 3 ̟ 8 + j 4 ̟ 6 + j 5 Ω 8 + j 6 Ω 6 ] .
Two different kinds of term have to be considered.

• Secular terms are here defined as terms being effective on large time scales. Their arguments do not contain any mean longitudes (j 1 = 0 and j 2 = 0).

• Resonant terms contain terms with 5λ 8λ 6 .

A first model we have made for the Titan-Iapetus gravitational interaction consisted in considering simply the fourth order expansion already made in [START_REF] Murray | Solar System Dynamics[END_REF] In Appendix B in this textbook, cosine arguments and amplitudes are already exposed according to their order in a generic way. Table B.1 at page 540 gathers all zeroth-order arguments and one can get our secular terms by plugging in j = 0. The same operation can be made for resonant terms, for which one has to plug in j = 5 in table B.16 at page 553. By doing so, we obtain 16 secular terms and 19 resonant terms for the direct part of the disturbing function. For indirect parts, only a single term comes into play. For the perturbation of Iapetus acting on Titan (Equation 3.28) one has only This model of the interaction between Titan and Iapetus seemed adequate for the study of the 5:1 resonance crossing, but several considerations on some approximation made during the derivation led us to consider more terms. Indeed, inclinations were assumed to be small so that only first, second, third and fourth powers of s = sin I 2 were kept. This can be done for Titan but the inclination of Iapetus is 15 degrees on average. This can not be considered small. This led to a redefinition of the order of the expansion. The first model was a fourth order expansion meaning that any term e a 1 8 e a 2 6 s a 3 8 s a 4 6 with

R E = - 3125
a 1 + a 2 + a 3 + a 4 > 4 (3.34)
was neglected. Here because s 8 is not negligible, we chose not to let a 3 play a role in the condition of the order. Subsequently, we have considered all possible terms where

a 1 + a 2 + a 4 ≤ 4 (3.35)
and because Titan stays close to the equator a 4 ≤ 2.

(3.36)

This last condition was used to limit the size of the expansion, as the number of terms and the size of their amplitudes got quite important.

Consequently, this new model included 31 secular terms and 61 resonant terms for the direct part.

Secular terms

In accordance with d'Alembert relation, secular terms will have the following constraint j 3 + j 4 + j 5 + j 6 = 0 (3.37) knowing that j 1 = j 2 = 0 already. Secular terms are composed of 31 cosine terms exposed Table 3.4 ( Φ 18,0,0 + e 2 6 Φ 18,0,1 + e 2 8 Φ 18,0,2 + s 2 6 Φ 18,0,3

R sec = 31 i=1 S i cos(φ i ). ( 3
+ s 2 8 Φ 18,1,0 + e 2 6 Φ 18,1,1 + e 2 8 Φ 18,1,2 + s Φ 18,1,3 +s 4 8 Φ 18,2,0 + e 2 6 Φ 18,2,1 + e 2 8 Φ 18,2,2 + s 2 6 Φ 18,2,3 ) 2Ω 8 -2̟ 6 e 2 6 s 2 8 ( Φ 19,0,0 + e 2 6 Φ 19,0,1 + e 2 8 Φ 19,0,2 + s 2 6 Φ 19,0,3 + s 2 8 Φ 19,1,0 + e 2 6 Φ 19,1,1 + e 2 8 Φ 19,1,2 + s Φ 19,1,3 +s 4 8 Φ 19,2,0 + e 2 6 Φ 19,2,1 + e 2 8 Φ 19,2,2 + s 2 6 Φ 19,2,3 ) 2Ω 8 + ̟ 8 -3̟ 6 e 3 6 e 8 s 2 8 Φ 20,0 + s 2 8 Φ 20,1 + s 4 8 Φ 20,2 3Ω 8 -Ω 6 -2̟ 8 e 2 8 s 6 s 3 8 1 -s 2 8 Φ 21,0 + s 2 8 Φ 21,1 + s 4 8 Φ 21,2 3Ω 8 -Ω 6 -̟ 6 -̟ 8 e 6 e 8 s 6 s 3 8 1 -s 2 8 Φ 22,0 + s 2 8 Φ 22,1 3Ω 8 -Ω 6 -2̟ 6 e 2 6 s 6 s 3 8 1 -s 2 8 Φ 23,0 + s 2 8 Φ 23,1 4Ω 8 -2Ω 6 -2̟ 8 e 2 8 s 2 6 s 4 8 Φ 24,0 + s 2 8 Φ 24,1 + s 4 8 Φ 24,2 4Ω 8 -2Ω 6 -̟ 8 -̟ 6 e 6 e 8 s 2 6 s 4 8 Φ 25,0 + s 2 8 Φ 25,1 + s 4 8 Φ 25,2 + s 6 8 Φ 25,3 4Ω 8 -2Ω 6 -2̟ 6 e 2

# Φ

Semi-major axis function 1

Φ 1,0,0 1 2 b 0 1/2 Φ 1,0,1 1 4 αD + 1 8 α 2 D 2 b 0 1/2 Φ 1,0,2 1 32 α 3 D 3 + 1 128 α 4 D 4 b 0 1/2 Φ 1,0,3 3 16 αD + 9 32 α 2 D 2 + 3 32 α 3 D 3 + 1 128 α 4 D 4 b 0 1/2 Φ 1,0,4 1 8 αD + 7 16 α 2 D 2 + 1 4 α 3 D 3 + 1 32 α 4 D 4 b 0 1/2 Φ 1,0,5 -1 2 αb 1 3/2 Φ 1,0,6 -1 4 α + 1 2 α 2 D + 1 8 α 3 D 2 b 1 3/2 Φ 1,1,0 -1 2 αb 1 3/2 Φ 1,1,1 -1 4 α + 1 2 α 2 D + 1 8 α 3 D 2 b 1 3/2 Φ 1,1,2 -3 32 α 3 D 2 + 1 16 α 4 D 3 + 1 128 α 5 D 4 b 1 3/2 Φ 1,1,3 -3 16 α + 3 4 α 2 D + 9 16 α 3 D 2 + 1 8 α 4 D 3 + 1 128 α 5 D 4 b 1 3/2 Φ 1,1,4 -1 8 α + α 2 D + 19 16 α 3 D 2 + 3 8 α 4 D 3 + 3 32 α 5 D 4 b 1 3/2 Φ 1,1,5 1 2 αb 1 3/2 + 3 4 α 2 5b 0 5/2 + b 2 5/2 Φ 1,1,6 1 4 α + 1 2 α 2 D + 1 8 α 3 D 2 b 1 3/2 + 9 8 α 2 + 9 8 α 3 D + 3 16 α 4 D 2 5b 0 5/2 + b 2 5/2 Φ 1,2,0 3 8 α 2 2b 0 5/2 + b 2 5/2 Φ 1,2,1 9 16 α 2 + 9 16 α 3 D + 3 32 α 4 D 2 2b 0 5/2 + b 2 5/2 Φ 1,2,2 9 64 α 3 D + 27 128 α 4 D 2 + 9 128 α 5 D 3 + 3 512 α 6 D 4 2b 0 5/2 + b 2 5/2 Φ 1,2,3 45 64 α 2 + 45 32 α 3 D + 45 64 α 4 D 2 + 15 128 α 5 D 3 + 3 512 α 6 D 4 2b 0 5/2 + b 2 5/2 Φ 1,2,4 27 32 α 2 + 81 32 α 3 D + 111 64 α 4 D 2 + 3 8 α 5 D 3 + 3 128 α 6 D 4 2b 0 5/2 + b 2 5/2 Φ 1,2,5 -3 4 α 2 6b 0 5/2 + b 2 5/2 + 15 16 α 3 D 29b 1 7/2 + b 3 7/2 Φ 1,2,6 -9 8 α 2 + 9 8 α 3 D + 3 16 α 4 D 2 6b 0 5/2 + b 2 5/2 + 45 16 α 3 + 15 8 α 4 D + 15 64 α 5 D 2 29b 1 7/2 + b 3 7/2 Φ 1,3,0 -5 16 α 3 9b 1 7/2 + b 3 7/2 Φ 1,3,1 -15 16 α 3 + 5 8 α 4 D + 5 64 α 5 D 2 9b 1 7/2 + b 3 7/2 Φ 1,3,2 -15 128 α 3 + 15 32 α 4 D + 45 128 α 5 D 2 + 5 64 α 6 D 3 + 5 64 α 7 D 4 9b 1 7/2 + b 3 7/2 Φ 1,3,3 -225 128 α 3 + 75 32 α 4 D + 225 256 α 5 D 2 + 15 128 α 6 D 3 + 5 1024 α 7 D 4 9b 1 7/2 + b 3 7/2 Φ 1,3,4 -45 16 α 3 + 5α 4 D + 305 128 α 5 D 2 + 25 64 α 6 D 3 + 5 256 α 7 D 4 9b 1 7/2 + b 3 7/2 Φ 1,3,5 15 16 α 3 33b 1 7/2 + b 3 7/2 Φ 1,3,6 45 16 α 3 + 15 8 α 4 D + 15 64 α 5 D 2 33b 1 7/2 + b 3 7/2
TABLE 3.5: Semi-major axis functions for the secular term # 1.

# Φ Semi-major axis function 2 Φ 2,0,0 1 2 -1 2 αD -1 4 α 2 D 2 b 1 1/2 Φ 2,0,1 -1 8 α 2 D 2 + 3 16 α 3 D 3 + 1 32 α 4 D 4 b 1 1/2 Φ 2,0,2 1 8 -1 8 α -11 16 α 2 D 2 -5 16 α 3 D 3 -1 32 α 4 D 4 b 1 1/2 Φ 2,0,3 1 2 α 2 D + 1 8 α 3 D 2 b 0 3/2 + b 2 3/2 Φ 2,1,0 1 2 α 2 D + 1 8 α 3 D 2 b 0 3/2 + b 2 3/2 Φ 2,1,1 1 8 α 2 D + 11 32 α 3 D 2 + 5 32 α 4 D 3 + 1 64 α 5 D 4 b 0 3/2 + b 2 3/2 Φ 2,1,2 3 4 α 2 D + 13 16 α 3 D 2 + 7 32 α 4 D 3 + 1 64 α 5 D 4 b 0 3/2 + b 2 3/2 Φ 2,1,3 -3 4 α 2 + 9 8 α 3 D + 3 16 α 4 D 2 11b 1 5/2 + b 3 5/2 + 1 2 α 2 D + 1 8 α 3 D 2 b 0 3/2 + b 2 3/2 Φ 2,2,0 -3 8 α 2 + 9 16 α 3 D + 3 32 α 4 D 2 5b 1 5/2 + b 3 5/2 Φ 2,2,1 -3 32 α 2 + 39 64 α 3 D + 39 64 α 4 D 2 + 21 128 α 5 D 3 + 3 256 α 6 D 4 5b 1 5/2 + b 3 5/2 Φ 2,2,2 -9 16 α 2 + 57 32 α 3 D + 141 128 α 4 D 2 + 27 128 α 5 D 3 + 3 256 α 6 D 4 5b 1 5/2 + b 3 5/2 Φ 2,2,3 3 4 α 2 + 9 8 α 3 D + 3 16 α 4 D 2 13b 1 5/2 + b 3 5/2 + 75 32 α 3 + 15 8 α 4 D + 15 64 α 5 D 2 29b 0 7/2 + 30b 2 7/2 + b 4 7/2 Φ 2,3,0 25 32 α 3 + 5 8 α 4 D + 5 64 α 5 D 2 9b 0 7/2 + 10b 2 7/2 + b 4 7/2 Φ 2,3,1 75 128 α 3 + 195 128 α 4 D + 235 256 α 5 D 2 + 45 256 α 6 D 3 + 5 512 α 7 D 4 9b 0 7/2 + 10b 2 7/2 + b 4 7/2 Φ 2,3,2 125 64 α 3 + 425 128 α 4 D + 185 128 α 5 D 2 + 55 256 α 6 D 3 + 5 512 α 7 D 4 9b 0 7/2 + 10b 2 7/2 + b 4 7/2 Φ 2,3,3 -75 32 α 3 + 15 8 α 4 D + 15 64 α 5 D 2 33b 0 7/2 + 34b 2 7/2 + b 4 7/2

# Φ

Semi-major axis function 3

Φ 3,0 3 16 -3 16 αD + 3 32 α 2 D 2 + 1 8 α 3 D 3 + 1 64 α 4 D 4 b 2 1/2 Φ 3,1 -15 64 α 3 D 2 + 3 32 α 4 D 3 + 1 128 α 5 D 4 b 1 3/2 + b 3 3/2 Φ 3,2 45 128 α 3 D + 99 256 α 4 D 2 + 3 32 α 5 D 3 + 3 512 α 6 D 4 b 0 5/2 + 4b 2 5/2 + b 4 5/2 Φ 3,

# Φ

Semi-major axis function 4

Φ 4,0 1 16 α 3 D 2 b 1 3/2 Φ 4,1 -1 16 α 3 D 2 b 1 3/2 + 3 16 α 3 D + 3 32 α 4 D 2 b 0 5/2 + 5b 2 5/2 Φ 4,2 15 64 α 3 + 15 32 α 4 D + 15 128 α 5 D 2 17b 1 7/2 + 13b 3 7/2 + 3 16 α 3 D + 3 32 α 4 D 2 b 0 5/2 + 6b 2 5/2 Φ 4,3 -15 64 α 3 + 15 32 α 4 D + 15 128 α 5 D 2 19b 1 7/2 +

# Φ

Semi-major axis function 5

Φ 5,0 -1 2 α 2 D + 1 8 α 3 D 2 b 0 3/2 Φ 5,1 9 2 α 2 + 27 4 α 3 D + 9 8 α 4 D 2 b 1 5/2 + 1 2 α 2 D + 1 8 α 3 D 2 b 0 3/2 Φ 5,2 -21 4 α 2 + 63 8 α 3 D + 21 16 α 4 D 2 b 1 5/2 + 75 16 α 3 + 15 4 α 4 D + 15 32 α 5 D 2 8b 0 7/2 + 7b 2 7/2 Φ 5,3 75 16 α 3 + 15 4 α 4 D + 15 32 α 5 D 2 9b 0 7/2 + 8b 2 7/2 + 315 32 α 4 + 175 32 α 5 D + 35 64 α 6 D 2 57b 1 9/2 + 13b 3 9/2 Φ 5,4 -315 64 α 4 + 175 64 α 5 D + 35 128 α 6 D 2 129b 1 9/2 + 29b 3 9/2
TABLE 3.9: Semi-major axis functions for the secular term # 5.

# Φ

Semi-major axis function 6

Φ 6,0 3 4 α + 1 2 α 2 D + 1 16 α 3 D 2 b 1 3/2 Φ 6,1 -3 4 α + 1 2 α 2 D + 1 16 α 3 D 2 b 1 3/2 + lef t( 15 8 α 2 + 15 16 α 3 D + 3 32 α 4 D 2 5b 0 5/2 + b 2 5/2 Φ 6,2 15 8 α 2 + 15 16 α 3 D + 3 32 α 4 D 2 6b 0 5/2 + b 2 5/2 + 225 64 α 3 + 45 32 α 4 D + 15 128 α 5 D 2 29b 1 7/2 + b 3 7/2 Φ 6,

# Φ

Semi-major axis function 7

Φ 7,0 -α 2 D + 1 4 α 3 D 2 b 0 3/2 Φ 7,1 9 2 α 2 + 27 4 α 3 D + 9 8 α 4 D 2 b 1 5/2 Φ 7
Φ 8,0,0 αb 1 3/2 Φ 8,0,1 1 2 α + α 2 D + 1 4 α 3 D 2 b 1 3/2 Φ 8,1,0 -αb 1 3/2 Φ 8,1,1 -1 2 α + α 2 D + 1 4 α 3 D 2 b 1 3/2 + 9 4 α 2 + 9 4 α 3 D + 3 8 α 4 D 2 2b 0 5/2 + b 2 5/2 Φ 8,2,0 3 2 α 2 2b 0 5/2 + b 2 5/2 Φ 8,2,1 -9 4 α 2 + 9 4 α 3 D + 3 8 α 4 D 2 2b 0 5/2 + b 2 5/2 + 45 8 α 3 + 15 4 α 4 D + 15 32 α 5 D 2 9b 1 7/2 + b 3 7/2 Φ 8,3,0 -15 8 α 3 9b 1 7/2 + b 3 7/2 Φ 8,

# Φ

Semi-major axis function 9

Φ 9,0 -α 2 D + 1 4 α 3 D 2 )b 2 3/2 Φ 9,1 3 2 α 2 + 9 4 α 3 D + 3 8 α 4 D 2 2b 1 5/2 + b 3 5/2 + α 2 D + 1 4 α 3 D 2 b 2 3/2 Φ 9,2 -3 2 α 2 + 9 4 α 3 D + 3 8 α 4 D 2 2b 1 5/2 +

# Φ

Semi-major axis function 10

Φ 10,0 -1 2 α 2 D + 1 8 α 3 D 2 b 0 3/2 Φ 10,1 9 2 α 2 + 27 4 α 3 D + 9 8 α 4 D 2 b 1 5/2 + 1 2 α 2 D + 1 8 α 3 D 2 b 0 3/2 Φ 10
Φ 11,0 3 4 α + 1 2 α 2 D + 1 16 α 3 D 2 b 1 3/2 Φ 11,1 -3 4 α + 1 2 α 2 D + 1 16 α 3 D 2 b 1 3/2 + 15 8 α 2 + 15 16 α 3 D + 3 32 α 4 D 2 5b 0 5/2 + b 2 5/2 Φ 11
Φ 12,0 -3 2 α + α 2 D + 1 8 α 3 D 2 b 1 3/2 Φ 12,1 15 4 α 2 + 15 8 α 3 D + 3 16 α 4 D 2 2b 0 5/2 + b 2 5/2 Φ 12,2 -225 32 α 3 + 45 16 α 4 D + 15 64 α 5 D 2 9b 1 7/2 + b 3 7/2

# Φ

Semi-major axis function 13

Φ 13,0 -9 4 α 2 + 27 8 α 3 D + 9 16 α 4 D 2 b 1 5/2 + 1 2 α 2 D + 1 8 α 3 D 2 b 0 3/2 Φ 13

# Φ

Semi-major axis function 14

Φ 14,0,0 1 2 αb 1 3/2 + 3 4 α 2 b 5/2 0 + 3 2 α 2 b 5/2 2 Φ 14,0,1 1 4 α + 1 2 α 2 D + 1 8 α 3 D 2 b 1 3/2 + 9 8 α 2 + 9 8 α 3 D + 3 16 α 4 D 2 b 0 5/2 + 2b 2 5/2 Φ 14,1,0 3 2 α 2 b 0 5/2 + 9 4 α 2 b 2 5/2 15α 3 b 1 7/2 + 15 4 α 3 b 3 7/2 Φ 14,1,1 Φ 14,0,1 Φ 14,
-3 4 α 2 + 9 8 α 3 D + 3 16 α 4 D 2 b 1 5/2 + 2b 3 5/2 + 1 2 α 2 D + 1 8 α 3 D 2 b 3/2 2 Φ 15,1 - 3 4 α 2 + 9 8 α 3 D + 3 16 α 4 D 2 2b 1 5/2 +
Φ 16,0 -1 32 α 4 D 3 + 1 192 α 5 D 4 b 2 3/2 Φ 16,1 9 64 α 4 D 2 + 5 64 α 5 D 3 + 1 128 α 6 D 4 b 1 5/2 + b 3 5/2 Φ 16
Φ 17,0,0 1 16 α 3 D 2 b 1 3/2 Φ 17,0,1 3 16 α 3 D 2 + 1 8 α 4 D 3 + 1 64 α 5 D 4 b 1 3/2 Φ 17,0,2 5 32 α 3 D 2 + 1 16 α 4 D 3 + 1 192 α 5 D 4 b 1 3/2 Φ 17,0,3 -1 16 α 3 b 1 3/2 + 3 16 α 3 D + 3 32 α 4 D 2 5b 0 5/2 + b 2 5/2 Φ 17,1,0 -3 16 α 3 D + 3 32 α 4 D 2 b 0 5/2 + b 2 5/2 Φ 17,1,1 -9 16 α 3 D + 27 32 α 4 D 2 + 9 32 α 5 D 3 + 3 128 α 6 D 4 b 0 5/2 + b 2 5/2 Φ 17,1,2 -15 32 α 3 D + 33 64 α 4 D 2 + 1 8 α 5 D 3 + 1 128 α 6 D 4 b 0 5/2 + b 2 5/2 Φ 17,1,3 15 32 α 3 + 15 16 α 4 D + 15 64 α 5 D 2 19b 1 7/2 + b 3 7/2 + 3 8 α 3 D + 3 16 α 4 D 2 3b 0 5/2 + b 2 5/2 Φ 17

# Φ

Semi-major axis function 18

Φ 18,0,0 -1 2 α 2 D + 1 8 α 3 D 2 b 0 3/2 Φ 18,0,1 -1 8 α 2 D + 11 32 α 3 D 2 + 5 32 α 4 D 3 + 1 64 α 5 D 4 b 0 3/2 Φ 18,0,2 -3 4 α 2 D + 13 16 α 3 D 2 + 7 32 α 4 D 3 + 1 64 α 5 D 4 b 0 3/2 Φ 18,0,3 9 2 α 2 + 27 4 α 3 D + 9 8 α 4 D 2 b 1 5/2 + 1 2 α 2 D + 1 8 α 3 D 2 b 0 3/2 Φ 18,1,0 3 2 α 2 + 9 4 α 3 D + 3 8 α 4 D 2 b 1 5/2 Φ 18,1,1 3 8 α 2 + 39 16 α 3 D + 39 16 α 4 D 2 + 21 32 α 5 D 3 + 3 64 α 6 D 4 b 1 5/2 Φ 18,1,2 9 4 α 2 + 57 8 α 3 D + 141 32 α 4 D 2 + 27 32 α 5 D 3 + 3 64 α 6 D 4 b 1 5/2 Φ 18,1,3 -6α 2 + 9α 3 D + 3 2 α 4 D 2 b 1 5/2 + 75 8 α 3 + 15 2 α 4 D + 15 16 α 5 D 2 7b 0 7/2 + 3b 2 7/2 Φ 18
Φ 19,0,0 3 4 α + 1 2 α 2 D + 1 16 α 3 D 2 b 1 3/2 Φ 19,0,1 5 32 α 3 D 2 + 1 16 α 4 D 3 + 1 192 α 5 D 4 b 1 3/2 Φ 19,0,2 3 8 α + 3 2 α 2 D + 9 8 α 3 D 2 + 1 4 α 3 D 3 + 1 64 α 5 D 4 b 1 3/2 Φ 19,0,3 -3 4 α + 1 2 α 2 D + 1 16 α 3 D 2 b 1 3/2 + 15 8 α 2 + 15 16 α 3 D + 3 32 α 4 D 2 b 0 5/2 + 5b 2 5/2 Φ 19,1,0 -15 8 α 2 + 15 16 α 3 D + 3 32 α 4 D 2 b 0 5/2 + b 2 5/2 Φ 19,1,1 -15 32 α 3 D + 33 64 α 4 D 2 + 1 8 α 5 D 3 + 1 128 α 6 D 4 b 0 5/2 + b 2 5/2 Φ 19,1,2 -45 16 α 2 + 45 8 α 3 D + 45 16 α 4 D 2 + 15 32 α 5 D 3 + 3 128 α 6 D 4 b 0 5/2 + b 2 5/2 Φ 19,1,3 15 4 α 2 + 45 8 α 3 D + 3 16 α 4 D 2 b 0 5/2 +
Φ 20,0 -5 8 α 2 D + 15 32 α 3 D 2 + 3 32 α 4 D 3 + 1 192 α 5 D 4 b 2 3/2 Φ 20,1 15 16 α 2 + 75 32 α 3 D + 9 8 α 4 D 2 + 11 64 α 5 D 3 + 1 128 α 6 D 4 b 1 5/2 + b 3 5/2 Φ 20

Resonant terms

We define the resonant angle φ = 5λ 8λ 6 . Indirect parts only produce resonant terms. Here for the extended model and as for the simple model, indirect parts only contain one term which are very similar to Equation 3.32 Equation 3.33.

R E = - 3125 384 e 4 8 1 -s 2 8 cos (5λ 8 -λ 6 -4̟ 8 ) (3.39) R I = - 125 384 e 4 8 1 -s 2 8 cos (5λ 8 -λ 6 -4̟ 8 ) (3.40)
The only addition to the simple model is the 1s 2 8 coefficient. Then the 61 terms coming from the direct part are listed in Table 3.36, Table 3.37 and Table 3.38.

R res = 61 i=1 S i cos(φ i ) (3.41) # Argument (φ i ) Amplitude (S i ) 1 φ -4ω 6 e 4 6 Φ 1,0 + s 2 8 Φ 1,1 + s 4 8 2Φ 1,2 2 φ -ω8 -3ω 6 e 8 e 3 6 Φ 2,0 + s 2 8 Φ 2,1 + s 4 8 Φ 2,2 3 φ -2ω 8 -
α 3 D 3 + 1 384 α 4 D 4 b 1 1/2 Φ 5,1 -3125 768 α + 295 96 α 2 D + 85 128 α 3 D 2 + 5 96 α 4 D 3 + 1 768 α 5 D 4 b 0 3/2 + b 2 3/2 Φ 5,2 5485 1024 α 2 + 845 256 α 3 D + 315 512 α 4 D 2 + 11 256 α 5 D 3 + 1 1024 α 6 D 4 5b 1 5/2 + b 3 5/2

# Φ

Semi-major axis function Φ 14,0,3

7 Φ 7,0 -25 4 α + 2α 2 D + 1 8 α 3 D 2 b 3 3/2 Φ 7,1 25 4 α + 2α 2 D + 1 8 α 3 D 2 b 3 3/2 + 99 8 α 2 + 27 8 α 3 D + 3 16 α 4 D 2 b 2 5/2 +
α + 3 4 α 2 D + 1 16 α 3 D 2 b 4s 3/2 Φ 13,0,1 -77 48 α + 23 48 α 2 D -17 64 α 3 D 2 -1 12 α 4 D 3 -1 192 α 5 D 4 b 4 3/2 Φ 13,0,2 -1519 32 α + 251 16 α 2 D -15 64 α 3 D 2 -5 16 α 4 D 3 -1 64 α 5 D 4 b 4 3/2 Φ 13,0,3 - 31 16 α + 3 4 α 2 D + 1 16 α 3 D 2 b 4 3/2 + 129 32 α 2 + 21 16 α 3 D + 3 32 α 4 D 2 5b 3 5/2 + b 5 5/2 Φ 13,1,0 -129 32 α 2 + 21 16 α 3 D + 3 32 α 4 D 2 b 3 5/2 + b 5 5/2 Φ 13,1,1 25 8 α 2 -5 64 α 3 D -99 128 α 4 D 2 -5 32 α 5 D 3 -1 128 α 6 D 4 b 3 5/2 + b 5 5/2 Φ 13,1,2 6063 64 α 2 + 1461 64 α 3 D -225 128 α 4 D 2 -9 16 α 5 D 3 -3 128 α 6 D 4 b 3 5/2 + b 5 5/2 Φ 13
25 4 α + 2α 2 D + 1 8 α 3 D 2 b 3 3/2 + 99 8 α 2 + 27 8 α 3 D + 3 16 α 4 D 2 5b 2 5/2 + b 4 5/2 Φ 14,1,0 99 8 α 2 + 27 8 α 3 D + 3 16 α 4 D 2 b 2 5/2 + b 4 5/2 Φ 14,1,1 -429 64 α 2 -411 64 α 3 D -543 128 α 4 D 2 -39 64 α 5 D 3 -3 128 α 6 D 4 b 3 5/2 + b 4 5/2 Φ 14,1,2 -7497 64 α 2 + 1113 64 α 3 D -543 128 α 4 D 2 -45 64 α 5 D 3 -3 128 α 6 D 4 b 2 5/2 + b 4 5/2 Φ 14,1,3 - 99 4 α 2 + 27 4 α 3 D + 3 8 α 4 D 2 3b 2 5/2 +
α + 5 4 α 2 D + 1 16 α 3 D 2 b 2 3/2 Φ 15,0,1 -85 32 α -95 16 α 2 D -213 64 α 3 D 2 -7 16 α 4 D 3 -1 64 α 5 D 4 b 2 3/2 Φ 15,0,2 -1355 48 α + 145 48 α 2 D -75 64 α 3 D 2 -1 6 α 4 D 3 -1 192 α 5 D 4 b 2 3/2 Φ 15,0,3 - 85 16 α + 5 4 α 2 D + 1 16 α 3 D 2 b 2 3/2 + 315 32 α 2 + 33 16 α 3 D + 3 32 α 4 D 2 5b 1 5/2 + b 3 5/2 Φ 15,1,0 -315 32 α 2 + 33 16 α 3 D + 3 32 α 4 D 2 b 1 5/2 + b 3 5/2 Φ 15,1,1 -315 64 α 2 + 1209 64 α 3 D + 891 128 α 4 D 2 + 3 4 α 5 D 3 + 3 128 α 6 D 4 b 1 5/2 + b 3 5/2 Φ 15,1,2 375 8 α 2 + 65 64 α 3 D -321 128 α 4 D 2 -9 32 α 5 D 3 -1 128 α 6 D 4 b 1 5/2 + b 3 5/2 Φ 15

# Φ

Semi-major axis function 24

Φ 24,0,0 -3 2 α 2 b 3 5/2 Φ 24,0,1 -3 4 α 2 + 9 4 α 3 D + 3 8 α 4 D 2 b 3 5/2 Φ 24,0,2 141 4 α 2 -9 4 α 3 D -3 8 α 4 D 2 b 3 5/2 Φ 24

Comparison between N-Body and Semi-analytic code.

Tests between the N-Body code, the first semi-analytic modelling and the extended one have been made. Direct comparisons show that elements of Titan and Iapetus behave similarly (Figure 3.2 and Figure 3.3). However, semi-analytic models seem to output a different evolution of Iapetus' eccentricity. Figure 3.4 shows a difference in frequency. A first reason for such a mismatch was thought to be the initial conditions used for N-Body codes. Such initial conditions should be fine tuned in order for the model to correctly average the N-Body code. Several initial conditions were tested "by hand" but no improvement was made. Another test starting with a low eccentricity for Iapetus was run (Figure 3.5). As semi-analytic models seem to better approximate the eccentricity in this simulation, the problem persisted. Codes were then thought to be implemented erroneously, as the whole subroutine for the semi-analytic coding was more the 15 000 lines long. Many efforts of checking and rechecking all the terms were made but the difference was still there. Besides, the frequency shift in eccentricity is associated to an increase in the period of the pericenter (Figure 3.6). If an implementation error of the J 2 effect or the Sun's perturbation were to be source of the shift then it would have also influenced the evolution of inclinations, both of which are matching the N-Body very well for both semi-analytic models.

These differences happen for both semi-analytic models which were implemented and verified independently. Therefore, we believe that a certain kind of mis-modelling associated with the type of expansion should be the source of the shift.

Legendre Expansion

In order to check the terms given by Equation C.23, Equation 3.28 and Equation 3.29, direct expansions were used. The direct part can be expanded in terms of Legendre polynomial (Appendix C). From Equation 3.26

1 |r 8 -r 6 | = 1 r 8 ∞ l=0 r 6 r 8 l P l (cos(Ψ)) (3.42)
Then by introducing the semi-major axes, one gets for the direct part

R D = ∞ l=0 α l r 6 a 6 l a 8 r 8 l+1 P l (cos(Ψ)) (3.43)
where r a and a r are given by Equation C.13 and Equation C.14. In order to perform all the calculations, we used Mathematica to expand Equation 3.43 to α 11 and average over the mean longitudes but keeping the resonant angles 5λ 8λ 6 , 10λ 8 -2λ 6 and 15λ 8 -3λ 6 intact. The whole expansion is too long to be exposed in this thesis, so we will explain a few features. We have counted 31 secular terms and 61 resonant terms, the same as in Table 3.4, Table 3.36, Table 3.37 and Table 3.38. For instance, the 27 th term in Table 3 where Φ are functions of the semi-major axis ratio only.

Φ 27,0 = 1 256 α 5 D 3 + 1 1024 α 6 D 4 b 2 5/2 (3.46) Φ 27,1 = - 15 512 α 5 D 2 + 5 256 α 6 D 3 + 5 2048 α 7 D 4 b 1 7/2 + b 3 7/2 (3

Summary

Two different models are built in this chapter. A first and simple one was constructed, implementing terms already given in [START_REF] Murray | Solar System Dynamics[END_REF] and the second extended the first one by changing the definition of the order of the resonance. Here, as the inclination of Iapetus is not negligible, s 8 was not constrained. However for each term, we truncated the series in s 8 when the difference between two consecutive orders of s 8 was under 10%. Eventually solar terms in λ ⊙ were added and also several terms in 10λ 8 -2λ 6 in later models.

Chapter 4

The Laplace Plane

The benefit of having semi-analytic equations is to be able to study the system in a more synthetic way. Now we have developed approximated equations for the motion of satellites, we can use them in order to compute analytic solutions of the equations of motion. In this chapter, we study the notion of the local Laplace Plane of a satellite. We have seen in the last chapter that the disturbing function associated with the Sun and averaged over the longitude of the Sun and the satellite can be written

R ⊙ =n 2 ⊙ a 2 s 1 32 8 + 12e 2 ⊙ + 15e 4 ⊙ + 6e 2 s 2 + 3e 2 ⊙ 1 -6s 2 s + 6s 4 s 1 -6s 2 ⊙ + 6s 4 ⊙ + cos(2Ω s -2Ω ⊙ ) 3 8 8 + 12e 2 ⊙ + 15e 4 ⊙ + 6e 2 s 2 + 3e 2 ⊙ s 2 s 1 -s 2 s s 2 ⊙ 1 -s 2 ⊙ + cos(Ω s -Ω ⊙ ) 3 8 8 + 12e 2 ⊙ + 15e 4 ⊙ + 6e 2 s 2 + 3e 2 ⊙ s s 1 -s 2 s 1 -2s 2 s s ⊙ 1 -s 2 ⊙ 1 -2s 2 ⊙ + cos(2ω s -2Ω ⊙ ) 15 8 e 2 s 2 + 3e 2 ⊙ 1 -s 2 s 2 s 2 ⊙ 1 -s 2 ⊙ + cos(4Ω s -2Ω ⊙ -2ω s ) 15 8 e 2 s 2 + 3e 2 ⊙ s 4 s s 2 ⊙ 1 -s 2 ⊙ + cos(3Ω s -Ω ⊙ -2ω s ) 15 4 e 2 s 2 + 3e 2 ⊙ s 3 s 1 -s 2 s s ⊙ 1 -s 2 ⊙ 1 -2s 2 ⊙ + cos(Ω s + Ω ⊙ -2ω s ) - 15 4 e 2 s 2 + 3e 2 ⊙ s s 1 -s 2 s 3/2 s ⊙ 1 -s 2 ⊙ 1 -2s 2 ⊙ + cos(2ω s -2Ω s ) 15 8 e 2 s 2 + 3e 2 ⊙ s 2 s 1 -s 2 s 1 -6s 2 ⊙ + 6s 4 ⊙
and we recall the average disturbing function of the flattening effect

R J 2 = G(M + m s )J 2 R 2 p a 3 (1 -e 2 ) 3/2 1 2 -3s 2 s 1 -s 2 s .
We form

R LP = R ⊙ + R J 2 (4.1)
Here the value of J 2 is J ′ 2,I given in Equation 2.24. We then plug R LP into the Lagrange Planetary Equations (Equation B.16)

di dt = - 1 na 2 E sin(i) ∂R LP ∂Ω - tan i 2 na 2 E ∂R LP ∂̟ dΩ dt = 1 na 2 E sin(i) ∂R LP ∂i .
where

ǫ J 2 = Gm S 1 + m i m S 16J 2 R 2 p 1 -e 2 i 3/2 (4.5)
and

ǫ ⊙ = Gm ⊙ 1 + m S m ⊙ 8a 2 i α 3 i,⊙ × 1 + 3 2 e 2 ⊙ + 15 8 e 4 ⊙ + 1 4 e 2 i 2 + 3e 2 ⊙ (3 -5 cos(2ω i )) (4.6) so that for m S m ⊙ << 1 and m i m S << 1 ǫ J 2 ǫ ⊙ ≈ 2 m S m ⊙ J 2 R 2 p a 3 ⊙ a 5 i F (e ⊙ , e i ).
(4.7)

F (e ⊙ , e i ) contains the dependency on both eccentricities

F (e ⊙ , e i ) = 1 1 -e 2 i 3/2 × 1 1 + 3 2 e 2 ⊙ + 15 8 e 4 ⊙ + 1 4 e 2 i 2 + 3e 2 ⊙ (3 -5 cos(2ω i ))
.

(4.8)

The first equation of 4.4 can be rewritten

tan (2i LP ) = sin(2i ⊙ ) ǫ J 2 ǫ ⊙ + cos (2i ⊙ )
.

(4.9)

Chapter 5

Sorting out results

In this chapter, we give the results of numerous numerical simulations of the 5:1 mean motion resonance crossing between Titan and Iapetus. Outcomes of the resonance are very different in nature but nevertheless we can separate them in three different ways. But first, we give the initial conditions chosen for the simulations.

Preliminary assessment

A dynamical system consists of a set of Ordinary Differential Equations, which in our case, were derived in the previous chapters, and a set of initial conditions. The latter were a source of error during the beginning of this work, as the initial ascending node of Iapetus was taken originally random. In chapter 4, we have shown that the orbital plane of Iapetus is forced around its Laplace plane, which has a 15 degree inclination with respect to the equator of Saturn and an ascending node equal to that of the Sun. In the numerical integrations, we have chosen to set the Sun's ascending node to 0. Therefore the tilt to the Laplace plane is computed using cos(J) = cos(i) cos(i LP ) + sin(i) sin(i LP ) cos(Ω) (5.1)

where i LP is the Laplace plane inclination set to 15.226 degrees.

Initial conditions

The equations of motion written in the previous chapter govern the orbital evolution of Titan and Iapetus through the resonance. Initial conditions were chosen so that Titan would migrate through the 5:1 mean motion resonance with Iapetus. We have chosen to perform 100 runs using different initial values for the mean anomaly of Iapetus so that we can study the resonance crossing statistically.

There were some concerns about the initial value of the semi-major axis for Titan. The semi-major axis ratio between Titan and Iapetus 1 was set to a value greater than 5 so that the 5:1 mean motion resonance would be crossed during the simulation. On top of this, we made sure that the initial ratio was far enough from the resonance so that all sub-resonances would be reached during a simulation. Choosing a value too high would have made the simulation starting directly with a resonant interaction between Titan and Iapetus. On the other hand, choosing a low 1 which we denote as α in this thesis. value would have made simulations much longer. Also, there is a fundamental difference between the N-Body and semi-analytic model concerning Iapetus' semimajor axis. For the N-body code, two different initial mean anomalies would give two different behaviours for the semi-major axis of Iapetus, as shown in Figure 5.1. The average value is different from a initial M to another, making the resonance longer to reach for simulation showing higher semi-major axes. This does not appear in the semi-analytic code as the semi-major axis only gets greatly disturbed when the resonance is encountered. Therefore, to make sure simulations start before the resonance, the value 8.131 × 10 -3 AU was chosen for the N-Body code (Table 5.1) whereas 8.1405 × 10 -3 AU was good enough for the semi-analytic code. We assume for most of our simulations that the eccentricity of Titan has not changed during the epoch of the resonance crossing, therefore it was set to 0.028, today's mean value. One can also give arguments for an evolution of Titan's average eccentricity with time. As explained in chapter 1, the tidal interaction between the planet and the satellite yields change in the satellite's eccentricity2 but because of the difficulty of making assessment on tidal parameters for bodies in general, we have assumed that the eccentricity did not change. We note that Fuller, Luan, and Quataert, 2016 does not argue about any change of eccentricity but emphasizes the semi-major axis time variation and the implication for the value of effective quality factors 3 .

Titan is also tied to its local Laplace plane chapter 4 which makes its maximum inclination around 1 degree. We have chosen this value as an initial condition (the initial ascending node is then 0 according to Equation 5.1).

Chaos

The crossing the 5:1 mean motion resonance is chaotic. Many resonant terms given in Table 3.36, Table 3.37 and Table 3.38 will be reached for the same value of the semi-major axis ratio as shown in Figure 5.2. For a general resonant term, one has 5n 8n 6 + j 3 ̟8 + j 4 ̟6 + j 5 Ω8 + j 6 Ω6 = 0

(5.2)

where the semi-major axis ratio, α, is

α = n 8 n 6 2/3
(5.3)

Here because of the motion of the orbital planes of Titan and Iapetus around their local Laplace plane, the rate of their ascending node is zero on average Ω6 = 0 (5.4) Ω8 = 0.

(5.5)

Therefore the value of α for which a resonant term oscillates is

α = 1 5 2/3 1 -j 3 ̟8 n 6 -j 4 ̟6 n 6 2/3
(5.6) Using Equation 5.6 for every term in Table 3.36, Table 3.37 and Table 3.38, one obtains the following positions of all sub-resonances

Simulations

Three different types of outcome of been produced :

• Capture and release from the resonance (Figure 5.3). Simple "kicks" or rapid resonant interactions are also observed in this type.

• Ejection of Iapetus Figure 5.16

• Smooth evolution of either the eccentricity or the inclination Figure 5.4.

The resonance acts on both the eccentricity and the inclination of Iapetus as we expected from the analytical expansions done in chapter 3. The outcome varies from one simulation to another due to chaos. The most interesting results are those showing capture into resonance and then release. In those cases both the eccentricity and the inclination behaved chaotically and their value after the release are unpredictable due to the stochastic nature of the dynamics. Figure 5.3 shows one of these simulations where the final elements are close to the values of Iapetus' eccentricity and tilt we see today. In this example we clearly distinguish a regular evolution before the resonance, a chaotic resonance crossing and finally regular dynamics after the resonance. In this specific example, the eccentricity grows up to 0.12 before decreasing and finally settling to around 0.05 and the tilt reached 8 degrees before settling at 4.5 degrees after the two bodies for released from resonance.

The second important result coming from simulations is the simple ejection of Iapetus out of the system of Saturn. In those cases, the resonance was strong enough to raise Iapetus eccentricity up to a certain threshold (around 0.25) which led to an automatic ejection (see Figure 5.16).

Finally, the third behaviour corresponds to a "smooth" evolution of the eccentricity or the inclination as depicted in Figure 5.4. We note that elements do not increase simultaneously but separately. For instance while the eccentricity increases, the inclination stays constant and vice-versa. A first guess makes us think that Iapetus was captured in one of the sub-resonances outlined in chapter 3 but a study of 19 resonant angles (chapter 3) did not show any of them entering into libration when the resonance has been reached. Arguments susceptible of leading to such behaviour would be the eccentricity-like arguments, for instance the angle 5λ 8λ 6 -4̟ 8 (Figure 5.5) slowly circulates as other angles. Therefore such evolution should come from the slow variations of several arguments acting together as for any chaotic evolution. The study is meant to evaluate the action of Titan on both the eccentricity and the inclination of Iapetus. As these simulations show the variation of only one element at a time, no further analyses were done.

N-Body code Q=20, Q=50, Q=75

As we started our study of the 5:1 mean motion resonance between Titan and Iapetus, we were first interested in low values of Q which are compatible with the theory outlined in Fuller, Luan, and Quataert, 2016. As for most of the work during this thesis, 100 simulations were done for every value of Q chosen. For Q = 20, Q = 50 and Q = 75, we observe that a good number of simulations show a "kick" in both eccentricity and inclination Figure 5.9. Others show a chaotic evolution of both elements which is usually short in time. Outcomes show a growth in eccentricity of a few percent (up to over 0.05 for a few simulations done with Q = 50 or Q = 75). On the other hand, Iapetus' orbital plane stays mainly close to its Laplace equilibrium. The tilt, which still gets excited by the resonance, rarely grows over 1 degree. As we can see with Figure 5.6, Figure 5.7 and Figure 5.8 high excitations of the eccentricity happens primarily for higher values of Q. Here simulations done with Q = 75 have post-resonance eccentricities in average higher than simulations done with Q = 50 or Q = 20. The same trend is true the inclination also. Dynamically, these first simulations suggest that a rapid crossing of the resonance is insufficient to excite Iapetus orbit, specially in the case of its inclination. FIGURE 5.9: Example of a simulation showing a fast and sudden change in the eccentricity and the inclination of Iapetus during the resonance. The crossing is obvious when the eccentricity "jumps" and settles to a higher value. The inclination shows a quick variation but comes down to its original value right after. The evolution is regular before and after this event. This simulation was done with Q = 20. In general, Q < 100 happens to be insufficient for the tilt to reach a value close to today's value (8 degrees). However, these simulations show that the eccentricity of Iapetus can be easily excited by the resonance. Also importantly, we have not seen any ejection or smooth evolution in this range of values. Q=100, Q=200, Q=400, Q=600 and Q=1500 Subsequently, our work brought us to run simulations with higher values of Q. Spanning from Q = 100 to Q = 1500, we have done 100 simulations for each value of Q, except for Q = 100, for which a total of 400 simulations were done. Contrary to lower values we have discussed in the previous subsection, chaos dominates the evolution of the orbit of Iapetus during the resonance. Generally, we have found that the resonant interaction between the two satellites is longer, giving time for the eccentricity and the inclination to reach higher values than in the previous simulations.

FIGURE 5.11: Post-resonant eccentricities and inclinations using Q=100. Out of 400 simulations, 5 end with the ejection of Iapetus, 1 shows a smooth evolution of its inclination and 5 show a smooth evolution of its eccentricity.

FIGURE 5.12: Post-resonant eccentricities and inclinations using Q=200. Out of 100 simulations, 3 end with the ejection of Iapetus and 15 show a smooth evolution of its eccentricity. Besides having unpredictable post-resonant elements, chaos makes the time during which both bodies are in resonance also unpredictable. For a given initial condition, one can not predict a priori how long the resonance will last. Usually, as for simulations done with a lower Q, Iapetus is released from the resonance with an excited eccentricity and inclination but the long resonant interaction between the two bodies leads to sometimes excessive eccentricities. We find that if the eccentricity of Iapetus exceeds a threshold of around 0.20 or 0.25 then Iapetus will get ejected from the system of Saturn later on in the simulation (??). When doing an overview of all simulations, it appears that the number of ejections is greater if Titan migrates slowly through the resonance (for high value of Q in other words). The same trend happens for smooth evolutions. In the first place, chaos makes the eccentricity vary rapidly on a relatively short time scale, then reaches 0.2 in a monotonic way. Once it has reached that value 20 millions years after the start of the simulation, the eccentricity seems to settle for a while before being excited again and growing drastically towards ejection (e>1). The settling of the eccentricity at around 0.2 illustrates the action of sub-resonances during the process.

The ejection and the smooth evolution scenario are unwanted and the number of simulations showing these behaviours increases with the quality factor. The rest of the simulations shows a release of the resonance followed by a regular dynamical evolution. In that range of quality factor many simulations have shown interesting postresonant elements. We recall that we are assessing the ability of Titan to excite Iapetus' orbit to around 0.03 in eccentricity and up to 8 degrees in the tilt to the Laplace plane with the 5:1 mean motion resonance. Figure 5.17 and Figure 5.18 show two simulations which are in good agreement with these target elements. FIGURE 5.17: Evolution of the eccentricity and the tilt with Q=400. This simulation is an example of the great range of values which can reach the eccentricity and the inclination during the resonance. Although a high value of 0.18 is reached, chaotic perturbations bring the eccentricity down rapidly before Iapetus leaves the resonance with a low 0.015 eccentricity. On the other hand the tilt reaches 10 degrees but does not undergo any drastic diminishing. Iapetus comes out of the resonance with a tilt to its local Laplace plane of about 7.5 degrees. FIGURE 5.18: Example of a simulation showing post-resonance eccentricity and inclination compatible with Iapetus' actual orbit. Here, the resonance is crossed 90 millions years after the start of the simulation. As Iapetus evolve chaotically during the crossing, its eccentricity raises over 0.2, while the tilt also increases to over 10 degrees. On such an orbit Iapetus is expected to get ejected, but at the end, before getting out of the resonance, the eccentricity decreases to 0.05 while the tilt stays well over the value we observe today. This simulation was done with Q=1500

Other "less successful" simulations show inclinations reaching reasonably good values, (Figure 5.3 is a good example of such simulations). The inclination settles at around 4.5 degrees after having reached a maximum of 8 degrees. The eccentricity after the release of the resonance is 0.05. This particular resonance agrees to a lesser extent with what we are looking for. Although post-resonant elements do not correspond to target elements, this simulation, also, shows that Titan, through the resonance, can increase Iapetus inclination to 8 degrees and beyond (as for Figure 5.17 and Figure 5.18). In the same vein, Figure 5.19 shows a 4.5 degrees increase in the tilt while keeping the eccentricity low. Iapetus comes out of the resonance with a 0.023 eccentricity but the post-resonant eccentricity value may have been between the minimum and the maximum reached during the chaotic evolution. An overview of post-resonant elements is given in 200,400,600,1500. i * is the tilt to the local Laplace plane of Iapetus and e its eccentricity after the satellite has been released from the resonance.

Semi-analytic code

The development of a semi-analytic code seemed necessary in order to run simulations with high values of Q. As short period variations were cancelled from the system, using the semi-analytic code in chapter 3 enabled us to make faster simulations to study the resonance with Q greater than 2000. Here, as for the N-Body code, we show the simulations done with the semi-analytic modelling. Histograms in Figure 5.21 show the post-resonant elements for simulations with Q=20. As for the N-Body code, the inclination was not highly excited by the resonance, whereas the eccentricity sometimes had values over 0.05. Many of simulations done with Q = 20 shows just a "kick", just as Figure 5.9, but others show a capture and a chaotic evolution of the eccentricity and the inclination.

For higher values, outcomes are similar to the ones obtained with the N-Body code. Chaos dominates the evolution of the eccentricity and the inclination and after the release values are comparable to the ones obtained with the N-Body code. However we have obtained only one smooth evolution (which appeared at Q=18000) and a greater number of ejections. As we analyse the overall post-resonant elements, it seems that the semi-analytic code outputs a lesser number of high inclinations than the N-Body code (Table 5.3 and Table 5.4). However, both methods do not have the same definition of an ejection. For the N-Body code, one calls an ejection when the eccentricity reaches 1 and beyond (basic definition of an open orbit in orbital dynamics). This definition can not hold for the semi-analytic approach, as the expansion of the disturbing function (chapter 3) is truncated to the fourth power in eccentricity. Therefore the equations of motion are not adequate when the eccentricity of Iapetus reaches high values. In our analysis software for the semi-analytic code, an ejection was considered when the eccentricity reached 0.25 because no simulations done with the N-Body code showed a resonance release with an eccentricity being over this limit. For the N-Body code, an eccentricity of 0.25 was a threshold over which Iapetus' was doomed to get ejected sooner or later.

Outputs

Number of simulations i * > 5 • + e < 0.05 to 200, 400, 600, 800, 1200, 1600, 2000, 5000, 18000 i * is the tilt to the local Laplace plane of Iapetus and e its eccentricity after the satellite has been released from the resonance.

0 i * > 5 • 3 i * > 4 • 17 i * >

Trends

An overall analysis of all simulations shows that the number of ejections increase with the quality factor. This is true for both sets of simulations (Figure 5.35). For the N-Body, more than a quarter of simulations done with Q=1500 shows an ejection of Iapetus and one can extrapolate that the number of ejections would get higher for higher values of Q. The same trend appears for the semi-analytic approach but with more ejections. Although both methods have equivalent outputs, the higher number of ejections given by the semi-analytic code was a source of concern. We can surely mention the different periods in pericentre motions done in tests for Figure 3.6 to be the source of such differences (due probably to errors in the implementation). Another reason could be the truncation of the mutual influence of the satellites to the fourth order in eccentricity, which could lead in a easier way to an ejection when the eccentricity of Iapetus gets relatively high. Also, as the semi-analytical code initially approximates the positions of the moons (Equation C.12), one can understand that probably a big number of terms using high order expansions in eccentricity are necessary for the model to approximate the evolution of mean longitudes, and subsequently the resonance. Terms neglected because they are small could become important when integrated over time. And the effect of chaos can only amplify the discrepancy between N-Body and semi-analytic modelling. Following [START_REF] Batygin | Analytical treatment of planetary resonances[END_REF], we understand that an analytic or semi-analytic model can not be considered at the same level as a direct N-Body simulation, but can only bring theoretical insights to the dynamical problem. Another important difference is the lack of simulations showing smooth evolutions of Iapetus' elements. The same reasons as before can apply but it could also be a consequence of the initial conditions which were chosen. Except for one simulation done with Q = 18 000, the semi-analytic program does not output any smooth evolution simulations. On the other hand, the N-Body code produces an increases number of these and the trend is similar to the number of ejections and it increases with Q (Table 5. We observe that the number of ejections using the semi-analytic model is greater than with the N-body. The reason is probably that the rates of mean longitudes are not properly modelled when we average out all the terms except the resonant ones and such biases appear during the resonance capture where mean longitudes play an important role.

Summary and conclusion

The important result which comes out of those simulations and their subsequent analysis is the constraint of the quality factor at the frequency of Titan, which determines the amount of energy dissipation due to tides. An upper value is determined by noticing that simulations showing ejections and smooth evolutions dominate the outcome of the resonance if Q is set to a value greater than around 2000. On the other hand, simulations done with a small quality factor show that chaotic evolution does not dominate, and the inclination of Iapetus can not reach any values close to the ones we observe today. Therefore, Titan could have possibly excited Iapetus' eccentricity and inclination with a past 5:1 mean motion resonance only if Q was between roughly 100 and 2000.

Chapter 6

Additional simulations and modelling 6.1 Additional simulations

First model results

When introducing the semi-analytic modelling in chapter 3, we first decided to implement secular and resonant terms already exposed in the appendixes of [START_REF] Murray | Solar System Dynamics[END_REF] As for the N-Body and the extended semi-analytic model, we have simulated the resonance crossing using several values of the quality factor.

Starting with a tilt

As the number of high tilts reached during simulation was low, we decided to start Iapetus orbital plane outside the Laplace equilibrium by starting with a few degrees of tilt. Results (Figure 6.1) show that the final inclination does not reach higher values in general. The higher the initial tilt, the lesser Iapetus' plane gets excited. Simulations done with an 8 degree tilt (today's value) show that the inclination is not excited (also shown in [START_REF] Ćuk | Titan-Hyperion Resonance and the Tidal Q of Saturn[END_REF].

Adding 10:2 terms

The difference in the number of ejections between semi-analytic and N-Body codes brought us to consider additional resonant terms. We have added to our software terms in 10λ 8 -2λ 6 . We have implemented 31 out of 61 possible terms and saw no additional differences in the outcomes. After the release of the resonance, the eccentricity was between 0 and 0.15 and the tilt climbed to a few degrees, sometimes reaching high values.

Quadruple precision

A run was done using quadruple precision numbers in order the check the sensitivity of the outcome to numerical errors (truncation, loss of precision ...). We have used the extended analytical model and results show no differences.

Axial precession of Saturn

An important dynamical feature which comes into play when the time scale considered is large is the displacement of the rotational axis of the planet. Generally, due to the action of the Sun and its satellites on its oblate figure, a planet polar axis undergoes periodic changes with a large range of frequencies. On the longest time scale, the main polar modification is a constant drift around the orbit normal of the Sun (or more generally around the invariable pole [START_REF] French | Geometry of the Saturn system from the 3 July 1989 occultation of 28 SGR and Voyager observations[END_REF] if orbital plane modifications are considered) called axial precession. At higher frequency, nutation appears as a polar wobble around the main obliquity value. For the Earth, the axial precession was first discovered by Hipparcus in the 2nd century B.C. and generally a first theoretical explanation was given Newton in his Principia. The earth spin axis makes a complete circle on the sky in around 26 000 years.

The reason for the axial precession is the torque of the Sun acting on the planet's oblate shape [START_REF] Bills | Free and forced obliquities of the Galilean satellites of Jupiter[END_REF])

T = 3Gm ⊙ r 3 (u.I × u) (6.1)
where G is the gravitational constant, m ⊙ the mass of the Sun, r Saturn's distance to the Sun, u is a unit vector pointing towards the Sun and I is the inertia tensor of Saturn. After averaging the torque over the rotation period and the orbital period of the body the precession equation becomes

ds dt = 3 2 n 2 σ (α * (n.s) + β * ) (s × n) (6.2)
where s and n are along the orbit pole and spin pole. In the case of Saturn, the time scale is much longer due to the larger mass of the planet and the larger distance from the Sun. As with the Earth, Saturn's satellites are believed to play a major role in its axial precession. However, the fact that Saturn is mostly composed of gas makes a theory for Saturn's spin axis difficult. Thanks to the gravity science conducted with the Cassini spacecraft, the spin axis rate of change has been measured with a good accuracy. [START_REF] French | Noncircular features in Saturn's rings IV: Absolute radius scale and Saturn's pole direction[END_REF] have measured a precession rate on the sky of ṅp = 0.207 ± 0.006 arcseconds per year using essentially Cassini orbiter data, but also previous measurements from Voyager.

Using a vectorial notation, the polar spin axis n p can be calculated in time using a time dependant rotation matrix

n p = R(t)   0 0 1   (6.3)
with initially R(0) = Id, the identity matrix, as initially, the reference frame is Saturn's centre with the z axis pointing towards the North Pole. In general, the rotation matrix of angle θ about a unit vector of direction u is

R θ ( u) =   cos(θ) + u x (1 -cos(θ)) u x u y (1 -cos(θ)) -u z sin(θ) u x u z (1 -cos(θ)) + u y sin(θ) u x u y (1 -cos(θ)) + u z sin(θ) cos(θ) + u 2 y (1 -cos(θ)) u y u z (1 -cos(θ)) -u x sin(θ) u x u z (1 -cos(θ)) -u y sin(θ) u y u z (1 -cos(θ)) + u x sin(θ) cos(θ) + u 2 z (1 -cos(θ))
  (6.4) In our case, the unit vector is the normal to the Sun's orbit, n ⊙ , we compute it using

n ⊙ = r ⊙ × v ⊙ |r ⊙ × v ⊙ | (6.5)
We recall in general that

r =    x = r [cos(Ω) cos(ω + f ) -sin(Ω) sin(ω + f ) cos(i)] y = r [sin(Ω) cos(ω + f ) + cos(Ω) sin(ω + f ) cos(i)] z = r [sin(ω + f ) sin(i)] (6.6)
where r is the distance to the Saturn, which for a conic trajectory is

r = a 1 -e 2 1 + e cos(f ) (6.7) and v =      ẋ = -na √ 1-e 2 {[sin(ω + f ) + e sin(ω)] cos(Ω) + [cos(ω + f ) + e cos(ω)] sin(Ω) cos(i)} ẏ = -na √ 1-e 2 {[sin(ω + f ) + e sin(ω)] sin(Ω) -[cos(ω + f ) + e cos(ω)] cos(Ω) cos(i)} ż = na √ 1-e 2 [cos(ω + f ) + e cos(ω)] sin(i)
(6.8) For us here, we assume that the Sun is on an unperturbed orbit, so its orbital plane stays fixed in time. The normal to its orbit is

n ⊙ =   sin(Ω ⊙ ) sin(i ⊙ ) -cos(Ω ⊙ ) sin(i ⊙ ) cos(i ⊙ )   (6.9)
Or by setting the ascending node to the origin, we have a simpler expression

n ⊙ =   0 -sin(i ⊙ ) cos(i ⊙ )   . (6.10)
Then plugging that into the rotational matrix, one obtains

R θ (n ⊙ ) =   cos(θ) -cos(i) sin(θ) -sin(i) sin(θ) cos(i) sin(θ) cos(θ) + sin(i) 2 (1 -cos(θ)) -sin(i) cos(i)(1 -cos(θ)) sin(i) sin(θ) -sin(i) cos(i)(1 -cos(θ)) cos(θ) + cos(i) 2 (1 -cos(θ))   .
(6.11) The angle θ is related to the time t through θ = ω p t (6.12)

where ω p is the frequency of the precession. Hence the orientation of Saturn's pole is given by

n p =   sin(i ⊙ ) sin(ω p t) -sin(i ⊙ ) cos(i ⊙ ) (1 -cos(ω p t)) cos(ω p t) + cos 2 (i ⊙ ) (1 -cos(ω p t))   (6.13)
We have implemented a modified J 2 force which takes into account the change in Saturn's pole direction. Results show that the orbital planes of Titan and Iapetus both follow the precession of Saturn (Figure 6.3) This result is similar to the dynamics explained in Goldreich, 1965b, but it was not believed a priori that it would apply trivially to our case. If Titan's orbital plane was to precess at a different frequency, then the Laplace plane of Iapetus would be changing in time.

However, the fact that both orbits follow the precessing equator of the planet make the dynamics similar to the one with no precession. Therefore it was adequate not to consider the precession of the planet in our simulations in general.

Hyperion and Jupiter

Initially, Hyperion and Jupiter were added to numerical models. But we found that Hyperion is not influenced by the mean motion resonance between Titan and Iapetus.

We have left aside the gravitational effect of the outer planets, but it is worth asking if it is a fair negligence. Jupiter evolves close to a 5:2 mean motion commensurability with Saturn and plays a major role in its secular motion around the Sun.

One of the main aspects of this perturbation is a secular change of Saturn's eccentricity with a period of around 50 000 years [START_REF] Murray | Solar System Dynamics[END_REF]. If we go back to equations 16 we see that the Laplace Plane equilibrium has a dependency on the planet's eccentricity (alias the Sun's eccentricity in our model). But although the Laplace Plane equilibrium will change at the same frequency as the Sun's eccentricity, the free inclination will still preserve a constant tilt to it.

Conclusion

In this thesis, we have run thousand of simulations in order to study a past 5:1 mean motion resonance between Titan and Iapetus. This study is in the continuation of results obtained in [START_REF] Lainey | Strong Tidal Dissipation in Saturn and Constraints on Enceladus' Thermal State from Astrometry[END_REF][START_REF] Lainey | New constraints on Saturn's interior from Cassini astrometric data[END_REF], where the dissipation inside Saturn was revealed to be high. Consequentially, saturnian satellites are migrating fast and on top of this, the quality factor of the planet Q is a frequency dependant tidal parameter. At the frequency of Enceladus, Tethys and Dione, its value was calculated to be around 1600, whereas for Rhea, Q = 300 fits better the observational data. The theory given in [START_REF] Fuller | Resonance locking as the source of rapid tidal migration in the Jupiter and Saturn moon systems[END_REF] gives values of Q in agreement with the values given by astrometric results and predicts a quality factor of 20 at the frequency of Titan. The condition for a 5:1 mean motion resonance between Titan and Iapetus to happen is that Q < 14000 at the frequency of Titan, therefore Fuller's theory predicts a crossing not long ago. As satellites are migrating rapidly, they should be young, and therefore the dynamical history of the saturnian system has to be reconsidered. [START_REF] Charnoz | Accretion of Saturns mid-sized moons during the viscous spreading of young massive rings: Solving the paradox of silicate-poor rings versus silicate-rich moons[END_REF] shows that creating the icy moons from massive rings accounts for their current positions and masses. However Titan should have played an important role in the early dynamics of the system, therefore we assumed during this thesis several values of Q, spanning from 20 to 18 000. Amongst tidal considerations, the study of a resonance with Titan was brought forward because of the unexplained nature of the orbit of Iapetus. The origins of its eccentricity of 0.03 and the tilt to its local Laplace plane of 8 degrees were still unclear despite the existence of theories suggesting that the tilt may be the natural consequence of the shape of the ancient circumplanetary disk [START_REF] Ward | Orbital inclination of Iapetus and the rotation of the Laplacian plane[END_REF] or the result of the gravitaional pull of a fifth planet during the early stages of the solar system [START_REF] Nesvorný | Excitation of the Orbital Inclination of Iapetus during Planetary Encounters[END_REF].

We were assessing the possibility that a past 5:1 mean motion resonance is responsible for Iapetus' eccentricity and tilt.

During this thesis we have implemented a N-Body code which includes a tidal force which makes Titan migrate through the resonance. Then, because of the time consuming nature of such methods, semi-analytic models were implemented, the most advanced one taking into account the non-negligible inclination of the orbit of Iapetus (chapter 3).

In chapter 4, the elements of the Laplace plane were then computed for a matter of clarification.

Finally the outcome of the resonance crossing were given in chapter 5. 100 simulations were produced for each value of Q because of the chaotic nature of the dynamics during the crossing. A clear trend appears : Titan is more likely to eject Iapetus out of the system if dissipation in the planet is low. Simulations show that if Q is set to be over 2000, then the fate of Iapetus is most likely an ejection of the system. At this stage, the N-Body code has also shown captures which lead to smooth growths of the eccentricity or the inclination. However both scenarios are unwanted and therefore Q is unlikely to be over 2000 at the frequency of Titan. On the other hand Q < 100 leads to weak resonance interaction between the two bodies as Titan rushes through the resonance without perturbing Iapetus much. It is then roughly between Q = 100 and Q = 2000 that simulations show the best agreement with Iapetus' orbit today. With those values, Iapetus gets released from the resonance with an eccentricity of a few percent (spanning 0 to 0.15) and the inclination also gets excited. We found that the eccentricity is easy to excite whereas the tilt to the Laplace plane needs a long resonant and chaotic interaction to reach values above 5 degrees. A low number of those simulations show a compatibility with Iapetus' current element. Out of 800 simulations 2 have shown a post resonance eccentricity of less than 0.05 with a tilt over 5 degrees. 11 show an tilt over 5 degrees but with a high eccentricity. Despite this low probability, the existence of such dynamical trajectories makes plausible the scenario of a recent excitation of Iapetus' elements by a 5:1 mean motion resonance with Titan, keeping in mind that a condition for it to happen is that 100 < Q < 2000 (6.14) at the frequency of Titan. These values indicate a high energy dissipation in Saturn and are in agreement with observations [START_REF] Lainey | New constraints on Saturn's interior from Cassini astrometric data[END_REF]. At this rate, the resonance would have been reached 40 to 800 millions years ago. This work was the subject of an article accepted in Astronomy and Astrophysics [START_REF] Polycarpe | Strong tidal energy dissipation in Saturn at Titan's frequency as an explanation for Iapetus orbit[END_REF] and a poster prize at CELMEC 2017 (http://adams.dm. unipi.it/~simca/celmecVII/prizes.html)

Further prospects

A subsequent work which can be undertaken in the future would concern the orientation of Iapetus' spin axis. We believe that if the chaotic perturbation with Titan happened not long ago, then there might be a dynamical trace of a chaotic orbital evolution of Iapetus on its spin axis. [START_REF] Davies | The control network of Iapetus[END_REF] show that the obliquity of Iapetus was low (below 2 degrees but with uncertainties) using data from the two Voyager spacecrafts. Also [START_REF] Bills | Forced Obliquity Variations for the Major Satellites of Saturn[END_REF] give an obliquity of roughly 0.001 • but we could not find the source of the value. However there are three decades separating Voyager's flybys and the mission of Cassini, therefore we believe that a new study of the orientation of Iapetus and a comparison to a dynamical model can bring more insights into the history of the satellite's motion.

Also, a fast orbital migration of Titan brings the 4:3 mean motion resonance with Hyperion into reconsideration. One should show that this resonance is compatible with a quality factor 100 < Q < 2000 at the frequency of Titan (like in [START_REF] Ćuk | Titan-Hyperion Resonance and the Tidal Q of Saturn[END_REF]. In addition, although we have assumed in this thesis that Titan had the same eccentricity as today, one can also speculate on its recent excitation. The resonance with Hyperion gives an upper limit to its eccentricity (around 0.06 or 0.07 [START_REF] Colombo | On the formation of the orbit-orbit resonance of Titan and Hyperion[END_REF] at the time of the capture. A few simulations of the 5:1 resonance with Iapetus show that a low eccentricity of Titan would avoid Iapetus' ejection whereas an eccentricity of 0.06 would perturb Iapetus even more. Therefore we are optimistic that Titan's past eccentricity can be constrained by the resonances of Titan, Hyperion and Iapetus.

Using Mathematica, a large number of their digits have been computed. The following is the decimal value of the nodes with 50 digits.

A.3.2 Roots specification

It is known that the nodes are the ones used in Radau quadrature, a Gaussian-like quadrature to evaluate numerically the integration of a function [START_REF] Hildebrand | Introduction to Numerical Analysis: 2Nd Edition[END_REF]. For a given regular function f , the evaluation of its integral can be made using

1 -1 f (t)dt = 2 n 2 f (-1) + n-1 i=1 W i f (x i ) + E (A.35)
where n is the number of abscissas where the function is evaluated, including -1 (here m = 8) and W i are weights which can be found in [START_REF] Hildebrand | Introduction to Numerical Analysis: 2Nd Edition[END_REF] x i are here the free nodes in the interval [-1, 1] and are solutions of .36) where P 7 and P 8 are Legendre polynomials of degree 7 and 8. The nodes computed in Equation A.34 can then be related to those solutions by x i = 2h i -1.

P 7 (x) + P 8 (x) 1 + x = 0 (A
During the implementation of the code, we have wondered if there existed any kind of shortcut to make the program more efficient. In our research, we wanted to shorten the evaluation of the functions involving the nodes by simplifying them. We wanted to find closed-form solutions for the roots but we only ended up proving that they are irrational and we managed to give a further relation involving all the roots. When multiplying Equation A.33 by 6435, we end up with a polynomial with integer coefficients 6435h 7 -24024h 6 + 36036h 5 -27720h 4 + 11550h3 -2520h 2 + 252h -8 (A.37)

Then the Rational root theorem implies that if a root of the polynomial was a rational number p q , then its numerator would be a factor of the constant term -8 p = {1, 2, 4, 8} A.38) and its denominator is a factor of the leading coefficient 6435. , 3, 5, 9, 11, 13, 15, 33, 39, 45, 55, 65, 99, 117, 143, 165, 195, 429, 495, 585, 715, 1287, 2145, 6435} (A.39) Knowing that the root is positive and less than 1 makes 90 possible rational candidates, but after plugging them in the polynomial, we found out that none of them was a root. Therefore these are irrational numbers. This property can be illustrated by the fact that there are no redundant patterns in the digits of the roots. in which t 0 is the epoch at which the body is at its pericentre. Its time derivative dM dt = n is called the mean motion and is related to the semi-major axis through the third law of Kepler

q = {1
n 2 a 3 = G(M + m). (B.3)
The true position of the body on its orbit is given by the true anomaly, which is not part of the set of orbital elements. To compute it from the mean anomaly, we first calculate the eccentric anomaly via Kepler's equation

B.2 Lagrange Planetary Equations

In the presence of a perturbing force, the trajectory of the body is expected to be of a different shape. However the conic still stands as a basic trajectory to depict perturbed orbits. During a perturbed evolution of a body, its position and velocity can be, at any time, related to a set of corresponding osculating orbital elements. At each moment, its instantaneous trajectory would be on a tangent to a conic. Lagrange's Planetary Equations govern the time rates of the osculating elements [START_REF] Murray | Solar System Dynamics[END_REF].

                                       da dt = 2 na ∂R ∂M de dt = Ecc 2 na 2 e ∂R ∂M -Ecc na 2 e ∂R ∂ω di dt = 1 na 2 Ecc sin(i) cos(i) ∂R ∂ω -∂R ∂Ω dΩ dt = 1 na 2 Ecc sin(i) ∂R ∂i dω dt = Ecc na 2 e ∂R ∂e - cos(i) na 2 Ecc sin(i) ∂R ∂i dM dt = n -Ecc 2 na 2 e ∂R ∂e -2 na ∂R ∂a (B.10)
where Ecc = √ 1e 2 . These classical elements defined above are not well suited to study the perturbed motion of a satellite. As shown in the next appendix, instead of choosing the angles (Ω, ω, M ), we prefer the angular set (Ω, ̟, λ) which appears in the expansion of the disturbing function. The new set of elements is defined the following way

               a = a e = e i = i Ω = Ω ̟ = Ω + ω λ = Ω + ω + M . ( B
= ∂R ∂λ de dt = -E(1-E) na 2 e ∂R ∂λ -E na 2 e ∂R ∂̟ di dt = - 1 na 2 E sin(i) ∂R ∂Ω - tan( i 2 ) na 2 E ∂R ∂λ + ∂R ∂̟ dΩ dt = 1 na 2 E sin(i) ∂R ∂i d̟ dt = E na 2 e ∂R ∂e + tan( i 2 ) na 2 E ∂R ∂i dλ dt = n -2 na ∂R ∂a + E(1-E) na 2 e ∂R ∂e + tan(i/2) na 2 E ∂R ∂i
.

(B.16)

One can see that those equations are not defined for circular orbits and for orbits with null inclination. This set of orbital elements show some singularities. Whenever one studies the evolution of the orbit of a body on a trajectory close to a circular orbit (e close to 0) or with a orbital plane close to the reference plane (i close to 0), those differential equations can bring wrong solutions when solved numerically or analytically. Therefore we use the following set of non-singular orbital elements

               a = a k = e cos (̟) h = e sin (̟) q = sin i 2 cos (Ω) p = sin i 2 sin (Ω) λ = λ (B.17)
and we derive here their corresponding Lagrange Planetary Equations. As before, we express the time derivative of the function R, now being a function of the set (a, k, h, q, p, λ)

∂R ∂i = 1 2 cos i 2 cos(Ω) ∂R ∂q + sin(Ω) ∂R ∂p (B.21b) ∂R ∂Ω = q ∂R ∂p -p ∂R ∂q (B.21c) ∂R ∂̟ = k ∂R ∂h -h ∂R ∂k (B.21d)
Finally, after all calculations have been made, the Lagrange Planetary Equations for non-singular elements are

                                             da dt = 2 na ∂R ∂λ dk dt = -1 na 2 k E E+1 ∂R ∂λ + E ∂R ∂h + h 2E q ∂R ∂q + p ∂R ∂p dh dt = -1 na 2 h E E+1 ∂R ∂λ -E ∂R ∂k -k 2E q ∂R ∂q + p ∂R ∂p dq dt = -1 4na 2 E ∂R ∂p -2q h ∂R ∂k -k ∂R ∂h -∂R ∂λ dp dt = 1 4na 2 E ∂R ∂q + 2p h ∂R ∂k -k ∂R ∂h -∂R ∂λ dλ dt = n + 1 na 2 -2a ∂R ∂a + E E+1 k ∂R ∂k + h ∂R ∂h + 1 2E q ∂R ∂q + p ∂R ∂p (B.22)
Application As shown in Appendix C, the disturbing function expressed in the classical orbital elements can be expanded as a sum of cosines. Now in order to make use of the non-singular elements, those cosines term have to be expressed in terms of the non-singular elements. First for the amplitudes, which are function of the eccentricities and inclinations, we simply invert the relations in Equation B.17

e = √ k 2 + h 2 s = q 2 + p 2 (B.23)
where s = sin i 2 . Then for the cosine itself, we have the general argument φ = j 1 λ ′ + j 2 λ + j 3 ̟ ′ + j 4 ̟ + j 5 Ω ′ + j 6 Ω (B.24) and using imaginary notation, we have

cos (φ) = ℜ exp I j 1 λ ′ + j 2 λ + j 3 ̟ ′ + j 4 ̟ + j 5 Ω ′ + j 6 Ω (B.25)
with I 2 = -1. Then by using the multiplicative property of the exponential function

cos (φ) = ℜ   exp I(j 1 λ ′ + j 2 λ) ×   k ′ + |j 3 | j 3 Ih ′ e ′   |j 3 | ×   k + |j 4 | j 4 Ih e   |j 4 | ×   q ′ + |j 5 | j 5 Ip ′ s ′   |j 5 | ×   q + |j 6 | j 6 Ip s   |j 6 |    (B.26)
and then by multiplying all factors together and separating the real part from the imaginary part, we finally get the general expression

e ′|j 3 | e |j 4 | s ′|j 5 | s |j 6 | cos(φ) = cos(j 1 λ ′ + j 2 λ)A(k, k ′ , h, h ′ , q, q ′ , p, p ′ )+ sin(j 1 λ ′ + j 2 λ)B(k, k ′ , h, h ′ , q, q ′ , p, p ′ ) (B.27)
where A and B are polynomial functions of the non-singular variables. For example, for a resonance term, we have

e ′ es ′ s cos(5λ ′ -λ -̟ ′ -̟ -Ω ′ -Ω) = (hh ′ pp ′ -kk ′ pp ′ -h ′ kp ′ q -hk ′ p ′ q -h ′ kpq ′ -hk ′ pq ′ -hh ′ qq ′ + kk ′ qq ′ ) cos(5λ ′ -λ)
+ (h ′ kpp ′ + hk ′ pp ′ + hh ′ p ′ qkk ′ p ′ q + hh ′ pq ′kk ′ pq ′h ′ kqq ′hk ′ qq ′ ) sin(5λ ′λ) (B.28) and for a secular term e ′2 e 2 cos 2̟ ′ -2̟ = h 2 h ′2h ′2 k 2 + 4hh ′ kk ′h 2 k ′2 + k 2 k ′2 .

(B.29)

For an external perturber, the external indirect part is

R E = - r i a i a j r j 2 cos(Ψ) (C.6)
where Ψ is the angular distance between the two satellites cos(Ψ) = r j .r i r j r i .

(C.7)

For a body orbiting on a conic, the position vector of a body expressed in terms its components is So that, after writing Equation C.7 as cos(Ψ) =

x j x i +y j y i +z j z i r j r i

, we obtain cos(Ψ) = 1s 2 is 2 j + s 2 i s 2 j cos (θ iθ j ) + s 2 i 1s 2 j cos (θ i + θ j -2Ω i ) + s 2 j 1s 2 i cos (θ i + θ j -2Ω j ) -2s i 1s 2 i s j 1s 2 j cos (θ i + θ j -Ω i -Ω j ) + 2s i 1s 2 i s j 1s 2 j cos (θ iθ j -Ω i + Ω j ) + s 2 i s 2 j cos (θ iθ j -2Ω i + 2Ω j ) (C.9)

where θ = Ω + ω + f and s = sin I 2 . Equation C.9 can also be found in [START_REF] Delaunay | Théorie du mouvement de la lune[END_REF]

C.2 Elliptical expansion

The conic's equation is (C.14) There cosine terms in Equation C.9 can be expanded using cos(θ iθ j ) = cos (5λ iλ j -4ω i ) e 4 i -9 128 cos (3λ i + λ j -4ω i ) e 4 i + cos (4λ iλ j -3ω i ) e 3 i -1 12 cos (2λ i + λ j -3ω i ) e 3 i e j cos (4λ i -3ω iω j ) e 3 i -1 12 e j cos (2λ i + 2λ j -3ω iω j ) e 3 i + 1 e j cos (2λ i -3ω i + ω j ) e 3 i + 4 3 e j cos (4λ i -2λ j -3ω i + ω j ) e 3 i -9 e 2 j cos (3λ i + λ j -2ω i -2ω j ) e 2 i -9 64 e 2 j cos (λ i + 3λ j -2ω i -2ω j ) e 2 i e j cos (3λ i -2ω iω j ) e 2 i -1 8 e j cos (λ i + 2λ j -2ω iω j ) e 2 i + e j cos (λ i -2ω i + ω j ) e 2 i + 9 8 e j cos (3λ i -2λ j -2ω i + ω j ) e 2 i + e 2 j cos (3λ i -3λ j -2ω i + 2ω j ) e 2 i + 1 64 e 2 j cos (λ iλ j -2ω i + 2ω j ) e 2 i -1 e 3 j cos (2λ i + 2λ jω i -3ω j ) e i -4 3 e 3 j cos (4λ jω i -3ω j ) e i + 1 e 3 j cos (2λ j + ω i -3ω j ) e i -1 8 e 2 j cos (2λ i + λ jω i -2ω j ) e i e 2 j cos (3λ jω i -2ω j ) e i + 1 8 e 2 j cos (λ j + ω i -2ω j ) e i + e j cos (ω iω j ) e i + 9 8 e 2 j cos (2λ i -3λ jω i + 2ω j ) e i + e 3 j cos (2λ i -4λ jω i + 3ω j ) e i + 4e 2 j e i + e i cos (2λ iλ jω i ) + e i e 2 je i cos (λ jω i ) -9 e 4 j cos (λ i + 3λ j -4ω j ) -1 12 e 3 j cos (λ i + 2λ j -3ω j ) + cos (λ i -3λ j + 2ω j ) + 4 3 e 3 j cos (λ i -4λ j + 3ω j ) + e 4 j cos (λ i -5λ j + 4ω j ) .

(C.15)

C.3 Legendre Expansion

First, let's derive an expansion for the gravitational perturbation of the Sun acting on one of the satellite. Here ⊙ will be used as a subscript for the elements of the Sun and s for the satellite. We assume that the Sun's motion around Saturn is unperturbed, so for the satellite, which is the internal body here, the equation of motion reads

C.4 Expansion using explicit expressions

In this thesis, we made an extended use of semi-analytic models using explicit developments appearing in [START_REF] Murray | Solar System Dynamics[END_REF] An amplitude term for a cosine in the direct part can be computed using × cos [j 1 λ 8 + j 2 λ 6 + j 3 ̟ 8 + j 4 ̟ 6 + j 5 Ω 8 + j 6 Ω 6 ] (C.23) which computes the amplitude of the argument φ = j 1 λ 8 + j 2 λ 6 + j 3 ̟ 8 + j 4 ̟ 6 + j 5 Ω 8 +j 6 Ω 6 . We note that one should use this relation twice in order to get the correct amplitude term. One with the argument φ and a second one with -φ because of the parity of the cosine function.

R D = imax i=0 2i! i! (-1) i 2 2i+1 α i
R φ = S + cos(φ) + S -cos(-φ).

(C.24)

The indirect parts are

R E = -κ m (1 -m)!
(1 + m)! F 1,m,p (I 6 )F 1,m,p ′ (I 8 )X 1,-j 2 -j 4 -j 2 (e 6 )X -2,j 1 +j 3 j 1 (e 8 )

× cos [j 1 λ 8 + j 2 λ 6 + j 3 ̟ 8 + j 4 ̟ 6 + j 5 Ω 8 + j 6 Ω 6 ] (C.25)

R I = -κ m (1 -m)!
(1 + m)! F 1,m,p (I 6 )F 1,m,p ′ (I 8 )X -2,-j 2 -j 4 -j 2 (e 6 )X 1,j 1 +j 3 j 1 (e 8 )

× cos [j 1 λ 8 + j 2 λ 6 + j 3 ̟ 8 + j 4 ̟ 6 + j 5 Ω 8 + j 6 Ω 6 ] . • α, a real number which is the ratio of the two semi-major axes.

D.2 Computation

One would compute a Laplace coefficient by using a numerical quadrature to evaluate the above integral. But it can actually be expanded to a series in α, the latter being comprised between 0 and 1. By Taylor expanding the integrand and then integrating with respect to Ψ, on can show that a Laplace coefficient can be expanded in a hypergeometric series [START_REF] Brouwer | Methods of celestial mechanics[END_REF]) [START_REF] Murray | Solar System Dynamics[END_REF] b j s (α) = 2ζ j s α j F (s, j, α 2 ) (D.4)

F being a series of the form F (s, j, α 2 ) = Hence, according to D'Alembert's ratio test, this series will converge to the value of the Laplace coefficient.

D.3 Relation between coefficients

Many useful relations exist between different Laplace coefficients and their derivatives. Here we are showing a few of them with their demonstration. The expression corresponds to a trigonometric series of the null function. Therefore we conclude that its Fourier coefficients are all equal to zero. For any j we then have The first term equals zero and the numerator in the integrand of the second term can be rewritten sin [(j + 1)Ψ] sin(Ψ) = 1 2 (cos(jΨ)cos [(j + 2)Ψ]) so that, after changing j + 1 to j, we have the relation For the last manoeuvre, we want to get rid of b j-1 s+1 . In order to do so, from Equation D.20, we express b j-1 s+1 as a function of the other two coefficients, which we plug into Equation D.22. 384 -1728α 2 + 2958α 4 -2259α 6 + 396α 8 -111α 10 -24α 12 α 6 (1α 2 ) 5 .

(D.30)

FIGURE 3 :

 3 FIGURE 3: From JPL Horizon

  FIGURE 1.3: Scheme depicting the tidal response of the primary. Here we are in the case of a satellite orbiting above the synchronous orbit. The tidal bulge, driven by the faster rotation of the primary is ahead of the direction of the satellite. The angle β represents the geometrical lag between the satellite's direction and the axis of the bulge.

FIGURE 1

 1 FIGURE 1.5: The gravitational pull of the satellite on a random point of the planet.

FIGURE 1 . 7 :

 17 FIGURE 1.7: Detailed scheme of Saturn's interior. The dissipation can either come from the non-elastic solid core or viscous friction in fluid element of the planet.

FIGURE 1 . 8 :

 18 FIGURE 1.8: Evolution of Titan's semi-major axis with the model of a constant tidal lag Equation1.17. Here Q was set to 10 -2 .

FIGURE 1

 1 FIGURE 1.9: Semi-major axis evolution of the satellites of Saturn under Fuller's model. No mutual gravitational influence between satellites was introduced.

FIGURE 1 .

 1 FIGURE 1.10: Quality factors of Saturn at the frequency of icy satellites according to the model in Fuller, Luan, and Quataert, 2016. This figure comes from the same article. LJ15 refers to the quality factor computed in Lainey et al., 2017

FIGURE 2 . 2 :

 22 FIGURE 2.2: Tilt to the Laplace plane for the orbital plane of Iapetus.It is computed using cos(J) = cos(i) cos(i LP ) + sin(i) sin(i LP ) cos(Ω -Ω ⊙ ) where the LP stands for the elements of the Laplace plane (chapter 4), Ω ⊙ is the ascending node of the Sun and J is the tilt.

FIGURE 2 . 3 :

 23 FIGURE 2.3: Motion of Titan's orbital plane over 10 000 years. Both the inclination and the longitude of the ascending node oscillate around the fixed values of the local Laplace plane with a period of around 705 years. The figure comes from a simulation done with all the major satellites of Saturn (Mimas, Enceladus, Tethys, Dione, Rhea, Titan, Hyperion, Iapetus), the Sun and the effect of the flattening of Saturn (Table 2.1). The angular elements of the Laplace plane are approximately 0.635 degrees for the inclination and 224.5 degrees for the ascending node. The initial conditions are those of the J2000 epoch taken from JPL ephemeris, Horizon (SAT 389.14).

FIGURE 3 . 1 :

 31 FIGURE 3.1: Evolution of Iapetus inclination under only the influence of the Sun. Simulations include a N-Body integration (in blue), an integration with the semi-analytic code comprising terms in Table 3.2 (in gold) and a simulation integration using in addition terms in Table3.3 (in green). These evolutions do not represent any real dynamics but are test-simulations to assess the reliability of the solar expansion.

FIGURE 3 . 2 :

 32 FIGURE 3.2: Evolution of Titan's elements using the N-Body code, the first expansion and the extended one. The two semi-analytical approach output the same evolution and approximate fairly well the evolution of Titan's eccentricity and inclination.

FIGURE 3 . 3 :

 33 FIGURE 3.3: Evolution of Iapetus' inclination using the N-Body code, the first expansion and the extended one. The two semi-analytical approach output the same evolution and approximate well the N-Body for the inclination of Titan.

FIGURE 3 . 4 :

 34 FIGURE 3.4: Iapetus eccentricity as a function of time.

  FIGURE 3.5: Iapetus eccentricity as a function of time.

FIGURE 3 . 6 :

 36 FIGURE 3.6: Iapetus argument of perichrone as a function of time.

FIGURE 5

 5 FIGURE 5.1: N-Body simulations with Q=100 for two different initial mean anomalies of Iapetus. On the left M 8,0 = 0 was chosen and on the right M 8,0 = π.

FIGURE 5 . 2 :

 52 FIGURE 5.2: Positions of resonant terms for the extended semianalytic model. Many are overlapping. The exact resonance is α r = 1 5 2/3 = 0.341 995 189 335 340

FIGURE 5 . 3 :

 53 FIGURE 5.3: Example figures of outputs obtained from our simulations. In this example we clearly distinguish a regular evolution before the resonance, a chaotic resonance crossing and finally regular dynamics after the resonance. In this specific example, the eccentricity grows up to 0.12 before decreasing and finally settling to around 0.05 and the tilt reached 8 degrees before settling at 4.5 degrees after the two bodies for released from resonance.

FIGURE 5 . 4 :

 54 FIGURE 5.4: Smooth evolutions of the eccentricity and the inclination coming from two different simulations. On the left, the eccentricity of Iapetus growths smoothly after both satellites enter the resonance, while the tilt stays constant. The right figure depicts the evolution of the inclination with respect to the Laplace plane. For this simulation the eccentricity, does not change.

FIGURE 5

 5 FIGURE 5.5: Evolution of the angle 5λ 8λ 6 -4̟ 8 during a smooth evolution of the eccentricity. The left figure is its computed value as a function of time and the right figure depicts its time rate. The speed shows rapid variation on short time scales but its average value is around -12 degrees per year.

FIGURE 5

 5 FIGURE 5.6: Post-resonant eccentricities and inclinations using Q=20.

FIGURE 5

 5 FIGURE 5.7: Post-resonant eccentricities and inclinations using Q=50

FIGURE 5

 5 FIGURE 5.8: Post-resonant eccentricities and inclinations using Q=75

FIGURE 5 .

 5 FIGURE 5.10: Example of a simulation showing a relatively long chaotic evolution of the eccentricity and the inclination of Iapetus during the resonance. The crossing is obvious when the eccentricity and the inclination vary rapidly and chaotically. As for all chaotic evolution, the post-resonance values reached are unpredictable. This particular simulation, done with Q = 75, shows a 0.06 eccentricity growth and a 1 degree growth in the tilt with respect to the local Laplace plane of Iapetus.

FIGURE 5 .

 5 FIGURE 5.13: Post-resonant eccentricities and inclinations using Q=400. Out of 100 simulations, 9 end with the ejection of Iapetus and 23 show a smooth evolution of its eccentricity.

FIGURE 5 .

 5 FIGURE 5.14: Post-resonant eccentricities and inclinations using Q=600. Out of 100 simulations, 19 end with the ejection of Iapetus, 2 show a smooth evolution of its inclination and 24 show a smooth evolution of its eccentricity.

FIGURE 5 .

 5 FIGURE 5.16: Simulation showing Iapetus being on an ejection trajectory. In the first place, chaos makes the eccentricity vary rapidly on a relatively short time scale, then reaches 0.2 in a monotonic way. Once it has reached that value 20 millions years after the start of the simulation, the eccentricity seems to settle for a while before being excited again and growing drastically towards ejection (e>1). The settling of the eccentricity at around 0.2 illustrates the action of sub-resonances during the process.

  Figure 5.11, Figure 5.12, Figure 5.13, Figure 5.14 and Figure 5.15 show the post-resonant elements obtained in the case where Iapetus was released from the resonance for respectively Q = 100, Q = 200, Q = 600 and Q = 1500.

FIGURE 5 .

 5 FIGURE 5.19: Example of a simulation done with Q=600. Here the eccentricity stays relatively low while the inclination grows to reach 4.5 degrees after the release.

FIGURE

  FIGURE 5.21: Q=20

FIGURE 5 .

 5 FIGURE 5.22: Simulation with the highest tilt done with Q = 20.Many of simulations done with Q = 20 shows just a "kick", just as Figure5.9, but others show a capture and a chaotic evolution of the eccentricity and the inclination.

FIGURE 5 .

 5 FIGURE 5.23: Post-resonant eccentricities and inclinations using Q=200. Out of 100 simulations, 22 end with the ejection of Iapetus.

FIGURE 5 .

 5 FIGURE 5.24: Post-resonant eccentricities and inclinations using Q=400. Out of 100 simulations, 22 end with the ejection of Iapetus.

FIGURE 5 .

 5 FIGURE 5.25: Post-resonant eccentricities and inclinations using Q=600. Out of 100 simulations, 36 end with the ejection of Iapetus.

FIGURE 5 .

 5 FIGURE 5.26: Post-resonant eccentricities and inclinations using Q=800. Out of 100 simulations, 46 end with the ejection of Iapetus.

FIGURE 5 .

 5 FIGURE 5.27: Simulation with a high inclination and a high eccentricity using Q=800.

FIGURE 5 .

 5 FIGURE 5.34: Post-resonant eccentricities and inclinations using Q=18000. Out of 100 simulations, 84 end with the ejection of Iapetus.

  FIGURE 5.35: Number of ejections as a function of the effective quality factor of the planet. Each dot is representative of a run of 100 simulations performed with the same value of Q. We observe that the number of ejections using the semi-analytic model is greater than with the N-body. The reason is probably that the rates of mean longitudes are not properly modelled when we average out all the terms except the resonant ones and such biases appear during the resonance capture where mean longitudes play an important role.

FIGURE 6

 6 FIGURE 6.1: Post-resonant eccentricities and inclinations using Q = 200 and an initial tilt of 5 degrees. Out of 100 simulations, 33 end with the ejection of Iapetus.

FIGURE 6

 6 FIGURE 6.2: Post-resonant eccentricities and inclinations using Q = 100 using quadruple floating precision numbers. Out of 100 simulations, 2 end with the ejection of Iapetus.

FIGURE 6 . 3 :

 63 FIGURE 6.3: Inclination of Iapetus with respect to two different planes : on the left the reference plane is the initial frame and on the right the reference plane is the plane of the precessing equator. The precession period of Saturn's pole was set to 150 thousand years. It is slow enough for the orbital planes of both satellites to follow the precession of the planet.

  (Ω) cos(ω + f )sin(Ω) sin(ω + f ) cos(I) y = sin(Ω) cos(ω + f ) + cos(Ω) sin(ω + f ) cos(I) z = sin(ω + f ) sin(I)

(

  f ) is a periodic function in M , therefore it can be expanded in its Fourier series e cos(M )+e 2 1 2 (1cos(2M ))+e 3 3 8 (cos(M )cos(3M ))+e 4 1 3 (cos(2M )cos(4M ))+O e 5 (4M )cos(2M ))+O e 5

4

  i + λ j -2ω j ) + e 2 i e je j cos (λ iω j ) + 5e 3 i e j4e i e j cos (2λ iω iω j ) +5e i e 3 j 4e i e j cos (2λ jω iω j ) 2 i e j + e j cos (λ i -2λ j + ω j ) + e j e i cos (2λ i -2λ jω i + ω j )

  ∇ rs R ⊙ (C.16)where the disturbing function R ⊙ is the one of an external satellite.R ⊙ = Gm ⊙ a ⊙ (R D,⊙ + αR E,⊙) rs r ⊙ is small, we Taylor expand Equation C.18. First we have the following expansion 1 P l are the Legendre polynomials which can be computed using the following relationP l (x)

  4n + 1)(sn)! 2 2n n!(2s -2n + 1)! s-2n m=0 κ m (s -2nm)! (s -2n + m)! × (-1) s-2n-m F s-2n,m,p (I 6 )F s-2n,m,p ′ (I 8 )

F

  and X are the inclination and eccenctricity functions which can be found in[START_REF] Ellis | The Disturbing Function in Solar System Dynamics[END_REF], we have for the function of inclinationF s-2n,m,p (I) = 1 2 s-2n (s -2n)! × min(p,[(s-2n-m)/2]) t=0 (2s -4n -2t)! (s -2nm -2t)! s -2n t × sin s-2n-m-2t (I) α = max(0, cb) and β = max(0, bc) and X a,b c,d are rational numbers called Newcomb operators defined in a recursive mannerX c -5d + 4 + 4b + a 4d X a,b c-1,d-1 + cd + b 2d j≥2 (-1) j 3/2 j X a,b c-j,d-j .a Taylor expansion of the ratio r r ′ around zero like in[START_REF] Murray | Solar System Dynamics[END_REF], one can expand around a fixed value close to the averaged value of the semi-major axes ratio. Denoting the latter A, one gets cos(Ψ) + A 2 (C.37)and is 2π-periodic in Ψ and can therefore be expanded in a trigonometric series. here the Laplace coefficients b j 1/2 . Continuing the expansion, the derivatives in Equation C.36 can be rewrittend n d( r r ′ ) n B( r r ′ = A) = d n dA n B(A), and exchanging the derivative with the summation sign, one has d n dA n B(A) expansions shown in Appendix C, one can show that the lowest power of the eccentricities in a ′ Therefore if we limit the expansion in eccentricity to the fourth power 1 then all power over 5 included can be neglected. By doing so, we change the direct part expansion to a finite series where N max denotes the order of the expansion in eccentricity.The expansion of the disturbing function[START_REF] Murray | Solar System Dynamics[END_REF] made the use of the expression 1 (1 -2α cos(Ψ) + α 2 ) s .(D.1)As this function is a 2π-periodic in the variable Ψ, it can be expanded in its Fourier series 1 (1 -2α cos(Ψ) + α 2 ) coefficients b j s are called Laplace coefficients and are functions of the semi-major axis ratio α. They can be calculated through the integral 2α cos(Ψ) + α 2 ) s dΨ.(D.3)It has three main parameters• s, which is a semi integer• j, an integer

+∞ n=0 FF

 n=0 n (s, j, α 2 ) (D.5) with coefficients F n following the sequence , in the recursion relation Equation D.6, that the numerator and denominator of the fraction are both O(n 2 ).

  First, using the integral expression in Equation D.3, one can see that because the integrand is an even function, one has, for any j to their definition in Equation D.2 we have, by changing the index s by s + 1,1 (1 -2α cos(Ψ) + α 2 ) s 1 1 -2α cos(Ψ) + α 2 cos(jΨ) cos(Ψ) = 1 2 (cos([j + 1]Ψ) + cos([j -1]Ψ)), we can rewrite the last term in the right-hand side as2α +∞ j=-∞ b j s+1 (α) cos(jΨ) cos(Ψ) = α +∞ j=-∞ b j s+1 (α) (cos([j + 1]Ψ) + cos([j -1]Ψ)) = α +∞ j=-∞ b j-1 s+1 (α) cos(jΨ) + +∞ j=-∞ b j+1 s+1 (α) cos(jΨ) (D.12) where we have also reordered the summation over j. We can expressed Equation D.11 using only one summation +∞ j=-∞ b j s -1 + α 2 b j s+1 + α b j+1 s+1 + b j-1 s+1 cos(jΨ) = 0. (D.13)

  2α cos(Ψ) + α 2 ) s dΨ (D.15) and by integrating by parts, splitting the integrand between cos [(j + 1)Ψ] 2α cos(Ψ) + α 2 ) s+1 .(D.16)

  for any j = 0. If we combine Equation D.14 and Equation D.17, we can derive the following important and useful relation enable us to express any Laplace coefficient as a function of others with lower j indexes and keeping the s index unchanged. We can derive another important equality be first changing j to j + 1 and s to s + 1 in the previous relation simply change this last equality by using j + 1 instead of j, + s) 1 + α 2 (js)(js + 1) + b j s+1 2α 2 (js)j 1 + α 2 2 α(js)(js + 1) s.(D.22)

Finally

  , after some arrangement and shifting s + 1 to s, we obtain the second useful relation analogous to Equation D.18 for lowering the s index.Equation D.18 and Equation D.23 are two relations with which one can express any Laplace coefficient b j s as a function of two main coefficients b 0 1/2 and b 1 1/2 . For a given coefficient b j s , we first recursively use Equation D.23 until we have expressed it as a combination of Laplace coefficient having s = 1/2. Then after using Equation D.18 as much as needed, a Laplace coefficient can be expressed as the direct part of the perturbing function in Equation C.23 also involves derivatives of Laplace coefficients, we show here that they can also be expressed functions of the principal Laplace coefficients. Let be the general Laplace coefficient be b j s . Its derivative with respect to α can be written as 2α cos(Ψ) + α 2 ) s . (D.25)Then one can swop the order of the derivation and the integration. One hasd dα cos(jΨ) (1 -2α cos(Ψ) + α 2 ) s = 2s cos(jΨ) (cos(Ψ)α)(1 -2α cos(Ψ) + α 2 ) s+1 . for any n ≥ 2. All the boxed equations derived above can be used to express any derivative of any Laplace coefficient as a function of only two principal coefficients. For example d 3 dα 3 b 4 3/2 =b 0 1/2 -192 + 840α 2 -1383α 4 + 1050α 6 + 57α 8 + 12α 10 α 5 (1α 2 )

TABLE 1 :

 1 Osculating orbital elements for the epoch 2 451 545 (January 1, 2000 and 12 hours). The reference plane is the center of the planet and the mean equator and the values are taken from JPL Horizon website.

TABLE 1

 1

.1: Values of Q needed for each satellite to migrate from the synchronous orbit to their current position in the age of the Solar System (4.5 billion years)

TABLE 2

 2

	.1: List of gravitational parameters for the planet Saturn.
	All values are taken from the latest version of the JPL ephemerides
	system HORIZON SAT 389.14 ftp://ssd.jpl.nasa.gov/pub/
	eph/satellites/nio/LINUX_PC/sat389xl.txt. The mass of
	Saturn is expressed in solar masses.

TABLE 2

 2

	.2: Contribution of each inner satellite to the flattening
	of the system. Values of semi-major axes are taken from the JPL
	Horizon website https://ssd.jpl.nasa.gov/?sat_elem and
	masses are the ones given from the SAT 389.14.

TABLE 3

 3

	.1: Main periods of the system considered. These values were
	computed from the simulation used for Figure 2.1 and Figure 2.3

TABLE 3

 3

.3: Terms having a dependence on the Sun's mean longitude (1 × λ ⊙ ). Note that no approximation was made concerning both inclinations.

TABLE 3

 3

.4: Table

containing

the secular part of the Titan-Iapetus interaction.

TABLE 3

 3

.6: Semi-major axis functions for the secular term # 2.

TABLE 3

 3

.7: Semi-major axis functions for the secular term # 3.

TABLE 3

 3

.8: Semi-major axis functions for the secular term # 4.

TABLE 3

 3

.10: Semi-major axis functions for the secular term # 6.

TABLE 3

 3

	#	Φ

.11: Semi-major axis functions for the secular term # 7.

TABLE 3

 3

.12: Semi-major axis functions for the secular term # 8.

TABLE 3 .

 3 13: Semi-major axis functions for the secular term # 9.

TABLE 3

 3

.14: Semi-major axis functions for the secular term # 10.

# Φ

Semi-major axis function 11

TABLE 3

 3

.15: Semi-major axis functions for the secular term # 11.

# Φ

Semi-major axis function 12

TABLE 3

 3

.16: Semi-major axis functions for the secular term # 12.

TABLE 3

 3

.17: Semi-major axis functions for the secular term # 13.

TABLE 3 .

 3 18: Semi-major axis functions for the secular term # 14.

	#	Φ	Semi-major axis function
		Φ 15,0
	15

TABLE 3

 3

.19: Semi-major axis functions for the secular term # 15.

# Φ

Semi-major axis function 16

TABLE 3

 3

	.20: Semi-major axis functions for the secular term # 16.

TABLE 3

 3

.21: Semi-major axis functions for the secular term # 17.

TABLE 3 .

 3 22: Semi-major axis functions for the secular term # 18.

	#	Φ	Semi-major axis function
	19

TABLE 3 .

 3 23: Semi-major axis functions for the secular term # 19.

	#	Φ	Semi-major axis function
	20

TABLE 3 .

 3 24: Semi-major axis functions for the secular term # 20.

	#	Φ	Semi-major axis function
	21	Φ 21,0 Φ 21,1 Φ 21,2	3 8 α 3 D + 3 16 α 4 D 2 b 0 5/2 32 α 3 + 75 16 α 4 D + 75 64 α 5 D 2 b 1 7/2 -75 64 α 4 + 315 315 64 α 5 D + 105 128 α 6 D 2 4b 0 9/2 + 3b 2 9/2

TABLE 3 .

 3 25: Semi-major axis functions for the secular term # 21.

	#	Φ	Semi-major axis function
	22	Φ 22,0 Φ 22,1	2 α 2 + 9 4 α 3 D + 3 8 α 4 D 2 b 1 5/2 -3 16 α 3 + 15 75 4 α 4 D + 15 32 α 5 D 2 3b 0 7/2 + 2b 2 7/2

TABLE 3

 3

		.26: Semi-major axis functions for the secular term # 22.
	#	Φ	Semi-major axis function
	23	Φ 23,0 Φ 23,1 -225 32 α 3 + 45 4 α 2 + 15 8 α 3 D + 3 16 α 4 D 2 b 2 5/2 -15 16 α 4 D + 45 64 α 5 D 2 3b 1 7/2 + 2b 3 7/2

TABLE 3

 3

	#	Φ	Semi-major axis function
	24	Φ 24,0 Φ 24,1 -135 64 α 3 + 135 75 64 α 3 + 75 32 α 4 D + 75 128 α 5 D 2 b 1 7/2 + 3 16 α 3 D + 3 32 α 4 D 2 b 0 5/2 32 α 4 D + 135 128 α 5 D 2 b 1 7/2 + 315 128 α 4 + 315 128 α 5 D + 105 256 α 6 D 2 9b 0 9/2 + 5b 2 9/2 Φ 24,2 135 64 α 3 + 135 32 α 4 D + 135 128 α 5 D 2 b 1 7/2 + 315 128 α 4 + 315 128 α 5 D + 105 256 α 6 D 2 12b 0 9/2 + 7b 2 9/2

.27: Semi-major axis functions for the secular term # 23.

TABLE 3 .

 3 28: Semi-major axis functions for the secular term # 24.

	#	Φ

TABLE 3 .

 3 29: Semi-major axis functions for the secular term # 25.

	#	Φ

TABLE 3 .

 3 30: Semi-major axis functions for the secular term # 26.

	#	Φ	Semi-major axis function
	27	Φ 27,0 Φ 27,1 Φ 27,2	1 256 α 5 D 3 + 1 1024 α 6 D 4 b 2 5/2 512 α 5 D 2 + 5 256 α 6 D 3 + 5 2048 α 7 D 4 b 1 7/2 + b 3 7/2 -15 1024 α 5 D + 105 35 2048 α 6 D 2 + 35 2048 α 7 D 3 + 35 24576 α 8 D 4 3b 0 9/2 + 8b 2 9/2 + 3b 4 9/2

TABLE 3 .

 3 31: Semi-major axis functions for the secular term # 27.

	#	Φ	Semi-major axis function
	28	Φ 28,0 Φ 28,1	128 α 4 D 2 + 5 128 α 5 D 3 + 1 256 α 6 D 4 b 1 5/2 -9 128 α 4 D + 15 45 32 α 5 D 2 + 35 256 α 6 D 3 + 5 512 α 7 D 4 b 0 7/2 + b 2 7/2
		Φ 28,2 -105 512 α 4 + 385 512 α 5 D + 525 1024 α 6 D 2 + 105 1024 α 7 D 3 + 35 6144 α 8 D 4 11b 1 9/2 + 3b 3 9/2

TABLE 3

 3

	.32: Semi-major axis functions for the secular term # 28.

TABLE 3 .

 3 33: Semi-major axis functions for the secular term # 29.

	#	Φ	Semi-major axis function
	30	Φ 30,0 Φ 30,1	32 α 2 + 75 64 α 3 D + 9 16 α 4 D 2 + 11 128 α 5 D 3 + 1 256 α 6 D 4 b 1 5/2 -15 128 α 3 + 735 525 128 α 4 D + 525 256 α 5 D 2 + 65 256 α 6 D 3 + 5 512 α 7 D 4 b 0 7/2 + b 2 7/2
		Φ 30,2 -735 128 α 4 + 735 128 α 5 D + 105 64 α 6 D 2 + 175 1024 α 7 D 3 + 35 6144 α 8 D 4 11b 1 9/2 + 3b 3 9/2

TABLE 3 .

 3 34: Semi-major axis functions for the secular term # 30.

	#	Φ	Semi-major axis function
	31	Φ 31,0 Φ 31,1 Φ 31,2	105 128 α 2 + 105 128 α 3 D + 63 256 α 4 D 2 + 7 256 α 5 D 3 + 1 1024 α 6 D 4 b 2 5/2 128 α 3 + 105 32 α 4 D + 105 128 α 5 D 2 + 5 64 α 6 D 3 + 5 2048 α 7 D 4 b 1 7/2 + b 3 7/2 -525 512 α 4 + 735 2205 256 α 5 D + 315 512 α 6 D 2 + 105 2048 α 7 D 3 + 35 24576 α 8 D 4 3b 0 9/2 + 8b 2 9/2 + 3b 4 9/2

TABLE 3

 3

.35: Semi-major axis functions for the secular term # 31.

TABLE 3 .

 3 36: First table containing the resonant part of the Titan-Iapetus interaction.

	# Argument (φ i )	Amplitude (S i )
	φ -ω8 -ω6 -3Ω 8 + Ω 6 φ -2ω 8 -3Ω 8 + Ω 6 φ + ω8 -ω6 -3Ω 8 -Ω 6 φ -3Ω 8 -Ω 6	e 8 e 6 s 6 s 3 8 1 -s 2 8 Φ 21,0 + s 2 8 Φ 21,1 + s 4 8 Φ 21,2 e 2 8 s 6 s 3 8 1 -s 2 8 Φ 22,0 + s 2 8 Φ 22,1 + s 4 8 Φ 22,2 e 6 e 8 s 6 s 3 8 Φ 23,0 + s 2 8 Φ 23,1 + s 4 8 Φ 23,2 8 1 -s 2 s 6 s 3 8 [ Φ 24,0,0 + e 2 6 Φ 24,0,1 + e 2 8 Φ 24,0,2 + 8 1 -s 2 s 2 8 Φ 24,1,0 + e 2 6 Φ 24,1,1 + e 2 8 Φ 24,1,2 +
	s 4 8 Φ 24,2,0 + e 2 6 Φ 24,2,1 + e 2 8 Φ 24,2,2 ] 8 1 -s 2 e 6 e 8 s 6 s 3 8 Φ 25,0 + s 2 8 Φ 25,1 + s 4 8 Φ 25,2 e 2 6 s 2 6 s 4 8 Φ 26,0 + s 2 8 Φ 26,1 + s 4 8 Φ 26,2 φ -ω8 -ω6 -4Ω 8 + 2Ω 6 e 6 e 8 s 2 φ -ω8 + ω6 -3Ω 8 -Ω 6 φ -2ω 6 -4Ω 8 + 2Ω 6 6 s 4 8 Φ 27,0 + s 2 8 Φ 27,1 + s 4 8 Φ 27,2 φ -2ω 8 -4Ω 8 + 2Ω 6 e 2 8 s 2 6 s 4 8 Φ 28,0 + s 2 8 Φ 28,1 + s 4 8 Φ 28,2 + s 6 8 Φ 28,3 φ + 2ω 8 -2ω 6 -4Ω 8 e 2 6 e 2 8 s 4 8 Φ 29,0 + s 2 8 Φ 29,1 + s 4 8 Φ 29,2 e 6 e 8 s 4 8 [ Φ 30,0,0 + e 2 6 Φ 30,0,1 + e 2 8 Φ 30,0,2 + s 2 6 Φ 30,0,3 + φ + ω8 -ω6 -4Ω 8 s 2 8 Φ 30,1,0 + e 2 6 Φ 30,1,1 + e 2 8 Φ 30,1,2 + s 2 6 Φ 30,1,3 ] s 4 8 Φ 30,0,0 + e 2 6 Φ 30,0,1 + e 4 6 Φ 30,0,2 + e 2 8 Φ 30,0,3 + e 4 8 Φ 30,0,4 + φ -4Ω 8 e 2 6 e 2 8 Φ 30,0,5 + s 2 6 Φ 30,0,6 + s 2 6 e 2 6 Φ 30,0,7 + s 2 6 e 2 8 Φ 30,0,8 +
		s 2 8 Φ 30,1,0 + e 2 6 Φ 30,1,1 + e 4 6 Φ 30,1,2 + e 2 8 Φ 30,1,3 + e 4 8 Φ 30,1,4 + e 2 6
	e 2 8 Φ 30,1,5 + s 2 6 Φ 30,1,6 + s 2 6 e 2 6 Φ 30,1,7 + s 2 6 e 2 8 Φ 30,1,8 e 6 e 8 s 4 8 [ Φ 32,0,0 + e 2 6 Φ 32,0,1 + e 2 8 Φ 32,0,2 + s 2 6 Φ 32,0,3 + s 2 8 Φ 32,1,0 + e 2 6 Φ 32,1,1 + e 2 8 Φ 32,1,2 + s 2 6 Φ 32,1,3 ] e 2 6 e 2 8 s 4 8 Φ 33,0 + s 2 8 Φ 33,1 + s 4 8 Φ 33,2 e 2 8 s 2 6 s 4 8 Φ 34,0 + s 2 8 Φ 34,1 + s 4 8 Φ 34,2 + s 6 8 Φ 34,3 φ + ω8 + ω6 -4Ω 8 -2Ω 6 e 6 e 8 s 2 φ -ω8 + ω6 -4Ω 8 φ -2ω 8 + 2ω 6 -4Ω 8 φ + 2ω 8 -4Ω 8 -2Ω 6 6 s 4 8 Φ 35,0 + s 2 8 Φ 35,1 + s 4 8 Φ 35,2 + s 6 8 Φ 35,3 φ + 2ω 6 -4Ω 8 -2Ω 6 e 2 6 s 2 6 s 4 8 Φ 36,0 + s 2 8 Φ 36,1 + s 4 8 Φ 36,2 + s 6 8 Φ 36,3 φ + ω8 -ω6 -5Ω 8 + Ω 6 e 6 e 8 s 6 s 5 8 Φ 33,0 + s 2 8 Φ 33,1 8 1 -s 2 φ -5Ω 8 + Ω 6 s 6 s 5 8 [ Φ 38,0,0 + e 2 6 + e 2 8 Φ 38,0,1 + 8 1 -s 2 s 2 8 Φ 38,1,0 + e 2 6 + e 2 8 Φ 38,1,1 +
	φ -ω8 + ω6 -5Ω 8 + Ω 6 φ + 2ω 8 -5Ω 8 -Ω 6	s 4 8 Φ 38,2,0 + e 2 6 + e 2 8 Φ 38,2,1 ] 8 1 -s 2 e 6 e 8 s 6 s 5 8 Φ 39,0 + s 2 8 Φ 39,1 + s 4 8 Φ 39,2 e 2 8 s 6 s 5 8 1 -s 2 8 Φ 40,0 + s 2 8 Φ 40,1 + s 4 8 Φ 40,2

TABLE 3 .

 3 37: Second table containing the resonant part of the Titan-Iapetus interaction.

	# Argument (φ i )	Amplitude (S i )
	φ + ω8 + ω6 -5Ω 8 -Ω 6 φ + 2ω 6 -5Ω 8 -Ω 6 φ + ω8 -ω6 -6Ω 8 + 2Ω 6 e 6 e 8 s 2 e 6 e 8 s 6 s 5 8 1 -s 2 8 Φ 41,0 + s 2 8 Φ 41,1 + s 4 8 Φ 41,2 e 2 6 s 6 s 5 8 1 -s 2 8 Φ 42,0 + s 2 8 Φ 42,1 + s 4 8 Φ 42,2 6 s 6 8 Φ 43,0 + s 2 8 Φ 43,1 + s 4 8 Φ 43,2 s 2 6 s 6 8 [ Φ 44,0,0 + e 2 6 + e 2 8 Φ 44,0,1 + φ -6Ω 8 + 2Ω 6 s 2 8 Φ 44,1,0 + e 2 6 + e 2 8 Φ 44,1,1 +
	s 4 8 Φ 44,2,0 + e 2 6 + e 2 8 Φ 44,2,1 ] φ -ω8 + ω6 -6Ω 8 + 2Ω 6 e 6 e 8 s 2 6 s 6 8 Φ 45,0 + s 2 8 Φ 45,1 + s 4 8 Φ 45,2 + s 6 8 Φ 45,3 φ + 3ω 8 -ω6 -6Ω 8 e 3 8 e 6 s 6 8 Φ 46,0 + s 2 8 Φ 46,1 φ + 2ω 8 -6Ω 8 e 2 8 s 6 8 Φ 47,0 + s 2 8 Φ 47,1 φ + ω8 + ω6 -6Ω 8 e 6 e 8 s 6 8 Φ 48,0 + s 2 8 Φ 48,1 φ + 2ω 6 -6Ω 8 e 2 6 s 6 8 Φ 49,0 + s 2 8 Φ 49,1 φ -ω8 + 3ω 6 -6Ω 8 e 8 e 3 6 s 6 8 Φ 50,0 + s 2 8 Φ 50,1 φ + 2ω 8 -7Ω 8 + Ω 6 e 2 8 s 6 s 7 8 1 -s 2 8 Φ 51,0 + s 2 8 Φ 51,1 φ + ω8 + ω6 -7Ω 8 + Ω 6 e 6 e 8 e 6 s 7 8 1 -s 2 8 Φ 52,0 + s 2 8 Φ 52,1 φ + 2ω 6 -7Ω 8 + Ω 6 e 2 6 s 6 s 7 8 1 -s 2 8 Φ 53,0 + s 2 8 Φ 53,1 φ + 2ω 8 -8Ω 8 + 2Ω 6 e 2 8 s 2 6 s 8 8 Φ 54,0 + s 2 8 Φ 54,1 + s 4 8 Φ 54,2 φ + ω8 + ω6 -8Ω 8 + 2Ω 6 e 6 e 8 s 2 6 s 8 8 Φ 55,0 + s 2 8 Φ 55,1 + s 4 8 Φ 55,2 φ + 2ω 6 -8Ω 8 + 2Ω 6 e 2 6 s 2 6 s 8 8 Φ 56,0 + s 2 8 Φ 56,1 + s 4 8 Φ 56,2 φ + 4ω 8 -8Ω 8 e 4 8 s 8 8 Φ 57,0 + s 2 8 Φ 57,1 + s 4 8 Φ 57,2 φ + 3ω 8 + ω6 -8Ω 8 e 6 e 3 8 s 8 8 Φ 58,0 + s 2 8 Φ 58,1 φ + 2ω 8 + 2ω 6 -8Ω 8 e 2 6 e 2 8 s 8 8 Φ 59,0 + s 2 8 Φ 59,1 φ + ω8 + 3ω 6 -8Ω 8 e 3 6 e 8 s 8 8 Φ 60,0 + s 2 8 Φ 60,1 φ + 4ω 6 -8Ω 8 e 4 6 s 8 8 Φ 61,0 + s 2 8 Φ 61,1

TABLE 3 .

 3 38: Third table containing the resonant part of the Titan-Iapetus interaction.

	# Φ	Semi-major axis function
	1	Φ 1,0 Φ 1,1 -1961 1045 384 + 229 96 αD + 43 64 α 2 D 2 + 7 96 α 3 D 3 + 1 384 α 4 D 4 b 5 1/2 768 α + 179 96 α 2 D + 57 128 α 3 D 2 + 1 24 α 4 D 3 + 1 768 α 5 D 4 b 4 3/2 + b 6 3/2 Φ 1,2 3393 1024 α 2 + 529 256 α 3 D + 219 512 α 4 D 2 + 9 256 α 5 D 3 + 1 1024 α 6 D 4 b 3 5/2 + 4b 5 5/2 + b 7 5/2

TABLE 3 .

 3 39: Semi-major axis functions for the resonant term # 1.

	# Φ	Semi-major axis function
	2	4 + 533 48 αD + 97 32 α 2 D 2 + 5 16 α 3 D 3 + 1 96 α 4 D 4 b 4 1/2 Φ 2,0 -51 Φ 2,1 1145 96 α + 103 12 α 2 D + 127 64 α 3 D 2 + 17 96 α 4 D 3 + 1 192 α 5 D 4 b 3 3/2 + b 5 3/2
		Φ 2,2 -1969 128 α 2 + 1205 128 α 3 D + 483 256 α 4 D 2 + 19 128 α 5 D 3 + 1 256 α 6 D 4 b 2 5/2 + b 4 5/2 + b 6 5/2

TABLE 3 .

 3 40: Semi-major axis functions for the resonant term # 2.

	# Φ	Semi-major axis function
	3	Φ 3,0 Φ 3,1 -2635 1407 64 + 307 16 αD + 163 32 α 2 D 2 + 1 2 α 3 D 3 + 1 64 α 4 D 4 b 3 1/2 128 α + 235 16 α 2 D + 211 64 α 3 D 2 + 9 32 α 4 D 3 + 1 128 α 5 D 4 b 2 3/2 + b 4 3/2 Φ 3,2 13545 512 α 2 + 2043 128 α 3 D + 795 256 α 4 D 2 + 15 64 α 5 D 3 + 3 512 α 6 D 4 b 1 5/2 + b 3 5/2 + b 5 5/2

TABLE 3 .

 3 41: Semi-major axis functions for the resonant term # 3.

	# Φ	Semi-major axis function
	4	Φ 4,0 Φ 4,1	389 24 + 697 48 αD + 121 32 α 2 D 2 + 17 48 α 3 D 3 + 421 96 α 4 D 4 b 2 1/2 1475 96 α + 265 24 α 2 D + 155 64 α 3 D 2 + 19 96 α 4 D 3 + 1 192 α 5 D 4 b 1 3/2 + b 3 3/2
		Φ 4,2 -2535 128 α 2 + 1525 128 α 3 D + 579 256 α 4 D 2 + 21 128 α 5 D 3 + 1 256 α 6 D 4 b 0 5/2 + b 2 5/2 + b 4 5/2

TABLE 3 .

 3 42: Semi-major axis functions for the resonant term # 4.

	# Φ	Semi-major axis function
	5	Φ 5,0	523 128 + 389 96 αD + 67 64 α 2 D 2 + 3 32

TABLE 3 .

 3 43: Semi-major axis functions for the resonant term # 5.

	# Φ	Semi-major axis function
	6	Φ 6,0 Φ 6,1 -31 31 16 + 3 4 αD + 1 16 α 2 D 2 b 4 3/2 16 + 3 4 αD + 1 16 α 2 D 2 b 4 3/2 -129 32 α + 21 16 α 2 D + 3 32 α 3 D 2 b 3 5/2 + 5b 5 5/2 Φ 6,2 129 32 α + 21 16 α 2 D + 3 32 α 3 D 2 b 3 5/2 + 6b 5 5/2 +
			855 128	α 2 +	15 8	α 3 D +	15 128	α 4 D 2	b 2 7/2 + 16b 4 7/2 + 13b 6 7/2
		Φ 6,3 -855 128 α 2 + 15 8 α 3 D + 15 128 α 4 D 2 b 2 7/2 + 18b 4 7/2 + 15b 6 7/2

TABLE 3 .

 3 44: Semi-major axis functions for the resonant term # 6.

TABLE 3

 3

	# Φ	Semi-major axis function
	8	Φ 8,0 Φ 8,1 -85 16 α + 5 4 α 2 D + 1 16 α 3 D 2 b 2 3/2 85 16 α + 5 4 α 2 D + 1 16 α 3 D 2 b 2 3/2 +
			315 32	α 2 +	33 16	α 3 D +	3 32	α 4 D 2	b 1 5/2 + 5b 3 5/2
		Φ 8,2	315 32	α 2 +	33 16	α 3 D +	3 32	α 4 D 2	b 1 5/2 + 6b 3 5/2 +
			1905 128	α 3 +	3 32	α 4 D +	15 128	α 5 D 2	b 0 7/2 + 16b 2 7/2 + 13b 4 7/2

.45: Semi-major axis functions for the resonant term # 7.

TABLE 3 .

 3 46: Semi-major axis functions for the resonant term # 8.

	# Φ	Semi-major axis function
	9	8 α + 3 2 α 2 D + 1 8 α 3 D 2 b 4 3/2 Φ 9,0 -31 Φ 9,1 129 16 α 2 + 21 8 α 3 D + 3 16 α 4 D 2 b 3 5/2 + 2b 5 5/2
		Φ 9,2 -855 64 α 3 + 15 4 α 4 D + 15 64 α 5 D 2 b 2 7/2 + 6b 4 7/2 + 3b 6 7/2

TABLE 3 .

 3 47: Semi-major axis functions for the resonant term # 9.

	# Φ	Semi-major axis function
	10	Φ 10,0 Φ 10,1 -99 25 2 α + 4α 2 D + 1 4 α 3 D 2 b 3 3/2 4 α 2 + 27 4 α 3 D + 3 8 α 4 D 2 b 2 5/2 + 2b 4 5/2 Φ 10,2 315 8 α 3 + 75 8 α 4 D + 15 32 α 5 D 2 b 1 7/2 + 6b 3 7/2 + 3b 5 7/2

TABLE 3 .

 3 48: Semi-major axis functions for the resonant term # 10.

	# Φ	Semi-major axis function
	11	8 α + 5 2 α 2 D + 1 8 α 3 D 2 b 2 3/2 Φ 11,0 -85 Φ 11,1 315 16 α 2 + 33 8 α 3 D + 3 16 α 4 D 2 b 1 5/2 + 2b 3 5/2
		Φ 11,2 -1905 64 α 3 + 45 8 α 4 D + 15 64 α 5 D 2 b 0 7/2 + 6b 2 7/2 + 3b 4 7/2

TABLE 3 .

 3 49: Semi-major axis functions for the resonant term # 11.

	# Φ	Semi-major axis function
	12	Φ 12,0 Φ 12,1 -1611 1145 96 α + 233 48 α 2 D + 13 64 α 3 D 2 -7 96 α 4 D 3 -1 192 α 5 D 4 b 5 3/2 64 α 2 + 505 64 α 3 D -3 128 α 4 D 2 -9 128 α 6 D 4 b 4 5/2 + b 6 5/2 64 α 5 D 3 -1 Φ 12,2 2645 64 α 3 + 1255 128 α 4 D + 285 512 α 5 D 2 -55 256 α 6 D 3 -5 512 α 7 D 4 b 3 7/2 + 3b 5 7/2 + b 7 7/2

TABLE 3

 3

	# Φ	Semi-major axis function
	Φ 13,0,0	31 16
	13

.50: Semi-major axis functions for the resonant term # 12.

TABLE 3

 3

.51: Semi-major axis functions for the resonant term # 13.

TABLE 3

 3

	# Φ	Semi-major axis function
	Φ 15,0,0	85 16
	15

.52: Semi-major axis functions for the resonant term # 14.

TABLE 3 .

 3 53: Semi-major axis functions for the resonant term # 15.

	# Φ	Semi-major axis function
	16	96 α + 275 48 α 2 D + 115 96 α 3 D 2 + 17 96 α 4 D 3 + 1 192 α 5 D 4 b 1 3/2 Φ 16,0 -295 Φ 16,1 845 64 α 2 + 895 64 α 3 D + 447 128 α 4 D 2 + 19 64 α 5 D 3 + 1 128 α 6 D 4 b 0 5/2 + b 2 5/2
		Φ 16,2 -2175 64 α 3 + 3355 128 α 4 D + 2805 512 α 5 D 2 + 105 256 α 6 D 3 + 5 512 α 7 D 4 4b 1 7/2 + b 3 7/2

TABLE 3 .

 3 54: Semi-major axis functions for the resonant term # 16.

	# Φ	Semi-major axis function
	17	512 α 2 + 275 48 α 2 D + 115 96 α 3 D 2 + 17 96 α 4 D 3 + 1 192 α 5 D 4 b 1 3/2 Φ 17,0 -5 Φ 17,1 243 8 α 2 + 9 8 α 3 D -9 16 α 4 D 2 b 4 5/2 +
			105 4		α 3 -	15 32	α 5 D 2	3b 3 7/2 + 7b 5 7/2
		Φ 17,2 -	315 8	α 3 +	45 64	α 5 D 2	2b 3 7/2 + 5b 5 7/2 +
			735 16	α 4 -	105 64	α 5 D -	105 128	α 6 D 2	3b 2 9/2 + 19b 4 9/2 + 13b 6 9/2
		Φ 17,3	735 16 α 4 -105 64 α 5 D -105 128 α 6 D 2 3b 2 9/2 + 20b 4 9/2 + 14b 6 9/2

TABLE 3

 3

	# Φ	Semi-major axis function
	18	Φ 18,0,0 Φ 18,0,1 Φ 18,0,2 -423 9 4 α 2 b 3 5/2 9 8 α 2 + 27 8 α 3 D + 9 16 α 4 D 2 b 3 5/2 8 α 2 -27 16 α 4 D 2 b 3 5/2 8 α 3 D -9 4 α 2 b 3 5/2 + 15 8 α 3 3b 2 7/2 + 7b 4 7/2 Φ 18,1,0 -9 Φ 18,1,1 -9 8 α 2 + 27 8 α 3 D + 9 16 α 4 D 2 b 3 5/2 +
			15 4	α 3 +	15 4	α 4 D +	15 32	α 5 D 2 3b 2 7/2 + 7b 4 7/2
		Φ 18,1,2	423 8	α 2 -	27 8	α 3 D -	9 16	α 4 D 2 b 3 5/2 +
			165 4	α 3 -	15 4	α 4 D -	15 32	α 5 D 2	3b 2 7/2 + 7b 4 7/2
		Φ 18,2,0	45 16 α 3 2b 2 7/2 + 5b 4 7/2 + 105 32 α 4 3b 1 9/2 + 19b 3 9/2 + 13b 5 9/2
		Φ 18,2,1	45 8	α 3 +	45 8	α 4 D +	45 64	α 5 D 2	2b 2 7/2 + 5b 4 7/2 +
			105 8	α 4 +	525 64	α 5 D +	105 128	α 6 D 2	3b 1 9/2 + 19b 3 9/2 + 13b 5 9/2
		Φ 18,2,2 -			495 8	α 3 -	45 8	α 4 D -	45 64	α 5 D 2	2b 2 7/2 + 5b 4 7/2 +
			525 8	α 4 -	525 64	α 5 D -	105 128	α 6 D 2	3b 1 9/2 + 19b 3 9/2 + 13b 5 9/2
		32 α 4 3b 1 9/2 + 20b 3 9/2 + 14b 5 9/2 Φ 18,3,0 -105 8 α 4 + 525 64 α 5 D + 105 128 α 6 D 2 3b 1 9/2 + 20b 3 9/2 + 14b 5 9/2 Φ 18,3,1 -105 Φ 18,3,2 525 8 α 4 -525 64 α 5 D -105 128 α 6 D 2 3b 1 9/2 + 20b 3 9/2 + 14b 5 9/2

.55: Semi-major axis functions for the resonant term # 17.

TABLE 3

 3

			.56: Semi-major axis functions for the resonant term # 18.
	# Φ	Semi-major axis function
	19	8 α 2 + 63 8 α 3 D + 9 16 α 4 D 2 b 2 5/2 Φ 19,0 -99 Φ 19,1 99 8 α 2 + 63 8 α 3 D + 9 16 α 4 D 2 b 2 5/2 +
			135 8	α 3 +	15 2	α 4 D +	15 32	α 5 D 2	3b 1 7/2 + 7b 3 7/2
		Φ 19,2 -	405 16	α 3 +	45 4	α 4 D +	45 64	α 5 D 2	2b 1 7/2 + 5b 3 7/2 +
			1365 32	α 4 +	945 64	α 5 D +	105 128	α 6 D 2	3b 0 9/2 + 19b 2 9/2 + 13b 4 9/2
		Φ 19,3	1365 32 α 4 + 945 64 α 5 D + 105 128 α 6 D 2 3b 0 9/2 + 20b 2 9/2 + 14b 4 9/2

TABLE 3

 3

	.57: Semi-major axis functions for the resonant term # 19.

TABLE 3

 3

			.58: Semi-major axis functions for the resonant term # 20.
	# Φ	Semi-major axis function
	21	4 α 2 + 27 4 α 3 D + 3 8 α 4 D 2 b 2 5/2 Φ 21,0 -99 8 α 3 + 75 8 α 4 D + 15 32 α 5 D 2 2b 1 7/2 + 3b 3 7/2 Φ 21,1 -315 Φ 21,2 -1365 8 α 4 + 1155 32 α 5 D + 105 64 α 6 D 2 b 0 9/2 + 4b 2 9/2 + 2b 4 9/2

TABLE 3

 3

		.59: Semi-major axis functions for the resonant term # 21.
	22	Φ 22,0 Φ 22,1 -1905 315 16 α 2 + 27 4 α 3 D + 3 8 α 4 D 2 b 2 5/2 64 α 3 + 45 8 α 4 D + 15 64 α 5 D 2 2b 0 7/2 + 3b 2 7/2 Φ 22,2 -15855 128 α 4 + 1365 64 α 5 D + 105 128 α 6 D 2 5b 1 9/2 + 2b 3 9/2

TABLE 3

 3

	# Φ	Semi-major axis function
	23	Φ 23,0 -81 4 α 2 + 3 4 α 3 D -3 8 α 4 D 2 b 4 5/2 Φ 23,1 105 32 α 5 D 2 2b 3 7/2 + 3b 5 7/2 4 α 3 -15 Φ 23,2 -735 8 α 4 -105 32 α 5 D -105 64 α 6 D 2 b 2 9/2 + 4b 4 9/2 + 2b 6 9/2

.60: Semi-major axis functions for the resonant term # 22.

TABLE 3

 3

.61: Semi-major axis functions for the resonant term # 23.

TABLE 3

 3

	# Φ	Semi-major axis function
	25	Φ 25,0 Φ 25,1 -135 33 4 α 2 + 21 4 α 3 D + 3 8 α 4 D 2 b 2 5/2 8 α 3 + 15 2 α 4 D + 15 32 α 5 D 2 2b 1 7/2 + 3b 3 7/2 Φ 25,2 1365 16 α 4 + 945 32 α 5 D + 105 64 α 6 D 2 b 0 9/2 + 4b 2 9/2 + 2b 4 9/2

.62: Semi-major axis functions for the resonant term # 24.

TABLE 3

 3

			.63: Semi-major axis functions for the resonant term # 25.
	# Φ	Semi-major axis function
	26	Φ 26,0	129 32 855 128	α 2 + α 3 +	21 16 15 8	α 3 D + α 4 D +	3 32 15 α 4 D 2 b 3 5/2 + 128 α 5 D 2 4b 2 7/2 + b 4 7/2
		Φ 26,1 -[	2565 128	α 3 +	45 8	α 4 D +	45 128	α 5 D 2	2b 2 7/2 + 4b 2 7/2 +
			7665 256	α 4 +	945 128	α 5 D +	105 256	α 6 D 2	4b 1 9/2 + 9b 3 9/2 + b 5 9/2 ]
		Φ 20,2	7665 256 α 4 + 945 128 α 5 D + 105 256 α 6 D 2 5b 1 9/2 + 12b 3 9/2 + 2b 5 9/2

TABLE 3

 3

	# Φ	Semi-major axis function
	27	Φ 27,0 -[ 315 16	99 8 α 3 + α 2 + 75 16	27 8 α 4 D + α 3 D + 15 64	3 16 α 5 D 2 α 4 D 2 b 2 5/2 + 4b 1 7/2 + b 3 7/2 ]
		Φ 27,1	945 16	α 3 +	225 16	α 4 D +	45 64	α 5 D 2	2b 1 7/2 + b 3 7/2 +
			1365 16	α 4 +	1155 64	α 5 D +	105 128	α 6 D 2	4b 0 9/2 + 9b 2 9/2 + b 4 9/2
		Φ 27,2 -1365 8 α 4 + 1155 32 α 5 D + 105 64 α 6 D 2 b 0 9/2 + 4b 2 9/2 + 2b 4 9/2

.64: Semi-major axis functions for the resonant term # 26.

TABLE 3 .

 3 65: Semi-major axis functions for the resonant term # 27.

	# Φ	Semi-major axis function
	28	Φ 28,0	315 32 1905 α 2 + 128 α 3 + 33 16 45 α 3 D + 16 α 4 D + 32 3 15 α 4 D 2 b 1 5/2 + 128 α 5 D 2 4b 0 7/2 + b 2 7/2
		Φ 28,1 -[	5715 128	α 3 +	135 16	α 4 D +	45 128	α 5 D 2	2b 0 7/2 + b 2 7/2 +
			15855 256	α 4 +	1365 128	α 5 D +	105 256	α 6 D 2	13b 1 9/2 + b 3 9/2 ]
		Φ 28,2	15855 256	α 4 +	1365 128	α 5 D +	105 256	α 6 D 2	17b 1 9/2 + 2b 3 9/2 +
			167265 1024	α 5 +	6615 256	α 6 D +	945 1024	α 7 D 2	21b 0 11/2 + 20b 2 11/2 + b 4 11/2
		Φ 28,3 -55755 1024 α 5 + 2205 256 α 6 D + 315 1024 α 7 D 2 75b 0 11/2 + 74b 2 11/2 + 5b 4 11/2

TABLE 3 .

 3 66: Semi-major axis functions for the resonant term # 28.

	# Φ	Semi-major axis function
	29	Φ 29,0 Φ 29,1 -65835 12255 512 α 2 + 57 32 α 3 D -45 64 α 4 D 2 + 3 512 α 6 D 4 b 5 5/2 1024 α 3 + 15 128 α 5 D 2 + 15 256 α 6 D 3 + 15 1024 α 7 D 4 b 4 7/2 + b 6 7/2 16 α 4 D -225 Φ 29,2 155855 4096 α 4 -385 256 α 5 D -945 1024 α 6 D 2 + 35 512 α 7 D 3 + 35 4096 α 8 D 4 3b 3 9/2 + 8b 5 9/2 + 3b 7 9/2

TABLE 3

 3

	# Φ	Semi-major axis function
	30	Φ 30,0,0 Φ 30,0,1 -351 81 16 α 2 + 3 16 α 3 D -3 32 α 4 D 2 b 4 5/2 128 α 2 -339 128 α 3 D -177 256 α 4 D 2 + 9 128 α 5 D 3 + 3 256 α 6 D 4 b 4 5/2 Φ 30,0,2 -8613 128 α 2 -213 256 α 4 D 2 + 15 128 α 5 D 3 + 3 256 α 6 D 4 b 4 5/2 128 α 3 D -447 Φ 30,0,3 -81 8 α 2 + 3 8 α 3 D + 3 16 α 4 D 2 b 4 5/2 +
			105 8	α 3 +	15 64	α 5 D 2	9b 3 7/2 + b 5 7/2
		Φ 30,1,0 -105 8 α 3 -15 64 α 5 D 2 b 3 7/2 + b 5 7/2 Φ 30,1,1 15 64 α 3 -645 512 α 5 D 2 + 75 256 α 6 D 3 + 15 512 α 7 D 4 b 3 7/2 + b 5 7/2 64 α 4 D -615 Φ 30,1,2 2625 16 α 3 -825 512 α 5 D 2 + 105 256 α 6 D 3 + 15 512 α 7 D 4 b 3 7/2 + b 5 7/2 64 α 4 D -1785 Φ 30,1,3 105 8 α 3 -15 64 α 5 D 2 11b 3 7/2 + 3b 5 7/2

.67: Semi-major axis functions for the resonant term # 29.

TABLE 3

 3

	.68: Semi-major axis functions for the resonant term # 30.

TABLE 3 .

 3 69: Semi-major axis functions for the resonant term # 31.

	# Φ	Semi-major axis function
	32	Φ 32,0,0 -33 16 α 2 + 21 16 α 3 D -3 32 α 4 D 2 b 2 5/2 Φ 32,0,1 33 128 α 2 -195 128 α 3 D -345 256 α 4 D 2 -33 128 α 5 D 3 -3 256 α 6 D 4 b 2 5/2 Φ 32,0,2 2499 128 α 2 + 1191 128 α 3 D -291 256 α 4 D 2 -39 256 α 6 D 4 b 2 5/2 128 α 5 D 3 -3 Φ 32,0,3 33 8 α 2 + 21 8 α 3 D + 3 16 α 4 D 2 b 2 5/2 +
			135 16	α 3 +	15 4	α 4 D +	15 64	α 5 D 2	9b 1 7/2 + b 3 7/2
		Φ 32,1,0	135 16 α 3 + 15 4 α 4 D + 15 64 α 5 D 2 b 1 7/2 + b 3 7/2
		Φ 32,1,1	405 128 α 3 + 675 64 α 4 D + 2715 512 α 5 D 2 + 195 256 α 6 D 3 + 15 512 α 7 D 4 b 1 7/2 + b 3 7/2
		Φ 32,1,2 -9225 128 α 3 + 1125 64 α 4 D -2625 512 α 5 D 2 -225 256 α 6 D 3 -15 512 α 7 D 4 b 1 7/2 + b 3 7/2 Φ 32,1,3 -135 16 α 3 + 15 4 α 4 D + 15 64 α 5 D 2 11b 1 7/2 + 3b 3 7/2

TABLE 3

 3

	# Φ	Semi-major axis function
	33	Φ 33,0 Φ 33,1 -20955 1575 512 α 2 + 327 64 α 3 D + 225 128 α 4 D 2 + 3 16 α 5 D 3 + 3 512 α 6 D 4 b 1 5/2 1024 α 3 + 345 16 α 4 D + 1485 256 α 5 D 2 + 135 256 α 6 D 3 + 15 1024 α 7 D 4 b 0 7/2 + b 2 7/2 Φ 33,2 100415 4096 α 4 + 9905 512 α 5 D + 2205 512 α 6 D 2 + 175 512 α 7 D 3 + 35 4096 α 8 D 4 11b 1 9/2 + 3b 3 9/2 + 2b 4 9/2

.70: Semi-major axis functions for the resonant term # 32.

TABLE 3

 3

			.71: Semi-major axis functions for the resonant term # 33.
	# Φ	Semi-major axis function
	34	Φ 34,0 Φ 34,1 -5775 128 α 3 -75 8 α 4 D + 75 128 α 5 D 2 b 4 7/2 5775 128 α 3 + 75 8 α 4 D + 75 128 α 5 D 2 b 4 7/2 +
			6405 256	α 4 +	735 128	α 5 D +	105 1024	α 6 D 2	5b 3 9/2 + 9b 5 9/2
		Φ 34,2	2135 256 α 4 + 245 128 α 5 D + 35 256 α 6 D 2 15b 3 9/2 + 28b 5 9/2
		Φ 34,3 -74025 1024 α 5 + 4725 256 α 6 D + 1575 1024 α 7 D 2 3b 2 11/2 + 14b 4 11/2 + 9b 6 11/2

TABLE 3 .

 3 72: Semi-major axis functions for the resonant term # 34.

	# Φ	Semi-major axis function
	35	Φ 35,0 Φ 35,1 -225 8 α 3 + 75 16 α 4 D -75 64 α 5 D 2 b 3 7/2 225 8 α 3 + 75 16 α 4 D + 75 64 α 5 D 2	2b 1 7/2 + 3b 3 7/2 +
			735 32	α 4 +	105 64	α 5 D -	105 128	α 6 D 2	5b 2 9/2 + 9b 4 9/2
		Φ 35,2	4725 256	α 5 -	315 512	α 7 D 2	15b 1 9/2 + 68b 3 9/2 + 43b 5 9/2
		Φ 35,3 -23625 256 α 5 -1575 512 α 7 D 3b 1 11/2 + 14b 3 11/2 + 9b 5 11/2

TABLE 3 .

 3 73: Semi-major axis functions for the resonant term # 35.

	# Φ	Semi-major axis function
	36	Φ 36,0 Φ 36,1 -825 128 α 3 + 75 16 α 4 D + 75 128 α 5 D 2 b 2 7/2 825 128 α 3 + 75 16 α 4 D + 75 128 α 5 D 2 b 2 7/2 +
			1995 256	α 4 +	525 128	α 5 D +	105 256	α 6 D 2	5b 1 9/2 + 9b 3 9/2
		Φ 36,2	665 256	α 4 +	175 128	α 5 D +	35 256	α 6 D 2	15b 1 9/2 + 28b 3 9/2 +
			9135 1024	α 5 +	945 256	α 6 D +	315 1024	α 7 D 2	15b 0 11/2 + 68b 2 11/2 + 43b 4 11/2
		Φ 36,3 -45675 1024 α 5 + 4725 256 α 6 D + 1575 1024 α 7 D 2 3b 0 11/2 + 14b 2 11/2 + 9b 4 11/2

TABLE 3

 3

	# Φ	Semi-major axis function
	37	16 α 3 + 15 4 α 4 D + 15 32 α 5 D 2 b 2 7/2 Φ 37,0 -75 Φ 37,1 315 32 α 4 + 175 32 α 5 D + 35 64 α 6 D 2 3b 1 9/2 + 4b 3 9/2

.74: Semi-major axis functions for the resonant term # 36.

TABLE 3

 3

		.75: Semi-major axis functions for the resonant term # 37.
	# Φ	Semi-major axis function
		Φ 38,0,0 Φ 38,0,1	15 8 α 3 b 1 7/2 45 8 α 3 + 15 4 α 4 D + 15 32 α 5 D 2 b 1 7/2
	38	16 α 4 3b 0 9/2 + 4b 2 9/2 Φ 38,1,0 -35 Φ 38,1,1 4725 128 α 5 + 945 64 α 6 D + 315 256 α 7 D 2 13b 1 11/2 + 5b 3 11/2

TABLE 3 .

 3 76: Semi-major axis functions for the resonant term # 38.

	# Φ	Semi-major axis function
	39	Φ 39,0 -75 16 α 3 + 15 4 α 4 D + 15 32 α 5 D 2 b 0 7/2 Φ 39,1 2205 32 α 4 + 1225 32 α 5 D + 245 64 α 6 D 2 b 1 9/2
		Φ 39,2 -2205 32 α 5 + 945 32 α 6 D + 315 128 α 7 D 2 5b 0 11/2 + 4b 2 11/2

TABLE 3

 3

	# Φ	Semi-major axis function
	40	32 α 3 + 15 16 α 4 D + 15 64 α 5 D 2 b 3 7/2 Φ 40,0 -15 Φ 40,1 105 64 α 4 + 105 64 α 5 D + 35 128 α 6 D 2 3b 2 9/2 + 4b 4 9/2
		Φ 40,2 -945 128 α 5 + 315 64 α 6 D + 315 512 α 7 D 2 3b 1 11/2 + 10b 3 11/2 + 5b 5 11/2

.77: Semi-major axis functions for the resonant term # 39.

TABLE 3 .

 3 78: Semi-major axis functions for the resonant term # 40.

	# Φ	Semi-major axis function
	41	Φ 41,0 Φ 41,1 -315 75 16 α 3 + 15 4 α 4 D + 15 32 α 5 D 2 b 2 7/2 32 α 4 + 175 32 α 5 D + 35 64 α 6 D 2 3b 1 9/2 + 4b 3 9/2 Φ 41,2 2205 64 α 5 + 945 64 α 6 D + 315 256 α 7 D 2 3b 0 11/2 + 10b 2 11/2 + 5b 4 11/2

TABLE 3 .

 3 79: Semi-major axis functions for the resonant term # 41.

	# Φ	Semi-major axis function
	42	256 α 3 + 45 16 α 4 D + 15 64 α 5 D 2 b 1 7/2 Φ 42,0 -225 Φ 42,1 735 64 α 4 + 245 64 α 5 D + 35 128 α 6 D 2 3b 0 9/2 + 4b 2 9/2
		Φ 42,2 -2205 64 α 5 + 315 32 α 6 D + 315 512 α 7 D 2 13b 1 11/2 + 5b 3 11/2

TABLE 3

 3

	# Φ	Semi-major axis function
	43	Φ 43,0 -[	75 32 315 α 3 + 64 α 4 + 15 8 175 α 4 D + 64 α 5 D + 15 64 128 α 5 D 2 b 2 7/2 + 35 α 6 D 2 6b 1 9/2 + b 3 9/2
		Φ 43,1	315 64	α 4 +	175 64	α 5 D +	35 128	α 6 D 2	9b 1 9/2 + 4b 3 9/2 +
			2205 64	α 5 +	945 64	α 6 D +	315 256	α 7 D 2	6b 0 11/2 + 11b 2 11/2 + b 4 11/2
		Φ 43,2 -11025 128 α 5 + 4725 128 α 6 D + 1575 512 α 7 D 2 3b 0 11/2 + 6b 2 11/2 + b 4 11/2

.80: Semi-major axis functions for the resonant term # 42.

TABLE 3 .

 3 81: Semi-major axis functions for the resonant term # 43.

	# Φ	Semi-major axis function
	Φ 44,0,0	15 16 α 3 b 1 7/2 + 35 32 α 4 6b 0 9/2 + b 2 9/2
	Φ 44,0,1	45 16	α 3 +	15 8	α 4 D +	15 64	α 5 D 2 b 1 7/2 +
	44	175 32	α 4 +	175 64	α 5 D +	35 128	α 6 D 2	6b 0 9/2 + b 2 9/2
	32 α 4 9b 0 9/2 + 4b 2 9/2 + 315 64 α 5 17b 1 11/2 + b 3 11/2 Φ 44,1,0 -35 Φ 44,1,1 -175 32 α 4 + 175 64 α 5 D + 35 128 α 6 D 2 9b 0 9/2 + 4b 2 9/2 +
		4725 128	α 5 +	945 64	α 6 D +	315 256	α 7 D 2	17b 1 11/2 + b 3 11/2
	Φ 44,2,0	1575 128 α 5 9b 1 11/2 + b 3 11/2
	Φ 44,2,1	23625 256 α 5 + 4725 128 α 6 D + 1575 512 α 7 D 2 9b 1 11/2 + b 3 11/2

TABLE 3 .

 3 82: Semi-major axis functions for the resonant term # 44.

	# Φ	Semi-major axis function
	45	128 α 5 + 4725 128 α 6 D + 1575 512 α 7 D 2 3b 0 11/2 + 6b 2 11/2 + b 4 11/2 Φ 45,0 -11025 Φ 45,1 315 64 α 4 + 175 64 α 5 D + 35 128 α 6 D 2 9b 1 9/2 + 3b 3 9/2 +
			2205 64	α 5 +	945 64	α 6 D +	315 256	α 7 D 2	6b 0 11/2 + 11b 2 11/2 + b 4 11/2
		Φ 45,2 -11025 128 α 5 + 4725 128 α 6 D + 1575 512 α 7 D 2 3b 0 11/2 + 6b 2 11/2 + b 4 11/2

TABLE 3

 3

			.83: Semi-major axis functions for the resonant term # 45.
	# Φ	Semi-major axis function
	46	128 α 4 D + 5 32 α 5 D 2 + 35 768 α 6 D 3 + 5 1536 α 7 D 4 b 4 7/2 Φ 46,0 -15 Φ 46,1 105 256 α 4 + 385 256 α 5 D + 525 512 α 6 D 2 + 105 512 α 7 D 3 + 35 3072 α 8 D 4 b 3 9/2 + b 5 9/2

TABLE 3

 3

	# Φ	Semi-major axis function
		Φ 47,0,0 Φ 47,0,1 Φ 47,0,2	5 64 α 3 + 5 32 α 4 D + 5 128 α 5 D 2 b 3 7/2 15 64 α 3 + 15 16 α 4 D + 45 64 α 5 D 2 + 5 32 α 6 D 3 + 5 512 α 7 D 4 b 3 7/2 25 128 α 3 + 5 8 α 4 D + 95 256 α 5 D 2 + 25 384 α 6 D 3 + 5 1536 α 7 D 4 b 3 7/2
	47	Φ 47,0,3 -	15 64	α 3 +	15 32	α 4 D +	15 128	α 5 D 2 b 3 7/2 +
				5 128	α 4 +	105 128	α 5 D +	35 256	α 6 D 2	13b 2 9/2 + b 4 9/2
		128 α 4 + 105 128 α 5 D + 35 256 α 6 D 2 b 2 9/2 + b 4 9/2 Φ 47,1,0 -105 128 α 4 + 525 64 α 5 D + 525 128 α 6 D 2 + 175 256 α 7 D 3 + 35 1024 α 8 D 4 b 2 9/2 + b 4 9/2 Φ 47,1,1 -525 256 α 4 + 1225 256 α 5 D + 1015 512 α 6 D 2 + 35 128 α 7 D 3 + 35 3072 α 8 D 4 b 2 9/2 + b 4 9/2 Φ 47,1,2 -735 Φ 47,1,3 105 32 α 4 + 105 32 α 5 D + 35 64 α 6 D 2 4b 2 9/2 + b 4 9/2

.84: Semi-major axis functions for the resonant term # 46.

TABLE 3 .

 3 85: Semi-major axis functions for the resonant term # 47.

	# Φ	Semi-major axis function
	48	Φ 48,0,0 -25 32 α 3 + 5 8 α 4 D + 5 64 α 5 D 2 b 2 7/2 Φ 48,0,1 -75 128 α 3 + 195 128 α 4 D + 235 256 α 5 D 2 + 45 256 α 6 D 3 + 5 512 α 7 D 4 b 2 7/2 64 α 3 + 425 128 α 4 D + 185 128 α 5 D 2 + 55 256 α 6 D 3 + 5 512 α 7 D 4 b 2 7/2 Φ 48,0,2 -125 Φ 48,0,3 75 32 α 3 + 15 8 α 4 D + 15 64 α 5 D 2 b 2 7/2 +
			315 64	α 4 +	175 64	α 5 D +	35 18	α 6 D 2	13b 1 9/2 + b 3 9/2
		Φ 48,1,0	315 64 α 4 + 175 64 α 5 D + 35 128 α 6 D 2 b 1 9/2 + b 3 9/2
		Φ 48,1,1	945 128 α 4 + 1505 128 α 5 D + 1295 256 α 6 D 2 + 385 512 α 7 D 3 + 35 1024 α 8 D 4 b 1 9/2 + b 3 9/2
		Φ 48,1,2	4725 256 α 4 + 5565 256 α 5 D + 3745 512 α 6 D 2 + 455 512 α 7 D 3 + 35 1024 α 8 D 4 b 1 9/2 + b 3 9/2
		Φ 48,1,3 -315 16 α 4 + 175 16 α 5 D + 35 32 α 6 D 2 4b 1 9/2 + b 3 9/2

TABLE 3 .

 3 86: Semi-major axis functions for the resonant term # 48.

	# Φ	Semi-major axis function
		Φ 49,0,0 Φ 49,0,1 -25 75 64 α 3 + 15 32 α 4 D + 5 128 α 5 D 2 b 1 7/2 128 α 3 + 5 8 α 4 D + 95 256 α 5 D 2 + 25 384 α 6 D 3 + 5 1536 α 7 D 4 b 1 7/2 Φ 49,0,2 225 64 α 3 + 75 16 α 4 D + 225 128 α 5 D 2 + 15 64 α 6 D 3 + 5 512 α 7 D 4 b 1 7/2
	49	Φ 49,0,3 -	225 64	α 3 +	45 32	α 4 D +	15 128	α 5 D 2 b 1 7/2 +
				735 128	α 4 +	245 128	α 5 D +	35 256	α 6 D 2	13b 0 9/2 + b 2 9/2
		128 α 4 + 245 128 α 5 D + 35 256 α 6 D 2 b 0 9/2 + b 2 9/2 Φ 49,1,0 -735 256 α 4 + 1225 256 α 5 D + 1015 512 α 6 D 2 + 35 128 α 7 D 3 + 35 3072 α 8 D 4 b 0 9/2 + b 2 9/2 Φ 49,1,1 -735 128 α 4 + 3675 128 α 5 D + 2205 256 α 6 D 2 + 245 256 α 7 D 3 + 35 1024 α 8 D 4 b 0 9/2 + b 2 9/2 Φ 49,1,2 -3675 Φ 49,1,3 735 32 α 4 + 245 32 α 5 D + 35 64 α 6 D 2 4b 0 9/2 + b 2 9/2

TABLE 3 .

 3 87: Semi-major axis functions for the resonant term # 49.

	# Φ	Semi-major axis function
	50	Φ 50,0 -175 128 α 3 + 245 128 α 4 D + 175 256 α 5 D 2 + 65 768 α 6 D 3 + 5 1536 α 7 D 4 b 1 7/2 Φ 50,1 735 32 α 4 + 735 32 α 5 D + 105 16 α 6 D 2 + 175 256 α 7 D 3 + 35 1536 α 8 D 4 b 1 9/2

TABLE 3 .

 3 88: Semi-major axis functions for the resonant term # 50.

	# Φ	Semi-major axis function
	51	64 α 4 + 105 64 α 5 D + 35 128 α 6 D 2 b 2 9/2 Φ 51,0 -105 Φ 51,1 945 256 α 5 + 315 128 α 6 D + 315 1024 α 7 D 2 4b 1 11/2 + 5b 3 11/2

TABLE 3 .

 3 89: Semi-major axis functions for the resonant term # 51.

	# Φ	Semi-major axis function
	52	Φ 52,0 Φ 52,1 -2205 315 32 α 4 + 175 32 α 5 D + 35 64 α 6 D 2 b 1 9/2 128 α 5 + 945 128 α 6 D + 315 512 α 7 D 2 4b 0 11/2 + 5b 2 11/2

TABLE 3 .

 3 90: Semi-major axis functions for the resonant term # 52.

	# Φ	Semi-major axis function
	53	Φ 53,0 -735 64 α 4 + 245 64 α 5 D + 35 128 α 6 D 2 b 0 9/2 Φ 53,1 19845 128 α 5 + 2835 64 α 6 D + 2835 1024 α 7 D 2 b 1 11/2

TABLE 3 .

 3 91: Semi-major axis functions for the resonant term # 53.

	# Φ	Semi-major axis function
	54	Φ 54,0	105 128 945 512	α 4 + α 5 +	105 128 315 256	α 5 D + α 6 D +	35 256 315 α 6 D 2 b 2 9/2 + 2048 α 7 D 2 8b 1 11/2 + b 3 11/2
		Φ 54,1 -	945 512	α 5 +	315 256	α 6 D +	315 2048	α 7 D 2	12b 1 11/2 + 5b 3 11/2 +
			17325 1024	α 6 +	17325 2048	α 7 D +	3465 4096	α 8 D 2	8b 0 13/2 + 13b 2 13/2 + b 4 13/2
		Φ 54,2	17325 4096 α 6 + 17325 2048 α 7 D + 3465 4096 α 8 D 2 10b 0 13/2 + 18b 2 13/2 + 3b 4 13/2

TABLE 3 .

 3 92: Semi-major axis functions for the resonant term # 54.

	# Φ	Semi-major axis function
	55	Φ 55,0		315 64 2205 α 4 + 256 α 5 + 175 64 945 α 5 D + 256 α 6 D + 128 35 315 α 6 D 2 b 1 9/2 + 1024 α 7 D 2 8b 0 11/2 + b 2 11/2
		Φ 55,1 -	2205 256	α 5 +	945 256	α 6 D +	315 1024	α 7 D 2	12b 0 11/2 + 5b 2 11/2 +
				17325 256	α 6 +	24255 1024	α 7 D +	3465 2048	α 8 D 2	21b 1 13/2 + b 3 13/2
		Φ 55,2	17325 256 α 6 + 24255 1024 α 7 D + 3465 2048 α 8 D 2 28b 1 13/2 + 3b 3 13/2

TABLE 3 .

 3 93: Semi-major axis functions for the resonant term # 55.

	# Φ	Semi-major axis function
	56	Φ 56,0		735 128 19845 α 4 + 256 α 5 + 245 128 2835 α 5 D + 128 α 6 D + 35 256	α 6 D 2 b 0 9/2 + 2835 2048 α 7 D 2 b 1 11/2
		Φ 56,1 -	37485 256	α 5 +	5355 128	α 6 D +	5355 2048	α 7 D 2 b 1 11/2 +
				31185 512	α 6 +	31185 2048	α 7 D +	3465 4096	α 8 D 2	13b 0 13/2 + 9b 2 13/2
		Φ 56,2	31185 512 α 6 + 31185 2048 α 7 D + 3465 4096 α 8 D 2 18b 0 13/2 + 13b 2 13/2

TABLE 3 .

 3 94: Semi-major axis functions for the resonant term # 56.

	# Φ	Semi-major axis function
	57	Φ 57,0 Φ 57,1 -315 35 2048 α 5 D + 105 4096 α 6 D 2 + 35 4096 α 7 D 3 + 35 49152 α 8 D 4 b 4 9/2 4096 α 5 + 315 1024 α 6 D + 945 4096 α 7 D 2 + 105 2048 α 8 D 3 + 105 32768 α 9 D 4 b 3 11/2 + b 5 11/2 Φ 57,2 3465 16384 α 6 + 3465 8192 α 7 D + 3465 16384 α 8 D 2 + 1155 32768 α 9 D 3 + 231 131072 α 10 D 4 5b 2 13/2 + 12b 4 13/2 + 5b 6 13/2

TABLE 3

 3

	# Φ	Semi-major axis function
	58	1024 α 4 + 385 1024 α 5 D + 525 2048 α 6 D 2 + 105 2048 α 7 D 3 + 35 12288 α 8 D 4 b 3 9/2 Φ 58,0 -105 Φ 58,1 2205 1024 α 5 + 4095 1024 α 6 D + 945 512 α 7 D 2 + 1155 4096 α 8 D 3 + 105 8192 α 9 D 4 b 2 11/2 + b 4 11/2

.95: Semi-major axis functions for the resonant term # 57.

TABLE 3 .

 3 96: Semi-major axis functions for the resonant term # 58.

	# Φ	Semi-major axis function
	59	Φ 59,0 Φ 59,1 -6615 2205 2048 α 4 + 3675 2048 α 5 D + 3045 4096 α 6 D 2 + 105 1024 α 7 D 3 + 35 8192 α 8 D 4 b 2 9/2 384 α 5 + 945 64 α 6 D + 38745 8192 α 7 D 2 + 2205 4096 α 8 D 3 + 315 16384 α 9 D 4 b 1 11/2 + b 3 11/2

TABLE 3 .

 3 97: Semi-major axis functions for the resonant term # 59.

	# Φ	Semi-major axis function
	60	256 α 4 + 735 256 α 5 D + 105 128 α 6 D 2 + 175 2048 α 7 D 3 + 35 12288 α 8 D 4 b 1 9/2 Φ 60,0 -735 Φ 60,1 6615 256 α 5 + 10395 512 α 6 D + 19845 4096 α 7 D 2 + 1785 4096 α 8 D 3 + 105 8192 α 9 D 4 b 0 11/2 + b 2 11/2

TABLE 3 .

 3 98: Semi-major axis functions for the resonant term # 60.

	# Φ	Semi-major axis function
	61	Φ 61,0 Φ 61,1 -33075 2205 1024 α 4 + 735 512 α 5 D + 315 1024 α 6 D 2 + 105 4096 α 7 D 3 + 35 49152 α 8 D 4 b 0 9/2 1024 α 5 + 4725 256 α 6 D + 14175 4096 α 7 D 2 + 525 2048 α 8 D 3 + 105 16384 α 9 D 4 b 1 11/2

TABLE 3

 3

.99: Semi-major axis functions for the resonant term # 61.

  .4 was found to be, for the Legendre expansion

		7875 2048	+ α 2 4 244 625 65 536	+ α 4 59 594 535 131 072	-s 2 8	1575 64	+ α 2 26 316 675 32 768	+ α 4 297 972 675 32 768
	+s 4 8	111 825 2048	+ α 2 247 037 175 65 536	+ α 4 9 356 341 995 131 072	-s 6 8	51 975 1024	+ α 2 143 468 325 16 384	+ α 4 19 266 061 815 65 536
	+s 8 8	17 325 1024	+ α 2 11 036 025 1024	+ α 4 46 270 899 675 65 536	-s 10 8	α 2 55 180 125 8192	+ α 4 4 197 157 965 4096
	+s 12 8	α 2 55 180 125 32 768	+ α 4 58 036 563 585 65 536	-s 14 8 α 4 55 180 125 32 768	+ s 16 8 α 4 2 749 862 115 32 768
		× α 6 e 4 8 s 4 8 cos (4Ω 8 -4̟ 8 )	(3.44)
	The expansion we use in the model gives
				R sec,27 = cos (4Ω 8 -4̟ 8 ) e 4 8 s 4 8 Φ 27,0 + s 2 8 Φ 27,1 + s 4 8 Φ 27,2	(3.45)

Table 5

 5

	.3

* > 5 • + e < 0.05 2 i * > 5 • 11 i * > 4 • 19 i * > 3

TABLE 5 .

 5 3: Overview of N-Body simulations done with Q set to 100,

TABLE 5 .

 5 4: Overview of semi-analytic simulations done with Q set

  Outcomes in post-resonance elements are similar to simulations given previously and the number of ejections is roughly equal to those in the extended model).

	Q	Number of ejections
	50	0
	100	5
	200	17
	400	19
	800	45
	1500	73
	5000	89

TABLE 6 .

 6 1: Number of ejections as a function of Q using the first semi-analytic model. For each value of Q, 100 simulations have been produced.

  This would give us the mean anomaly after plugging the eccentric anomaly on the right-hand side. But what we need is the contrary. However it is known that Kepler's equation is transcendental, meaning that the reverse relation between E and M doesn't exist explicitly. However many series approximation can be used, among which the Fourier expansion of E -M[START_REF] Murray | Solar System Dynamics[END_REF], which gives

		M = E -e sin(E).	(B.4)
		E = M +	∞	b s (e) sin(sM )	(B.5)
				s=1
	where	b s (e) =	2 s	J s (se)	(B.6)
	and J s the Bessel function of the first kind which can be computed using the
	series (Watson, 1995)
		J s (x) =	∞ m=0	s+2m m!(s + m)! (-1) m x 2	.	(B.7)
	However, as series usually give approximate results, one would try to solve
	Kepler's equation numerically by finding the zero of the function
		f (E) = E -e sin(E) -M	(B.8)
	having plugged in e and M initially. Finally the true anomaly is is related to
	the eccentric anomaly through Gauss' Formula
		f = 2 arctan		1 + e 1 -e	tan	E 2	.	(B.9)

  To complete the calculation we need to express the partial derivatives of R with respect to new variables. To do so, we start by stating that R was first a function of the old set (a, e, i, Ω, ω, M ) so that its differential with respect to time can be written

																												.11)
	Let's compute their time variations. We see from Equation B.10 that
														d̟ dt	=	E na 2 e	∂R ∂e	+	tan(i/2) na 2 E	∂R ∂i
	and							dλ dt	= n -	2 na		∂R ∂a	+	E(1 -E) na 2 e	∂R ∂e	+	tan(i/2) na 2 E	∂R ∂i
		dR dt	=	∂R ∂t	+	∂R ∂a	da dt	+	∂R ∂e	de dt		+	∂R ∂i	di dt		+	∂R ∂Ω	dΩ dt	+	∂R ∂ω	dω dt	+	∂R ∂M	dM dt	(B.12)
	Then we express the same differential stating that R is now a function of (a, e, i, Ω,
	̟, λ)
			dR dt		=	∂R ∂t	+	∂R ∂a	da dt	+		∂R ∂e	de dt		+	∂R ∂i	di dt	+	∂R ∂Ω	dΩ dt	+	∂R ∂̟	d̟ dt	+	∂R ∂λ	dλ dt	(B.13)
	Using the definition (Equation B.11), one can rewrite Equation B.13
	dR dt	=	∂R ∂t	+	∂R ∂a	da dt	+	∂R ∂e	de dt	+	∂R ∂i	di dt	+	dΩ dt		∂R ∂Ω	+	∂R ∂̟	+	∂R ∂λ	+	dω dt	∂R ∂̟	+	∂R ∂λ	+	∂R ∂λ	dM dt
																												(B.14)

  At this point, we can plug derivatives in Equation B.16 into Equation B.19, but we still need to have expressions of partial derivatives of R with to the old elements as a function of partials taken with the new ones. Plugging Equation B.19 into Equation B.18, we have

	dR dt	=	∂R ∂a	da dt	+	de dt			cos(̟)	∂R ∂k	+ sin(̟)	∂R ∂h
								+	di dt	1 2	cos	i 2					cos(Ω)	∂R ∂q	+ sin(Ω)	∂R ∂p
								+	dΩ dt			q	∂R ∂p	-p	∂R ∂q	(B.20)
								+	d̟ dt		k	∂R ∂h	-h	∂R ∂k
								+	∂R ∂λ	dλ dt
	so that after comparing Equation B.13 with Equation B.20, relations between partial
	derivatives are
			∂R ∂e		= cos(̟)	∂R ∂k	+ sin(̟)	∂R ∂h	(B.21a)
	dR dt	=	∂R ∂t	+	∂R ∂a	da dt		+	∂R ∂k	dk dt	+	∂R ∂h	dh dt	+	∂R ∂q	dq dt	+	∂R ∂p	dp dt	+	∂R ∂λ	dλ dt	(B.18)
	Now, we take their derivative with respect to time
				dk dt	=	de dt	cos (̟) -e	d̟ dt	sin (̟)	(B.19a)
				dh dt	=	de dt	sin (̟) + e	d̟ dt	cos (̟)	(B.19b)
					dq dt	=	1 2	di dt		cos	i 2		cos(Ω) -sin	i 2	sin(Ω)	dΩ dt	(B.19c)
					dp dt	=	1 2	di dt		cos	i 2		sin(Ω) + sin	i 2	cos(Ω)	dΩ dt	(B.19d)

https://ssd.jpl.nasa.gov/?sat_elem

Explicit equations of motion are outlined in chapter

A synchronous orbit is when Ωp = n

m is changed to m * and M to mS

Note that for a satellite beneath the synchronous orbit, the right-hand side will have a minus sign, like for Phobos (cf.[START_REF] Efroimsky | Physics of bodily tides in terrestrial planets and the appropriate scales of dynamical evolution[END_REF]

We also call it the one-body problem for the choice of frame centred on the planet.

Harmonics up to J6 have to be included in order to compute the direct dynamics of inner moons, as in Figures 2,

and

in the introduction

and eventually every harmonic : 10λ8 -2λ6, 15λ8 -3λ6 ...

Although J3, J4 and higher harmonics are important for closer satellites

We recall that tides raised on a planet by a satellite tend to increase the satellite eccentricity whereas tides on the satellite decrease it.

The author has repeated in private communications that its theory does not involve the eccentricity of satellites

See Everhart, 1974 for a precise computation

We report a mistake in[START_REF] Everhart | Implicit Single-Sequence Methods for Integrating Orbits[END_REF] on the value of the first node h1 (denoted h2 in the article).

the order of all expansions in that thesis

Acknowledgements

We have h1 = 0.05626 25605 26922 14646 56521 instead of h1 = 0.05626 25605 36922 14646 56521 91032. The eleventh decimal is a 2 instead of a 3.

Part II

The 5:1 Mean Motion Resonance Appendix A

Numerical tools

Since the beginning of the numerical era, many numerical schemes were developed. Celestial mechanics offers a large number of ordinary differential equations to solve. Equations may be solved analytically, only making use of theorems coming from the science of dynamical systems. But since the development of digital computers, the search of generating correct numerical solutions has been growing. Now such solutions are usually produced along with an approximate analytical solution. The goal would be to compute step by step the evolution of dynamical states of a system (such as positions and velocities), given by differential equations coming from physical theories, starting with a set of initial conditions For this thesis, we have been solving equations intensively. For solving the equations of motion, we have used the method of Gauss-Radau, introduced in [START_REF] Everhart | Implicit Single-Sequence Methods for Integrating Orbits[END_REF]Everhart, 1985, then improved in [START_REF] Rein | IAS15: a fast, adaptive, high-order integrator for gravitational dynamics, accurate to machine precision over a billion orbits[END_REF]. Amongst many other numerical methods, this numerical integrator is known to be precise and quick for solving problems from celestial mechanics. Although, many public codes were at our disposition, we have decided to create our own software to do our work. This has not been done from scratch; we have followed the implementation coming from already built codes which we have found or that a sympathetic soul has given us. Reasons for building our own code was for understanding what I was doing and to learn correctly how to implement a code, but also to avoid the uncomfortable use of a "black-boxes" during my thesis. Debugging a code of several thousand lines could be very costly in time. Here was created a software which was compared to the benchmark.

A.1 Gauss-Radau

A.1.1 Presentation

Amongst many numerical integrators for ODEs that was created, the one we have chosen is the one described by Edgar Everhart in Everhart, 1985. We present it here. Laws of Newtonian mechanics make use of second order differential equations. Indeed, Newton's second law relates the acceleration of a body to the external forces applying on it. Therefore, in general, we have an differential equation of the form ẍ = F (x, ẋ, t) (A.1) in which the term on the left side of the equation is the acceleration and the right side term corresponds to the total force acting on the body. Like numerous numerical integration schemes such the types of Runge-Kutta, the Gauss-Radau enables one to compute the dynamical state of the system at a step T starting from an initial state, by evaluating the total force applied on a body at different substeps h n . The property of the Gauss-Radau integrator is that the substeps correspond to nodes in Gaussian quadrature theory. The latter is used to integrate functions precisely. By chosing those substeps one can reach a high order scheme.

First, we expand Equation A.1 in the time variable t, which appears implicitly in the equations of motion. For a 15 th order Gauss-Radau integrating scheme, we expand the acceleration to the seventh power in time

Then integrating with respect to time, we obtain the following expressions concerning the position and the velocity

(A.4) Now, we introduce the coefficients B i = A i T i as well as the substeps h n = tn T . The acceleration at the substeps h n is rewritten as

and the intermediate positions and velocities

and

Now we introduce once more the node h n and a set of coefficients G i in the equation for the acceleration, or the "force"

where we denote r ij = 1 h i -h j . We have expressed ẍ in two ways. Using Equation A.8 and Equation A.5

where coefficients c ij are computed using the following relations .11) Inverse relations are also needed

(A.13)

A.2 Algorithm

The algorithm starts with the initial conditions x 0 and ẋ0 from which we first compute the force

During a step, we only compute F 0 at the beginning of the first sequence. For the first step, the B and G coefficients are set to 0. Then the algorithm alternates between the evaluation of the force at a substep h n and the updating of the B n coefficient.

Explicitly the code will work as follows

• At substep h 1 , the position and velocity, x 1 and ẋ1 , are predicted via Equation A.6 and Equation A.7. Here, for the first step, B coefficients start at zero, but for further steps B will already have attributed values, computed from previous step (see further).

• The force F 1 at substep h 1 is computed using the values of predicted positions and velocities x 1 and ẋ1 .

• G 1 is upgraded from Equation A.9a. Here the previous value of G 1 was zero.

• B 1 is computed from Equation A.10a using the upgraded value of G 1 Explicitly B 1 = G1

Now B 1 has been updated from the evaluation of the force function at h 1 the algorithm continues on the node h 2

• x 2 and ẋ2 are predicted at h 2 with Equation A.6 and Equation A.7. Here we emphasize that computations are made using the upgraded value of B 1 .

• F 2 is calculated using x 2 and ẋ2 .

• G 2 is calculated using Equation A.9b. In here F 2 has just been calculated at the previous and G 1 is the value calculated from the last substep.

• B 1 and B 2 are updated using the updated G 2 in Equation A.10a and Equation A.10b.

Then the same scheme is used for the substeps h 3 , h 4 , h 5 , h 6 and h 7 . At each substep, x i , ẋi are predicted using Equation A.6 and Equation A.7 where values of the set (B j ) are those computed from the previous substep. Then G i is updated from the force F i and then B 1 , B 2 , ..., B i using Equation A.10. Once the force is computed and all B i upgraded at the substep h 7 , one has a complet set B j , j = 1, 7 which will be used to compute the position and the velocity and the time T using the correction equations

The algorithm previously outlined corresponds to one sequence. It is not sure whether the position and the velocity calculated at Equation A.15 and Equation A.16 are yet the final values which we are looking for. In general, another sequence is needed, in which the B values computed in the previous sequence are used to start again the algorithm from h 1 . After the end of the next sequence, another and more precise set of B values are plugged in the correction equations, Equation A.15 Equation A.16.

According to the pionneer Everhart, 1985, six sequences are needed for the first step (initialising B values to zero) and then two are enough to obtain accurate values for x(T ) and ẋ(T ). Here we have followed the indication given in [START_REF] Rein | IAS15: a fast, adaptive, high-order integrator for gravitational dynamics, accurate to machine precision over a billion orbits[END_REF], in which the number of sequences are not predefined but determined during the algorithm. The position and velocity at time T in Equation A.15 and Equation A.16 will get refined sequence after sequence and the sequences are stopped when these values have converged numerically. However, except for very rare cases where a third sequence is needed, two sequences are enough to obtain converged values of x(T ) and ẋ(T ). Also, it is worth mentioning that one would avoid large numbers of force function calls. Counting the evalution of F 0 at the beginning of the first sequence and 7 evaluations inside each sequence, the total number of evaluation per step is

where s denotes here the number of sequences. The Gauss-Radau algorithm is a 15 th order method and one would want to keep the number of function calls down to 15 for this type of integrator. More function evaluations will ruin the benefit of using this specific integrator compared to the use of (for instance) Explicit Runge-Kutta integrators which can be set to be of order n for n force function evaluation per step [START_REF] Stoer | Introduction to numerical analysis[END_REF] Once a first step has been computed and the dynamical state is upgraded to x(T ) and ẋ( T ), a second step can be started. But before evaluating the force at h 1 , one needs to provide good B values to start the first sequence. Going back to Equation A.2 and using another time step T ′ , the acceleration at

Then, defining Q = T ′ T , the ratio between the new time step and the previous one .19) and knowing that

and expanding the parenthesis, starting B values for the next step are given by

Anticipated correction Further refinement can be made before starting a new step. Between the beginning and the end of a step, B values are subjected to some change in order for x(T ) and ẋ(T ) to converge numerically. As explained in Everhart, 1985 and[START_REF] Rein | IAS15: a fast, adaptive, high-order integrator for gravitational dynamics, accurate to machine precision over a billion orbits[END_REF][START_REF] Rein | IAS15: a fast, adaptive, high-order integrator for gravitational dynamics, accurate to machine precision over a billion orbits[END_REF], this change will stay relatively constant from one step to another. Therefore, one can proceed to a further anticipated correction

where B start , are the B values used at the start of a new step , B new are those computed previously from Equation A.21 and B old start -B old end is the change computed between the beginning and the end of the previous step. G values are computed using Equation A.12. Contrary to the first step, the second and subsequent steps will use these relatively good starting values. This will make the convergence of x(T ) and ẋ(T ) much quicker, using usually only two sequences.

Time step

This Gauss-Radau integration scheme is meant to use variable time steps during the whole integration. A chosen time step is correct if it is, a priori, below the lowest time scale of the system. For instance, in the frame of the one-body problem 1 the dynamical time scale corresponds to the orbital period. In order for the integrator to work, one needs a time step which is a fraction of that dynamical scale.

A.3 Gauss-Radau nodes

A.3.1 Derivation

We following here the derivation conducted in [START_REF] Everhart | Implicit Single-Sequence Methods for Integrating Orbits[END_REF]. The position at time t in Equation A.4 was expanded to the order 9 in time. To reach an order 15 integrator, we first expand it to the 15 th power .23) the difference between Equation A.23 and Equation A.4 is

and similarly, the difference in velocity is

We want ∆x = 0 and ∆ ẋ = 0 and such a condition happens if all D i equal zero.

We have 12 equations for 7 unknowns, but we also have some redundancy. Here Coefficients c 7i are actually the elementary symmetric polynomials of h 1 , h 2 , h 3 , h 4 , h 5 , h 6 and h 7 which means that they verify the Vieta's formulas between the roots and the coefficients of a polynomial. Explicitly

The nodes are therefore the roots of the following polynomial

A.4 Numerical implementation

A.4.1 A few changes

A few changes have been made regarding the implementation of other source code [START_REF] Chambers | Mercury -A New Software Package for Orbital Integrations[END_REF], [START_REF] Bolmont | Mercury-T: A new code to study tidally evolving multi-planet systems. Applications to Kepler-62[END_REF]. Any variable involving the value of the Gauss-Radau nodes have been computed beforehand with a precision of 50 digits using Mathematica and the MPFUN 4 arbitrary-precision software for Fortran. This avoids any risk of round-off errors in the prior computation of the c ij and d ij coefficients. Their values have been stored in a separate module as parameters of the system. These high precision numbers permitted us to produce simulations in quadruple precision when needed, instead of double-precision.

A.4.2 Evaluation of the perturbing acceleration

Amongst the flattening force and the tidal acceleration, satellites undergo a third body perturbation from both the Sun and another satellite Equation 2.9

The direct two-body acceleration on the left-hand part of the equation is calculated directly, but we have opted for a more refined implementation to compute the thirdbody acceleration,

The right hand part of the equation can be rewritten

which makes problematic the direct computation of

Indeed, if this condition arises, subtracting 1 from the ratio r 3 ij r 3 j may lead to loss of significant digits, because of the close values of the two terms subtracted [START_REF] Stoer | Introduction to numerical analysis[END_REF]. In general, for two numbers a and b, their subtraction will lead the a high loss a numerical precision when the numbers are very close

where ǫ a-b is the relative error associated to resulting operation. Its expression

shows that when a and b are close to equal, the relative error gets very high and could lead miscalculations in the code.

In order to avoid any loss of precision due to catastrophic cancellations, we present a substitute to the direct computation of the perturbing acceleration which is exposed in [START_REF] Battin | Astronautical Guidance. Electronic sciences series[END_REF]. We first define

4 http://www.davidhbailey.com/dhbpapers/mpfun2015.pdf and the function f so that

The third-body acceleration will then have the form .46) q ij still presents a singularity but one can regularise its expression by making suitable manoeuvres

Then recalling that r ij = r jr i before simplifying, we finally arrive at a regularised expression

From there, we compute f (q ij ) with

and by expanding the numerator, we finally obtain a regularised form for the perturbing acceleration

Lagrange Planetary Equations in Non-singular coordinates

As we have seen in chapter 2, the equations of motion for a perturbed body of mass m around a primary of mass M is Equation 2.17

where R stands for a general disturbing function. For instance, it can represent the action of another body on m, or higher harmonics of the gravitational potential of the primary.

B.1 Generalities on orbital elements

In the absence of perturbation, the trajectory of m will be a conic which can be described with a set of six elements

• a : the semi-major axis. It is half of the long axis (in the case of an ellipse) of the conic and is directly related to the two-body mechanical energy per unit mass (denoted C) of the body orbiting the primary. One has

• e : the eccentricity. It relates to the shape of the conic. A circular orbit will have a null eccentricity, whereas an eccentricity greater than 1 depicts an open orbit, an hyperbola. Then for values between 0 and 1, one has an ellipse, with distinct pericentre and apocentre. The boundary value of 1 stands for a parabolic orbit.

• i : the inclination. The angle between the orbit's normal and the normal to the chosen reference plane.

• Ω : the longitude of the ascending node. Denotes the angular position where the orbital plane crosses the reference plane in an upward way. It is calculated with respect to a reference direction, usually chosen arbitrarily.

• ω : the argument of the pericentre. The angle between the longitude of the ascending node and the pericentre.

• M : the mean anomaly. A fictitious angle related to the position of the body on its orbit. Its variation in time is constant and is given by

Appendix C

Disturbing Function

C.1 Introduction

The disturbing function of a body j acting on i was introduced in chapter 2 in the expression of the perturbed acceleration of a body i Equation 2.17

It represents the third body perturbation of a body j acting on i. We recall here its expression Equation 2.16

and the same expression holds for the perturbation of i acting on j.

Now, before making any expansion, one needs to specify the position of each body. We choose , without lose of generality, that the body j orbits further away than the body i, meaning that as they orbit around Saturn

This would be, in our case, the situation for Titan and Iapetus, or any satellite of Saturn and the Sun. To go further, we rewrite those functions according to Murray and Dermott, 2000

where α is the semi-major axis ratio (Equation 3.27). Here both functions share what we call the direct part

And depending on the respective position of the satellites we have different forms of the indirect part. An inner perturber would be associated with the following internal indirect part

All intermediary summation indexes are defined below. We first have

and then

• m min = 0 if s, j 5 are both even or both odd

• m min = 1 if s, j 5 are neither both even nor both odd

N max denotes the order of the expansion. For the first and simple semi-analytic model, N max was set to 4. But in our extended version, N max was set to higher values in order to obtain a higher expansion in s 8 . All the above series were implemented in Mathematica but a few changes were made in order to compute high orders without adding unnecessary computation time. First, in the inclination function sin(I) and cos(I) were replaced with 2s √ 1s 2 and 1 -2s 2 respectively. For Titan approximations were to make sure that no powers higher than 2 in inclination appeared. Then, the order of eccentricities was set to 4. Therefore, the infinite series in Equation C.28 were replaced with the summation making sure that no higher power than 4 of eccentricity was considered. Then several value of N max were tried.

D.4 Implementation

Now that the mathematics have been given, we will deal with their implementation.

We have seen that any derivative of any Laplace coefficient can be expressed as a function a two principal coefficients. Therefore at the end of the sequence for the b j 1/2 coefficients, b 20 1/2 is miscalculated. Therefore, we decided to compute b j 1/2 coefficients using the series in Equation D.4.