Titan’s rotation

Titan’s rotation

A 3-dimensional theory
B. Noyelles University of Namur - Department of Mathematics - Rempart de la Vierge 8 - 5000 Namur - Belgium IMCCE, Paris Observatory, UPMC, Univ. Lille 1, CNRS UMR 8028 - 77 avenue Denfert Rochereau - 75014 Paris - France    A. Lemaître University of Namur - Department of Mathematics - Rempart de la Vierge 8 - 5000 Namur - Belgium    A. Vienne IMCCE, Paris Observatory, UPMC, Univ. Lille 1, CNRS UMR 8028 - 77 avenue Denfert Rochereau - 75014 Paris - France LAL / University of Lille - 1 impasse de l’Observatoire - 59000 Lille - France
Received / Accepted
Key Words.:
Celestial Mechanics – Planets and satellites: individual: Titan
offprints: B. Noyelles, e-mail: noyelles@imcce.fr



We study the forced rotation of Titan seen as a rigid body at the equilibrium Cassini state, involving the spin-orbit synchronization.


We used both the analytical and the numerical ways. We analytically determined the equilibrium positions and the frequencies of the 3 free librations around it, while a numerical integration associated to frequency analysis gave us a more synthetic, complete theory, where the free solution split from the forced one.


We find a mean obliquity of 2.2 arcmin and the fundamental frequencies of the free librations of about 2.0977, 167.4883, and 306.3360 years. Moreover, we bring out the main role played by Titan’s inclination on its rotation, and we suspect a likely resonance involving Titan’s wobble.


1 Introduction

Since the terrestrial observations of Lemmon et al. (Lemmon93 ()), the rotation of Titan, Saturn’s main satellite, has been assumed to be synchronous or nearly synchronous. This has been confirmed by Lemmon et al. (Lemmon95 ()) and by Richardson et al. (Richardson04 ()) with the help of Voyager I images. In this last work, Titan’s rotation period is estimated at days, whereas its orbital period is days.

The spin-orbit synchronization of a natural satellite is very common in the solar system (such as for the Moon and the Galilean satellites of Jupiter) and is known as a Cassini state. This is an equilibrium state that has probably been reached after a deceleration of the spin of the involved body under dissipative effects, like tides.

Recently, Henrard and Schwanen (Henrard04 ()) have given a 3-dimensional elaborated analytical model of the forced rotation of synchronous triaxial bodies, after studying the librations around the Cassini state. This model has been successfully applied by Henrard on the Galilean satellites Io (Henrard05i ()) and Europa (Henrard05c ()), seen as rigid bodies. Such studies require knowing some parameters of the gravitational field of the involved bodies, which cannot be considered as spheres. Another analytical study has been performed for Mercury by D’Hoedt and Lemaître (DHoedt04 ()), for the case of a spin-orbit resonance.

Since the first fly-bys of Titan by the Cassini spacecraft, we have a first estimation of the useful parameters, more particularly Titan’s and (Tortora et al. Tortora06 ()), so a similar study of Titan’s rotation can be made. In this paper, we propose a study of Titan’s forced rotation, where Titan is seen as a rigid body. The originality of this study over Henrard’s previous studies is that we use both the analytical and the numerical tools and compare our results.

2 Expressing the problem

Titan is here considered as a triaxial rigid body whose principal moments of inertia are written respectively as , , and , with .

2.1 The variables

Our variables and equations have already been used in previous studies; see for instance Henrard and Schwanen (Henrard04 ()) for the general case of synchronous satellites, Henrard (Henrard05i ()) for Io, and Henrard (Henrard05c ()) for Europa.

We consider 3 reference frames : the first is centered on Titan’s mass barycenter and is in translation with the inertial reference frame used to describe the orbital motion of the Saturnian satellites in the TASS1.6 theory (see Vienne & Duriez Vienne95 ()). This is a cartesian coordinate system whose origin is the center of Saturn, and it refers to the equatorial plane of Saturn and the node of this plane with the ecliptic at J2000. The second frame is linked to Titan’s angular momentum, and the third one is rigidly linked to Titan. In this last frame, Titan’s matrix of inertia is written as


We first use Andoyer’s variables (see Andoyer Andoyer26 () and Deprit Deprit67 ()), which are based on two linked sets of Euler’s angles. The first set locates the position of the angular momentum in the first frame , while the second locates the body frame in the second frame tied to the angular momentum (see Fig. 1).

Figure 1: The Andoyer variables (reproduced from Henrard Henrard05i ()).

