Titan’s rotation
Key Words.:
Celestial Mechanics – Planets and satellites: individual: TitanAbstract
Context:
Aims:
We study the forced rotation of Titan seen as a rigid body at the equilibrium Cassini state, involving the spinorbit synchronization.
Methods:
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.
Results:
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.
Conclusions:
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 spinorbit 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 3dimensional 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 spinorbit resonance.
Since the first flybys 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
(1) 
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).
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
(2) 
and
(3) 
We also deduce from the definitions of the angles and (the wobble):
(4) 
from which we can easily deduce
(5) 
and consequently
(6) 
As a result, the Hamiltonian of the free rotation in the modified Andoyer’s variables is
(7) 
with
(8) 
and
(9) 
2.3 Perturbation by Saturn
Considering the parent body Saturn as a point mass , the gravitational potential of the perturbation can be written as
(10) 
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
(11) 
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
(12) 
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
(13) 
coming from the development of (see Brouwer & Clemence Brouwer61 ()):
(14) 
from which the mean anomaly disappears after averaging. The perturbing potential now reads
(15) 
and we set , so that we can write
(16) 
with
(17) 
and
(18) 
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
(19) 
with
(20) 
where is Saturn’s radius, and its .
Finally, the Hamiltonian of the problem reads
(21) 
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 spinorbit 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
(22) 
with
(23) 
(24) 
and
(25) 
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:
(26) 
with
(27) 
and
(28) 
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 flybys 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 flyby 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.
Parameter  Numerical value 

km  
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
(29) 
and
(30) 
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 SaturnTitan, 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
(31) 
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 timedependent. 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. J4 Callisto, S6 Titan or S8 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 

At the exact equilibrium, we have , , , and . These two last equations give
(32) 
and
(33) 
with
(34) 
Since Titan’s ascending node oscillates around a fixed value, we have . A numerical resolution of (32) and (33) gives
(35) 
(36) 
hence,
(37) 
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 ()):
(38) 
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
(39) 
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 

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.
(40) 
With these new variables, the Hamiltonian (39) can be written as
(41) 
with
(42) 
(43) 
(44) 
(45) 
(46) 
(47) 
The numerical values of these coefficients are gathered Table 5.
Parameter  Numerical Value 

We can now introduce the last following set of polar canonical coordinates:
with
(48) 
(49) 
(50) 
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
(51) 
with
(52) 
(53) 
(54) 
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 Nbody 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) 

u  
v  
w 
Proper Modes  Henrard  SONYR  this paper 

Io  
u  days  days  days 
v  days  days  days 
w  days  days  
Europa  
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 ()).
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 AdamsBashforthMoulton order predictorcorrector 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 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
(55) 
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 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 700yearperiodic 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 quasiperiodic 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 longperiod terms, the total timeinterval 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 shortperiod terms from the longperiod 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 shortperiod terms are detected, but with a wrong frequency. When a frequency is too high, it is detected as
(56) 
in analyzing the first sample, and as
(57) 
in analyzing the second one, where and are (a priori unknown) integers. We have
(58) 
and
(59) 
If we now define as the closest integer to the real (i.e. ), we have
(60) 
and finally
(61) 
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 shortperiod term.
N  Amp.  Phase ()  T (y)  Amp.  Phase ()  T (y) 

1  
2  
3 
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.
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 verylongperiod arguments due to the solar perturbation that are assumed to be constant on an ephemeridestimescale.
Proper  Frequency  Period  Cause 
Mode  
days  Rhea  
days  Titan  
days  Iapetus  
years  
years  
years  
years  
years  
years  
years  Sun  
years  
years  
years 
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.
N  Amp.  Phase ()  T (y)  Ident.  Cause 

1  constant  
2  
3  days  unknown 
Note \thetheorem
The series are in cosine. The fourth column gives the orbital parameters associated to each identified term.
N  Amp.  Phase ()  T (y)  Ident.  Cause 

1  constant  
2  
3  
4  
5  
6  
7  
8  Sun 
N  Amp.  Phase ()  T (y)  Ident.  Cause 

1  
2  
3  
4  
5  
6 
N  Amp.  Phase ()  T (y)  Ident.  Cause 

1  
2  days  unknown  
3  
4  
5  days  unknown  
6  
7  
8  days  unknown  
9 
N  Amp.  Phase ()  T (y)  Ident.  Cause 

1  
2  
3  
4  
5  
6 