New insights on Titan’s interior from its obliquity
Abstract
We constructed a 6degrees of freedom rotational model of Titan as a 3layer body consisting of a rigid core, a fluid global ocean, and a floating ice shell. The ice shell exhibits partiallycompensated lateral thickness variations in order to simultaneously match the observed degreetwo gravity and shape coefficients. The rotational dynamics are affected by the gravitational torque of Saturn, the gravitational coupling between the inner core and the shell, and the pressure coupling at the fluidsolid boundaries. Between and of our model Titans have an obliquity (due to a resonance with the year periodic annual forcing) that is consistent with the observed value.
The shells of the successful models have a mean thickness of to km, and an ocean of 250 km thickness. Our simulations of the obliquity evolution show that the Cassini obliquity measurement is an instantaneous one, and does not represent a mean value. Future measurements of the time derivative of the obliquity would help to refine the interior models. We expect in particular a variation of roughly 7 arcmin over the duration of the Cassini mission.
1 Introduction
The Cassini spacecraft, in orbit around Saturn since July 2004, has allowed huge progress on modelling of the internal structure and the rotational dynamics of Titan. An internal ocean is consistent with the measurements of the tidal Love number (Iess et al., 2012) and was theoretically predicted by Lunine & Stevenson (1987), this prediction being supported by several following studies, e.g. (Grasset & Sotin, 1996; Grasset et al., 2000; Tobie et al., 2005; Fortes et al., 2007). A comparison between the shape of Titan (Zebker et al., 2009) (Tab.1) and its gravity field (Iess et al., 2010) (Tab.2) suggests either variations in the thickness of a floating ice shell (Nimmo & Bills, 2010; Hemingway et al., 2013) or lateral variations in the shell’s density (Choukroun & Sotin, 2012).
Parameter  Value 

Subplanetary equatorial radius  km 
Along orbit equatorial radius  km 
Polar radius  km 
Mean radius  km 
SOL1  SOL2  