The canonical set of Andoyer’s variables consists of the three angular variables and their conjugated momenta defined by the norm of the angular momentum and two of its projections:

Unfortunately, these variables present two singularities: when (i.e., the angular momentum is colinear to , there is no wobble), and are undefined, and when (i.e., when Titan’s principal axis of inertia is perpendicular to its orbital plane), and are undefined. That is why we use the modified Andoyer’s variables:

where is Titan’s mean orbital motion , , and . With these new variables, the singularity on has been dropped.

2.2 The free rotation

To describe the dynamics of the system, we should consider the free rotation and the perturbations by other bodies. The Hamiltonian of the free body rotation is also the kinetic energy of the rotation where is the instantaneous rotation vector, the angular momentum vector with respect to the center of mass, and the scalar product of the vector and , where and are respectively defined as




We also deduce from the definitions of the angles and (the wobble):


from which we can easily deduce


and consequently


As a result, the Hamiltonian of the free rotation in the modified Andoyer’s variables is






2.3 Perturbation by Saturn

Considering the parent body Saturn as a point mass , the gravitational potential of the perturbation can be written as


where is the density inside the volume of the body and the distance between the point mass and a volume element inside the body. Using the usual expansion of the potential in spherical harmonics (see for instance Bertotti and Farinella Bertotti90 ()), we find


where and are respectively the longitude and the latitude of Saturn’s barycenter of mass in Titan’s frame, and the distance between this Saturn’s barycenter of mass and the origin of the frame (Titan’s barycenter of mass). If we limit the expansion of (11) to the second order terms and drop the term , which does not produce any effect on the rotation, we have


where , , and are the coordinates of Saturn’s center of mass in Titan’s frame (so we have ). Here, depends on the time since Titan’s motion around Saturn is not circular, but we can introduce , the mean value of , since is Saturn’s mean semimajor axis and its mean eccentricity. (They correspond respectively to Titan’s mean semimajor axis and mean eccentricity in a Saturnian frame.)

We use the formula


coming from the development of (see Brouwer & Clemence Brouwer61 ()):


from which the mean anomaly disappears after averaging. The perturbing potential now reads


and we set , so that we can write






where and are respectively Titan’s mass and radius.

As Henrard (Henrard05c ()) did for Jupiter, we also take Saturn’s oblateness into account. The perturbing potential due to Saturn’s oblateness reads




where is Saturn’s radius, and its .

Finally, the Hamiltonian of the problem reads


3 Analytical study

We intend to use the Hamiltonian (21) to analytically determine the equilibrium position of Titan in the Cassini state related to the spin-orbit synchronization and the 3 frequencies of the free librations around this equilibrium, using the method explained in (Henrard & Schwanen Henrard04 ()). For this analytical study, we consider that Titan has a circular orbit around Saturn, whose inclination on Saturn’s equatorial plane is given by only one periodic term extracted from TASS1.6 ephemerides (Vienne & Duriez Vienne95 ()). This implies that the ascending node of Titan oscillates around a fixed value, so it cannot disappear after averaging the equations. That is why the analytical solutions of (Henrard and Schwanen Henrard04 ()) cannot be used directly, and we first must check that they become unchanged without averaging the ascending node. The true orbital eccentricity of Titan is about , but the opportunity to neglect it will be discussed later, after comparison with the numerical study.

This way, the vector locating Saturn’s barycenter is colinear to






where , and are respectively Titan’s mean inclination, argument of the node and mean longitude in the inertial frame of the ephemerides. (The subscript refers to the fact that Titan is Saturn’s sixth satellite.)

To obtain the coordinates , , and of Saturn is the reference frame bound to Titan , 5 rotations are to be performed:






Table 1 gives the values of the physical and dynamical parameter that we use, and Table 2 gathers the computed values of the corresponding parameters used in the Hamiltonian (21).

Parameters Values References
TASS1.6 Vienne95 ()
TASS1.6 Vienne95 ()
TASS1.6 Vienne95 ()
km IAU 2000 Seidelmann02 ()
Pioneer & Voyager Campbell89 ()
Pioneer & Voyager Campbell89 ()
km IAU 2000 Seidelmann02 ()
Pioneer, Voyager + IERS 2003
Cassini Tortora06 ()
Cassini Tortora06 ()
Note \thetheorem