–  
Cassini observed Titan’s rotation as well. The most recent measurements suggest the expected synchronous rotation (Meriggiola & Iess, 2012) and a pretty high obliquity of at the mean date March 2007, already detected by (Stiles et al., 2008). If we assume that the rotation of Titan has reached its most probable dynamical equilibrium state, i.e. Cassini State 1, then this obliquity is not consistent with a rigid Titan (Noyelles et al., 2008; Bills & Nimmo, 2008, 2011). However, the presence of an internal ocean can lead to a resonant process raising the obliquity of Titan (Baland et al., 2011), making the high obliquity a possible signature of a global subsurface ocean.
In this paper, we simulate the rotation of Titan, considering both the internal structure and all the dynamical degrees of freedom. Our Titan is a 3layer body composed of a rigid inner core, a global ocean and rigid shell with a variable thickness. For each of the 2 rigid layers, we simulate at the same time the longitudinal motion, the orientation of the angular momentum, and of the figure polar axis. The dynamics of these 2 layers will be affected by the gravitational pull of Saturn, the pressure coupling at the interface with the ocean and the gravitational coupling between them. The pressure coupling is modelled after Baland et al. (2011) and the gravitational coupling after Szeto & Xu (1997). In calculating the torques, we take into account variations in the thickness of the ice shell (Nimmo & Bills, 2010) consistent with the gravity and topography constraints. We then identify interior structures for which the predicted rotation state is consistent with the observations, before simulating the expected behavior of the obliquity of Titan.
Our model confirms the conclusion of Baland et al. (2011) that the unexpectedly high obliquity of Titan could be due to a resonance with the periodic annual forcing. We go further, however, in showing that the obliquity is predicted to be timevariable (Fig 9): a prediction which analysis of Cassini radar observations (Bills et al., 2013) should be able to test.
2 The equations of the problem
The approach that we follow below is a generalization of the scheme adopted by Baland et al. (2011). There are two important innovations in our approach. First, we consider the threedimensional orientation of the shell and of the core, so that we can simultaneously treat both obliquity (Bills & Nimmo, 2011; Baland et al., 2011) and also longitudinal librations (Van Hoolst et al., 2013; Richard et al., 2014), as well as determining the magnitude of the usually neglected polar motion. Second, we explicitly take into account the rigidity and spatial variations in thickness of the ice shell, which are indicated by Titan’s topography and gravity (Nimmo & Bills, 2010; Hemingway et al., 2013) and which affect the resulting torques.
Although our model is quite complicated, it does neglect some potentially significant effects. Most notably, we do not consider the effects of either atmospheric torques or torques due to flow in the subsurface ocean. The low viscosity of water argues against the latter being important, but in some cases tidal forcing can lead to strong flows (Noir et al., 2009; Cebron et al., 2012). We defer consideration of this topic to future work. As discussed below, we also neglect the potential effect of restoring torques due to elastic deformation of the ice shell (Goldreich & Mitchell, 2010; Richard et al., 2014).
2.1 Parameterization of the problem
We simulate the orientation of both the rigid inner core and the rigid shell (or crust). For that, we need 2 sets of Euler angles, respectively for the core and for the shell to represent the orientation of the principal axes of inertia of the considered layer in an inertial reference frame (Fig.1).
These Euler angles present a virtual singularity. If the quantity is null, then the angles and are not uniquely defined, but their sum is. In practice, it appears that when is small enough, then numerical uncertainties can erronously suggest an erratic behavior. We bypassed this problem by using the following cartesianlike coordinates:
When is null, then and are both null and the system does not present any singularity.
We also need to represent the angular momentum of each of these rigid layers. For that we use as variables the components of the associated rotation vector ; this yields for the shell and for the core. , and are the principal moments of inertia of the layer under consideration; we have for the core:
(1)  
(2)  
(3) 
and similar formulae for the shell. is the density of the core, being the classical writings for the cartesian coordinates, in the reference frame of the principal axes of inertia .
2.2 Kinematic equations controlling the Euler angles
To determine the relations linking the Euler angles of a layer to the components of the angular momentum, we have to keep in mind the geometry of the problem (Fig.1). As explained for instance in (Fowles & Cassiday, 1999), the rotation vector represents 3 successive rotations:

a rotation around the zaxis (here ) of an angle ,

then a rotation around the new, but not final, xaxis of angle ,

and finally a rotation around the final zaxis, i.e. , of an angle .
This reads mathematically:
(4) 
with
(5) 
and
(6) 
This yields
(7)  
(8)  
(9) 
and
(10)  
(11)  
(12) 
The equations (10) and (12) illustrate the virtual singularity we mentioned above. These formulae are the same as the ones present in (Black et al., 1995; Harbison et al., 2011) but are given with a different sign in numerous other studies, e.g. (Williams et al., 2001), probably because of a different sign convention. We now get straightforwardly:
(13)  
(14)  
(15) 
The virtual singularity has nearly disappeared. The only numerical problem that could remain is due to in the Eq.(13) and (14), so we replace it by its Taylor expansion for .
As shown in (Baland et al., 2011), the dynamical equations reduce to
(16)  
(17) 
the relevant torques being:

: gravitational torque of Saturn on the core,

: pressure coupling of the ocean at the coreocean boundary,

: gravitational torque of the shell on the core,

: gravitational torque of Saturn on the shell,

: pressure coupling of the ocean at the shellocean boundary,

: gravitational torque of the core on the shell.
Baland et al. (2011) have shown that the sum of the torques acting on the ocean is null if the fluid is in hydrostatic equilibrium.
Since we work in the noninertial reference frames of the principal axes of inertia of the core , and of the shell , we must add in Eq.(16) and in Eq.(17).
We now detail the torques affecting the rotation.
2.3 The gravitational pull of Saturn
For a rigid triaxial body whose principal moments of inertia are , the gravitational torque of a perturber is classically given by:
(18) 
where is the gravitational constant, the mass of the perturber, and the vector locates the perturber from the center of mass of the triaxial body. In our case, this is the TitanSaturn vector. A derivation of this torque is proposed in (Murray & Dermott, 2000), inspired from (MacMillan, 1936; Ramsey, 1937, 1940).
We here treat the inner core and the shell as 2 independent triaxial rigid bodies, and we have:
(19)  
(20) 
We use the TASS1.6 ephemerides (Vienne & Duriez, 1995) to express the TitanSaturn vector . These ephemerides are presented under a quasiperiodic form, in which every sinusoidal contribution due to any perturber like Saturn’s oblateness or the Solar attraction is explicitly expressed. They are given in the inertial reference frame defined by the equatorial plane of Saturn at J2000 and the node of this plane with the ecliptic at the same date. As an example, in this reference frame, the quantities related to the inclination of Titan and its ascending (Titan is denoted as S6 Titan) are
(21) 
the quantities and being given in the Tab.3.
i  Amplitude  Phase  Period  

(rad/y)  (years)  
0  0  –  
1  
2  
3  
4 
Similar tables exist for the mean motion, the mean longitude, and the quantity relevant to the eccentricity and the longitude of the pericenter , i.e. .
The TASS1.6 ephemerides are probably not the most accurate we could get, but they have the huge advantage of being presented in a quasiperiodic form giving explicitly the different contributions affecting Titan’s orbit. From the periods of the sinusoidal quantities indexed from 1 to 4 in the Eq.(21) and the Tab.3, we can say that corresponds to the motion of the ascending node induced by the oblateness of Saturn, is a perturbation by Iapetus, and correspond to the Solar orbital perturbation, the orbital period of Saturn around the Sun being years. It is impossible to have such a decomposition with the JPL HORIZON ephemerides since they are given over a too short timespan ( years) with respect to the relevant periods. This is why we choose to use TASS1.6, Baland et al. (2011) having made the same choice.
The choice of the reference frame is not straightforward. To get an obliquity having a straightforward physical meaning, it is often advisable to use the Laplace Plane, that minimizes the variations of the inclination. Unfortunately, there are in the literature several inconsistent definitions of the Laplace Plane, since there are several ways to minimize the variations of the inclination (over which time interval should we minimize? how do we measure the variations of the inclination?…). Noyelles (2009) suggests using the constant term in the quasiperiodic decomposition of the inclination to define the inertial reference frame. This is a kind of averaging of the orbital plane, that is very close to the Laplace Plane. In particular, this choice avoids a problem of apparent erratic behavior of the rotation pole that could happen for a satellite orbiting far off its parent planet, when the rings’ plane is chosen as the reference plane. The reason is that these satellites usually have a significantly inclined orbit because of the Solar perturbation.
This is why our reference frame is obtained from the reference frame of the ephemerides after 2 rotations: a rotation of , and a rotation of , these numbers being derived from the line in the Tab.3. A very easy way to implement these rotations is just to drop from Eq.(21).
After a straightforward calculation we get from TASS1.6 the TitanSaturn vector in the inertial reference frame that we just defined, and then the vector after 3 rotations of the Euler angles:
(22)  
(23) 
At this time we can simulate the rotation of the inner core and the rigid shell as two independent rigid bodies. We successfully checked our code in comparing with the results given by the numerical code based on an Hamiltonian formalism that we used in several previous studies (Noyelles et al., 2008; Noyelles, 2009, 2010; Noyelles et al., 2011).
2.4 The gravitational coupling between the inner core and the shell
This subsection is based on Szeto & Xu (1997), in which the gravitational coupling between the shell and the core of the Earth is estimated. This gravitational coupling is due to the misalignment of the principal axes of inertia of these 2 layers. We now need to be more specific as to their structure:

The inner core is a triaxial ellipsoid with a constant density , its radii being denoted . Its surface is the coreocean boundary.