The mean values of Titan’s mean motion , eccentricity and inclination come from TASS1.6 theory (Vienne & Duriez Vienne95 ()), the radii come from the IAU 2000 recommendations (Seidelmann et al. Seidelmann02 ()), Titan’s mass and Saturn’s come from the Pioneer and Voyager space missions (Campbell & Anderson Campbell89 ()). These two values are those used in TASS1.6 theory, we choose to keep them in order to remain coherent. The mass of Saturn has been derived from the fly-bys of the Pioneer and Voyager space missions, but the published value is given in solar masses. That is why we also indicate IERS 2003 as a reference, which gives us the solar mass. Titan’s J2 and C22 come from the fly-by T11 of the Cassini space mission (Tortora et al. Tortora06 ()), but unfortunately no value for is available yet. We can only hypothesize that it should be included between and , as the case for the Galilean satellites of Jupiter.

Table 1: Physical and dynamical parameters.
Parameter Numerical value
Table 2: Values used in the Hamiltonian (21) that have been computed from the physical and orbital parameters given Table 1.

3.1 Equilibrium

We consider here that the system is exactly at the Cassini state. This implies that:

  • The axis of least inertia, , points to the center of mass of Saturn, so we have , as the mean longitude of Saturn in the frame .

  • The ascending node of the frame (associated to the angular momentum) in the inertial frame has the same precession rate as the ascending node of Saturn in the same inertial frame, i.e. , is the argument of the ascending node of Saturn.

  • There is no wobble, so the angular momentum is colinear with Titan’s axis of highest inertia . This implies , so and .

We have




so it is convenient to introduce this new set of canonical variables:

where represents the angle between the axis of least inertia of Titan and the direction Saturn-Titan, and is the difference between the two ascending nodes. At the exact equilibrium, these two angles should be zero.

This way, and also assuming (i.e., neglecting Titan’s orbital eccentricity), the Hamiltonian (21) becomes


with the mean longitude disappearing after averaging, except of course in the variable. The term has to be added because the canonical transformation we use is time-dependent. The Hamiltonian (31) has been computed with Maple software, and the analytical expressions of the coefficients and are the same as in Henrard and Schwanen (Henrard04 ()). This means that, assuming that the ascending node of the orbit of the considered body circulates or not does not change the expressions of and , the formulae given in Henrard and Schwanen (Henrard04 ()) can be applied to bodies whose node does not circulate, e.g. J-4 Callisto, S-6 Titan or S-8 Iapetus. The analytical expressions of these coefficients are recalled in App.A, while Table 3 gives their numerical values in our context.

Parameter Numerical Value
Table 3: Numerical values of and .

At the exact equilibrium, we have , , , and . These two last equations give






Since Titan’s ascending node oscillates around a fixed value, we have . A numerical resolution of (32) and (33) gives




the asterisk meaning "at the equilibrium".

In the orbital model we use, is constant at rad, which is exactly the value of we get. Such an accuracy of 11 digits is given to indicate the numerical equality of the two values, but it has no real physical meaning, so Titan’s mean obliquity (measured with respect to its orbital inclination) should be nearly zero. This is confirmed by this formula, given by Henrard and Schwanen (Henrard04 ()):


which becomes when the mean value of the precession rate of the line of nodes is zero.

3.2 The fundamental frequencies of the free librations

Since the equilibrium has been found, the Hamiltonian is centered in order to study the behavior of the system near the equilibrium. We introduce a new set of canonical variables:

As a translation, this transformation is canonical. In these variables, the main part of the Hamiltonian of the problem is quadratic. Its quadratic part is named and we have


The analytical expressions of the coefficients and are recalled in App.B, and their numerical values are gathered in Table 4. The reader should be aware that these coefficients are similar to those in Henrard and Schwanen (Henrard04 ()), but different from the ones used by Henrard for Io (Henrard05i ()) and Europa (Henrard05c ()) where other variables are used.

Parameter Numerical Value
Table 4: Numerical values of the coefficients and

We now introduce the following new set of canonical variables:

with and conveniently chosen so as to untangle the variables and . It can be easily checked that this transformation is canonical, because it preserves the differential form; i.e.


With these new variables, the Hamiltonian (39) can be written as




The numerical values of these coefficients are gathered Table 5.

Parameter Numerical Value
Table 5: Numerical values of the coefficients of the Hamiltonian after the variables have been untangled.

We can now introduce the last following set of polar canonical coordinates:



The purpose of this last canonical transformation is to show the free librations around the exact Cassini state. The arguments of these free librations are , , and , and the amplitudes associated are proportional to , , and respectively. We can easily check that this transformation is canonical because we have . We can now write