The shell (or crust) has a constant density as well. Its shape can be described by two concentric and coaxial triaxial ellipsoids. The radii of the outer one are denoted , they correspond to the observed shape of Titan (Zebker et al., 2009), while the inner one is aligned with the outer one, its radii are denoted . This inner edge of the crust is the shellocean boundary. It is important to define ,, as well as , and to allow for the possibility of lateral shell thickness variations, which Nimmo & Bills (2010) argued are likely to exist.
Following Szeto & Xu (1997), the gravitational torque acting on the inner core due to the shell reads:
(24) 
where points in this subsection to the position of a mass element of the core, and is the potential of the shell given by (Szeto & Xu, 1997):
(25) 
where

is a constant that does not affect the result,

,

,

is the polar flattening of the shellocean boundary,

is the polar flattening of the surface of Titan,

is the equatorial ellipticy of the shellocean boundary,

is the equatorial ellipticity of the surface of Titan,

is a Legendre polynomial,

is a Legendre associated function,

and are respectively the colatitude and the east longitude of the mass element involved, in the reference frame of the principal axes of inertia of the shell .
After some algebra (see App.A) we get
(26) 
where , and () are the elements of the transition matrix between the coordinates in the reference frame of the shell and the ones in the reference frame of the core , i.e.
(27) 
An analogous calculation gives us the torque of the inner core acting on the shell:
(28) 
2.5 Influence of the ocean
Baland et al. (2011) have shown that if the ocean is in hydrostatic equilibrium, then the pressure torque can be expressed as additional terms in the gravitational torques acting on the two rigid layers. We can split the ocean into two parts: a top and a bottom one, their boundary being spherical. The resulting pressure torque integrated over a spherical boundary is null, so the radius of this sphere has no influence on the results (an outcome we verified numerically).
We call , , , , and the principal moments of inertia, respectively of the top ocean in the reference frame of the shell and of the bottom ocean in the reference frame of the core . Instead of writing the contribution of the ocean as independent torques and , it is more appropriate to alter the other torques as , , , and . And we have:
(29) 
(30) 
(31) 
and
(32) 
with
(33)  
(34) 
being the constant density of the ocean.
We have now the whole equations of the problem, consisting of 12 variables , , , , , , , , , , , , the first 6 describing the orientation of the core, and the last 6 the orientation of the shell. The components of the rotation vector in the reference frame of the principal axes of inertia of the considered layer are obtained from the angular momentun and division by the appropriate moment of inertia.
As a summary, we here gather these equations. We have:
(35)  
(36)  
(37)  
(38)  
(39)  
(40)  
for the core, and
(41)  
(42)  
(43)  
(44)  
(45)  
(46)  
for the shell.
3 A numerical solution
A numerical solution of the equations (35) to (46) is here appropriate since we want to include complete ephemerides and 6 dynamical degrees of freedom.
3.1 Numerical integration of the equations
The numerical integrations are performed with the AdamsBashforthMoulton 10th order predictor corrector integrator (see e.g. (Hairer et al., 2003)), with a tolerance of and a step size of day. This corresponds to of the orbital period of Titan.
The rotation of Titan is expected to be at a dynamical equilibrium. Such equilibriums are known as Cassini States (Cassini, 1693; Colombo, 1966; Peale, 1969) for rigid bodies. The expected state for the natural satellites of the giant planets is Cassini State 1 since it is the most stable. In our case of a 3layer Titan, we initially assume an analogous state in which the inner core and the crust are close to the location of Cassini State 1 if they were only interacting with Saturn. This state corresponds to a synchronous rotation, a small obliquity and a small polar motion. As a consequence, the angular momentum of the shell and the core should approximate and where is the mean motion, or orbital frequency, of Titan, and the spin angles of the core and of the shell should be always close to the orbital mean longitude of Titan .
Because of the effects that are neglected in the theory of the Cassini States, especially the couplings between the different layers and the perturbations considered in the orbital motion of Titan, it is very difficult to derive analytically the optimal initial conditions for the Euler angles and the rotation vector. In practice, our initial conditions are usually close enough to the optimum state that oscillations round the equilibrium result, these free oscillations having an arbitrary amplitude due to the choice of the initial conditions, and a proper frequency whose value depends on the parameters of the system, here the interior of Titan. These free oscillations pollute the analysis of the solutions in acting as a noise, for this reason we wish their amplitude to be as small as possible. For that, we refine numerically the initial conditions thanks to an iterative algorithm based on the frequency analysis.
The basic idea is that since the orbital motion of Titan can be given under a quasiperiodic form, and that the rotation of Titan is not expected to be chaotic, then this rotation can be expressed under a quasiperiodic form as well. A complex variable of the problem that does not diverge can read as a sum of a converging trigonometric series like
(47) 
where are constant complex amplitudes, and constant frequencies, with
(48) 
the bullet meaning that the coefficients have been numerically determined. A detailed description of the algorithm is given in Appendix B. In the case of a real variable, Eq.48 becomes
(49) 
or
(50) 
where the amplitudes are now real, and the are real phases, previously included in the complex amplitudes in Eq.48.
The frequency analysis algorithm we use is based on NAFF (see (Laskar, 1993) for the method, and (Laskar, 2005) for the convergence proofs), with a refinement suggested by (Champenois, 1998) consisting in iterating the process to improve the accuracy of the determination. The frequencies have 2 origins: they might be either forcing frequencies, present in the orbital motion of Titan, or free frequencies, due to the departure from the exact equilibrium. The amplitude associated with the latter should be as small as possible. To get the appropriate initial conditions we use an iterative algorithm (Noyelles et al., 2014), consisting in:

A first numerical integration of the equations of the system, with initial conditions conveniently chosen,

Frequency analysis of the solution and identification of the contributions depending on the free modes,

Evaluation of the free modes at the origin time of the numerical simulation, and removal from the initial conditions,
then the process is iterated until convergence. This algorithm has already been successfully applied in problem of rotational dynamics (Dufey et al., 2009; Noyelles, 2009; Robutel et al., 2011), in dynamics of exoplanetary systems (Couetdic et al., 2010), and in the analysis of groundtrack resonances around Vesta (Delsate, 2011).
We know of at least 3 alternative methods to reduce the amplitude of the free librations:

Bois & Rambaux (2007) propose to fit the mean initial conditions in order to locate the spinorbit system at its center of libration,

Peale et al. (2007) add a damping in the equations that reduces the amplitude of the free librations. The damping must be slow enough, i.e. adiabatic, to not alter significantly the location of the equilibrium,

Yseboodt & Margot (2006), in the framework of a numerical integration of the spin and of the orbit of Mercury, start from a simple SunMercury system in which the equilibrium is very easy to determine analytically, and slowly switch on the planetary perturbations in order to create an adiabatic devitation of the equilibrium without creation of any free libration. In our case, this would require us to simultaneously integrate the orbit of Titan, rather than use the existing ephemerides.
All of these methods, including ours, give accurate results when appropriately used.
In this study, we simulated the rotation of thousands of model Titans (see Sec.4) and did not refine the initial conditions for all of them. In practice, we did it for just a few of them, and got initial conditions that we considered to be good enough for the remainder.
3.2 Outputs
Our set of variables describes all the dynamical degrees of freedom of the core and the shell, so we are able to express any observable of the rotation. Our outputs are, for these two rigid layers:

the longitudinal librations,

the obliquity,

and the polar motion.
The longitudinal librations of the shell are usually considered as the most significant output since they can reveal a global fluid layer (Van Hoolst et al., 2008). There are at least two ways to define them: the tidal librations and the physical librations. The tidal librations represent the longitudinal misalignment between the directions of the long axis of the layer under consideration (the shell or the core) and the SaturnTitan direction. The physical significance of these librations is that they control the amount of tidal stress and heating arising. They are given by
(51)  
(52) 
where is the unit vector tangent to the trajectory of Titan around Saturn:
(53) 
and the unit vector normal to the orbit:
(54) 
and being respectively the position vector of Titan, and its velocity.
The physical librations