The numerical results are gathered in Table 6, so the periods associated to the 3 free librations around the equilibrium state are respectively , , and years. Table 7 gives an application of the formulae given in this paper to the Galilean satellites of Jupiter Io and Europa, and we make a comparison with the analytical results of Henrard (Henrard05i () and Henrard052 ()) and the numerical results of Rambaux and Henrard (Rambaux05 ()) obtained with the SONYR model (Rambaux & Bois Rambaux04 ()), which is a relativistic N-body model. The small differences between our results and Henrard’s analytical results come from Henrard neglecting , , , and for Io, and and for Europa.

Proper Modes () (period in years)
Table 6: The free librations around the equilibrium state.
Proper Modes Henrard SONYR this paper
u days days days
v days days days
w days days
u days days days
v years years years
w years years
Note \thetheorem

The results labelled "Henrard" come from (Henrard Henrard05i ()) for Io and (Henrard Henrard052 ()) for Europa, while the results labelled "SONYR" come from (Rambaux & Henrard Rambaux05 ()).

Table 7: Comparison of the periods of the free librations of Io and Europa given by different models.

4 Numerical study

To check the reliability of our previous results and to go further in the study of Titan’s forced rotation, we used the numerical tool. This allowed us first to obtain a solution for the rotation of Titan and then to describe it by frequency analysis and to split the free from the forced solutions.

4.1 Numerical integration

We integrated the 6 equations coming from the Hamiltonian (21) over 9000 years, i.e. between -4500 and 4500 years, the time origin being J1980. In these equations, and come from TASS1.6 ephemerides. We recall that , , and are the coordinates of the barycenter of mass of Saturn in the frame rigidly linked to Titan (see Sect. 2.3).

We used the Adams-Bashforth-Moulton -order predictor-corrector integrator, with a constant timestep year, i.e. day. We considered that the shortest significant fundamental period of the system is given by , i.e. days .

Variable Expression
Note \thetheorem

These conditions have been arbitrarily chosen near the Cassini state, with , , and respectively the values of Titan’s mean longitude, argument of the ascending node, and its instantaneous angular rate, at t=-4500 years, given by TASS1.6. is the initial value of the obliquity on the same date.

Table 8: Initial conditions chosen for the numerical integration, at t=-4500 years.

Table 8 gathers the initial conditions we used. These conditions were arbitrarily chosen near the Cassini state. It implies that, by choosing these initial conditions, we supposed that Titan is at the Cassini state. In these initial conditions, the initial value of is defined as


This equation is very similar to (38). We used it to be sure that the system is near the equilibrium. We did not want to start at the exact equilibrium but very close in order to be able to detect the 3 free librations that we studied in the previous section. However, we should keep in mind that the frequencies computed in Table 6 are in fact limits of the frequencies of the free librations when their amplitudes tend to zero. Thus, too high amplitudes of free librations would alter the frequencies too much. In that way, their comparison with the expected fundamental frequencies of the free librations might be difficult, so their identification as these fundamental frequencies could become doubtful.

(a) (b)
(c) (d)
Figure 2: Numerical simulation of Titan’s obliquity over 9000 years, the time origin being J1980= 2444240 JD. Here the behavior of the variables , , , and is being displayed. Explanations are in the text.

Figure 2 gives plots of some significant data resulting from the numerical integration. Figure 2a shows the behavior of the variables (modulus of the plotted value) and , and Figure 2b shows the behavior of , i.e. the difference between the two nodes. We can see that this angle is oscillating around , as predicted by the theory. We can also visually detect a period of about years, and will see later that it is a forced component due to the behavior of Titan’s orbital ascending node. Figure 2c shows the behavior of the wobble , we obtained it from the variables and . Finally, Figure 2d shows the “obliquity” . In this last panel, we can see the same 700-year-periodic contribution detected in . Unfortunately, looking at these plots does not give information on the free and the forced components of the solutions. That is why we used the frequency analysis technique.

4.2 Analysis of the solutions

We use the frequency analysis to describe the solutions given by the numerical integration, i.e. to give a quasi-periodic representation of these solutions. Such a technique has already been used often to describe the orbital motion of planets (see Laskar Laskar88 ()) or natural satellites (cf. for instance Vienne & Duriez Vienne95 () or very recently Lainey et al. Lainey06 ()).

One of the main difficulties with this kind of problem is that we have two timescales for the periods of the terms that appear in the synthetic representations: Titan’s orbital period is roughly days, while the period of its pericenter is about years. In order to correctly detect the long-period terms, the total time-interval used to analyze the solution should be about as long as the longest period expected, here years, and the timestep shorter than half the shortest period expected (about 5 days). Thus, data over 3200 years should be represented with a timestep of 2.5 days, but this would require about 500000 points. This would take a very long computation time, but fortunately some alternative techniques exist to solve this problem.

The most common technique is the use of a digital filter that splits the short-period terms from the long-period ones (see Carpino et al. Carpino87 ()). However, this technique might alter the signal. Another technique has been used here, which consists in using two samples of data with very close timesteps, as explained in Laskar (Laskar04 ()). More precisely, for each variable, we extracted two samples of 65536 data from the results of the numerical integration, one point every 848 for the first sample and one point every 864 for the other one. As a result, the first sample represents the solutions for years with a timestep days, and the second represents the solutions for years with a timestep days.

These two timesteps are far too large to detect contributions with a period of about days. In fact, the short-period terms are detected, but with a wrong frequency. When a frequency is too high, it is detected as


in analyzing the first sample, and as


in analyzing the second one, where and are (a priori unknown) integers. We have




If we now define as the closest integer to the real (i.e. ), we have


and finally


where Eq.(60) requires that and are close enough, i.e. . In our case, the highest frequency that we can detect with this method is , so we can detect every term with a period longer than days, while analyzing only one sample would give periods longer than about 100 days (i.e. 2 timesteps). Such accuracy is enough to detect Titan’s orbital period.

Table 9 is an example of the decomposition of a variable (here , with the norm of the angular momentum) with the two timesteps. The algorithm used for determining each frequency is taken from (Laskar et al. Laskar92 ()) and has been iteratively applied to refine each frequency, as described in (Champenois Champenois98 ()). In this table, term 1 is a constant part, while the second one is clearly the free libration associated to the proper mode . The slight difference between the obtained and the expected periods should be partly due to the associated amplitude not being null, and partly to the approximations used in our analytical model (i.e. no eccentricity and a constant inclination). However, we can see that the two determinations give very different results for term 3, so we can infer that it is in fact a short-period term.

N Amp. Phase () T (y) Amp. Phase () T (y)
Note \thetheorem

The origin of phases is here the origin of the frequency analysis, i.e. 4499.99344 years before J1980. The series are given in cosine.

Table 9: Decomposition of the solution for with the two timesteps, i.e. days (left) and days (right).

Applying (61) we find days. That is very near to Titan’s orbital period, so this term should be an integer combination of Titan’s orbital period and other contribution(s).

It is now interesting to identify the periodic terms contained in the solutions associated to the considered variables. Table 10 gives the proper modes that are expected. They should appear in the quasiperiodic decompositions of the solutions as parts of integer combinations, so integer combinations of the frequencies of the proper modes are performed to identify each term of the decompositions. We do not use the phases because they are uncertain in the Titan ephemerides given by Vienne & Duriez (Vienne95 ()). The reason is that the given phases are in fact integer combinations of the phases coming from the identified proper modes and very-long-period arguments due to the solar perturbation that are assumed to be constant on an ephemerides-timescale.

Proper Frequency Period Cause
days Rhea
days Titan
days Iapetus
years Sun
Note \thetheorem

The modes to (first part of the Table) come from Vienne & Duriez (Vienne95 ()), while the second part contains the free librations around the Cassini state. These terms have been evaluated from the solutions given by our numerical integration. The fourth column gives the orbital parameter to which the proper mode is linked, being the eccentricity of the satellite , and the sine of its semiinclination. The subscripts are for Rhea, for Titan, and for Iapetus.

Table 10: Proper modes of the system.
N Amp. Phase () T (y) Ident. Cause
1 constant
3 days unknown
Note \thetheorem

The series are in cosine. The fourth column gives the orbital parameters associated to each identified term.

Table 11: Quasiperiodic decomposition of the variable .
N Amp. Phase () T (y) Ident. Cause
1 constant
8 Sun
Table 12: Quasiperiodic decomposition of the variable . The series are in cosine.
N Amp. Phase () T (y) Ident. Cause
Table 13: Quasiperiodic decomposition of the complex variable .
N Amp. Phase () T (y) Ident. Cause
2 days unknown
5 days unknown
8 days unknown
Table 14: Quasiperiodic decomposition of the variable . The series are in sine.
N Amp. Phase () T (y) Ident. Cause