Librational response of a deformed 3layer Titan perturbed by nonkeplerian orbit and atmospheric couplings
Abstract
The analyses of Titan’s gravity field obtained by Cassini space mission suggest the presence of an internal ocean beneath its icy surface. The characterization of the geophysical parameters of the icy shell and the ocean is important to constrain the evolution models of Titan. The knowledge of the librations, that are periodic oscillations around a uniform rotational motion, can bring piece of information on the interior parameters.
The objective of this paper is to study the librational response in longitude from an analytical approach for Titan composed of a deep atmosphere, an elastic icy shell, an internal ocean, and an elastic rocky core perturbed by the gravitational interactions with Saturn. We start from the librational equations developed for a rigid satellite in synchronous spinorbit resonance. We introduce explicitly the atmospheric torque acting on the surface computed from the Titan IPSL GCM (Institut Pierre Simon Laplace General Circulation Model) and the periodic deformations of elastic solid layers due to the tides. We investigate the librational response for various interior models in order to compare and to identify the influence of the geophysical parameters and the impact of the elasticity.
The main librations arise at two wellseparated forcing frequency ranges: low forcing frequencies dominated by the Saturnian annual and semiannual frequencies, and a high forcing frequency regime dominated by Titan’s orbital frequency around Saturn. At low forcing frequency, the librational response is dominated by the Saturnian gravitational torque and the atmospheric torque has a small effect. In addition, the libration amplitude in that case is almost equal to the magnitude of the perturbation. The modulation of the gravitational torque amplitude at the orbital frequency with periodic deformation induces longperiod terms in the librational response which contain information on the internal structure. At high forcing frequency the libration depends on the inertia of the layers and the elasticity can strongly reduce its amplitude at orbital frequency. For example, the amplitude of diurnal libration for oceanic models goes from about meters if the icy shell is purely rigid to meters when the elasticity is included, i.e. a reduction of about . For models without ocean, diurnal libration goes from meters in a rigid case to meters for an elastic case, a very low reduction due to the weak deformation of an entirely solid satellite compared to the deformation of a thin icy shell. Oceanic models with elastic solid layers have the same order of libration amplitude than the oceanless models, which makes more challenging to differentiate them by the interpretation of librational motion.
keywords:
Titan, libration, elasticity, dynamic, orbital perturbation, atmospheric torquecomma,sort
Introduction
Titan exhibits a rich interior structure due to its large mean radius of 2575 kilometers. The recent measurements of the gravity field by Iess et al. (2010, 2012) reveal that Titan’s moment of inertia (MoI) is as low as . A large panel of internal structures, made of a low density core surrounded by icy layers, are consistent with this range of MoI and the mean density value of kg m deduced by Iess et al. (2010). In addition, the presence of a global internal ocean has been suggested from several techniques. First, Lorenz et al. (2008) deduced the putative internal ocean through the response of the rotational motion of the surface to the atmospheric coupling. A second piece of evidence has been obtained from the determination of the nonzero obliquity by Stiles et al. (2008) from the radar images of Cassini. They determined an obliquity equal to 0.3 degree that is larger than the degree obtained for a rigid Titan in a Cassini state (e.g. Bills and Nimmo (2008)). Then Bills and Nimmo (2008) and Baland et al. (2011) suggested that this deviation is related to the influence of an internal ocean on the obliquity. Additional arguments in favor of this internal ocean have been obtained by electrical (Béghin et al., 2012), topographical (Nimmo and Bills, 2010) and gravitational (Iess et al., 2012) analyses of Titan from CassiniHuygens data.
Titan’s rotational motion has been measured by Stiles et al. (2008, 2010) that used the Cassini spacecraft’s radar images in order to follow landmarks at the surface. They obtained a nearsynchronous rotation for Titan with a drift rate of degree per day, i.e. degree per year (Stiles et al., 2010). This approach has been recently revisited by Meriggiola and Iess (2012) that determined a synchronous spinorbit motion within a residual of about 0.02 degree per year. The main advantage of their approach is the introduction of Titan’s librational motion in the reduction process. The librations describe the oscillations around the uniform rotational motion. Here we focus on the librations in longitude that correspond to the oscillations of the body principal axis projected onto the equatorial plane of the satellite. In the case of Titan, they have two distinct origins. The first one comes from the gravitational torque exerted by Saturn on the dynamical figure of Titan. The second one results from the coupling between the surface and the dense atmosphere of Titan, the atmospheric torque. Their amplitudes are modified by the interior parameters of Titan. However, Goldreich and Mitchell (2010) suggested that the elastic behavior of the icy shell is not negligible and may strongly reduce the librational motion. Such prediction has been confirmed by Van Hoolst et al. (2013) that computed the libration in longitude at the orbital frequency of Titan. In that case, Titan’s surface will deform instead of rotate since the ocean figure should always point toward the planet and exert an elastic torque on the icy shell.
The objective of this paper is to determine the librational reponse of an elastic Titan at various forcing frequencies resulting from its orbital motion. By including the atmospheric coupling, we want to decorrelate surface forcing from the internal geophysical properties related to the internal ocean and perturbation of Titan’s orbit. We investigate the wide spectrum of Titan’s librations on contrary to Van Hoolst et al. (2009, 2013) that focused on the diurnal frequency (Titan’s orbital frequency) and its harmonics, and on the Saturnian semiannual frequency (twice Saturn’s orbital frequency) for the atmospheric coupling. The main interest is that in the orbital motion, there are librations at the Saturnian semiannual frequency that comes from the interaction of Saturn with the Sun. Such librations have the same frequency than the main component of the atmospheric torque. In addition, the orbital frequency spectrum can be in, or close to, resonance with some proper frequencies as shown for the Galilean satellites (Rambaux et al., 2011). Finally, the periodic variation of the gravitational torque amplitude at orbital frequency provides longperiod terms in the libration with amplitudes dependent on the interior model.
In this paper, the recent interior models of Titan are also used (e.g. Fortes (2012); CastilloRogez and Lunine (2010); McKinnon and Bland (2011)) in order to compute the values of the proper frequencies and to discuss the differences in the amplitude of librations. The elasticity is investigated by computing the radial deformations of surfaces due to the tides and responsible for the gravitational torques amplitude variations.
Finally, we use the 3D atmospheric model from Lebonnois et al. (2012) that predicts an atmospheric torque smaller than in the Tokano and Neubauer (2005) paper used in all previous studies (Lorenz et al., 2008; Karatekin et al., 2008; Van Hoolst et al., 2009, 2013; Goldreich and Mitchell, 2010).
In the first part of the paper, the internal structure models selected for the librational computation are described. The properties of these models depend on the history and energy sources of Titan (e.g. Tobie et al. (2012)). We selected a broad range of possible internal scenarios for Titan in order to characterize the impact of the geophysical parameters on the librations. The atmospheric torques of Charnay and Lebonnois (2012) (called here CH12) and of Tokano and Neubauer (2005) (called here TO05) are described and discussed in Sect. 2. The orbit of Titan is then analyzed in Sect. 3 by using the frequency analysis method providing the spectrum of the orbital motion. In Sect. 4, the elasticity is introduced in the librational equations through periodic torques amplitudes and the librations are analytically determined for the rigid and elastic cases. Finally, the behavior of the libration angle at different forcing frequency ranges is analyzed and the influence of the geophysical parameters of the different interior models is discussed in Sect. 5.
1 Interior models
Titan’s internal structure has been revealed by accurately tracking the trajectory of Cassini spacecraft approaching Titan during six flybys (Iess et al., 2010, 2012). They measured the gravity field and its variations allowing to infer information on the density profile of the satellite. Since the inversion between the gravity field and density profile is not unique, only a range of models can be determined. The addon assumption that Titan is close to the hydrostatic equilibrium led Iess et al. (2010) to obtain an estimation of the moment of inertia (MoI) of Titan between and (where and are the mass and mean radius of Titan, respectively). Such small MoI implies an increase of density towards Titan’s center and requires a lowdensity core to match the mean density of kg m (Iess et al., 2010).
In parallel to these gravity measurements, two categories of thermal and chemical models have been developed to determine Titan’s internal structure (e.g. CastilloRogez and Lunine (2010, 2012); McKinnon and Bland (2011); Fortes (2012) and see the review of Tobie et al. (2012)). In the first category, as developed by Fortes (2012), the models described essentially the upper layers (dense and light ocean, comparison of solid models with pure water ice or methane clathrate layers), while in the second category the models focused on the inner core composition assuming the presence of a global ocean layer (e.g. CastilloRogez and Lunine (2012); McKinnon and Bland (2011)).
CastilloRogez and Lunine (2012) built interior models made of anhydrous silicate core surrounded by hydrated rock and water ice while McKinnon and Bland (2011) focused on hydrated silicates surrounded by mixture of rock and ice. These models of inner core allow the presence of a small iron core, and are compatible with an internal ocean as deduced by Iess et al. (2012).
Six different models reported in Table 1 have been selected, coming from Fortes (2012) (the models called F1, F2 and the solid model F3), CastilloRogez and Lunine (2010) (models CA10 and FE10) and McKinnon and Bland (2011) (model MC11). The FE10 model is similar to the CA10 model with an additional small inner iron core as suggested by CastilloRogez and Lunine (2010). The icy shell thickness has been taken equal to kilometers for each model. Such value corresponds to the upper bound obtained by topographical model of (Nimmo and Bills, 2010). The lower bound has been obtained by Béghin et al. (2012) by using the Schumann’s resonance in Titan’s atmosphere. The icy shell thickness is an essential parameter for Titan’s models with a rigid shell (Van Hoolst et al., 2009) but the elasticity strongly diminish its influence (see Section 5). Fortes (2012) used oceans with a bottom mean radius of kilometers (models F1 and F2). To be able to compare the influence of the ocean density and the inner core structure on the libration, the same ocean bottom mean radius is used for the CA10, MC11 and FE10 models. For a given set of solid layers sizes and densities, the ocean density is then adjusted to conserve Titan’s MoI and mass.
F1  F2  F3  
Layer  (kg m)  R (km)  (kg m)  R (km)  (kg m)  R (km) 
Icy shell  930.9  2575.5  930.9  2575.5  932.8  2575.5 
Ocean  1023.5  2475  1281.3  2475     
HP ice 1          1193.3  2429 
HP ice mantle  1272.7  2225  1350.9  2225  1268.7  2304 
Silicate mantle  1338.9  2163      1338.9  2185 
Silicate core  2542.3  2116  2650.4  1984  2584.1  2055 
(MR)  0.3414  0.3413  0.3413  
(MR)  0.0357  0.0357    
(MR)  0.2320  0.2133    

CA10  MC11  FE10  
Layer  (kg m)  R (km)  (kg m)  R (km)  (kg m)  R (km) 
Icy shell  932.8  2575.5  932.8  2575.5  932.8  2575.5 
Ocean  1317.6  2475  1315.5  2475  1224.1  2475 
HP ice mantle  1350.9  2225  1350.9  2225  1350.9  2225 
Rocky mantle  2520  2000  2200  2075  2520  2000 
Silicate core  3400  900  3300  1300  3400  900 
Iron core          6500  500 
(MR)  0.3402  0.3378  0.3336  
(MR)  0.0358  0.0358  0.0358  
(MR)  0.2095  0.2073  0.2096 
2 Atmospheric torque
Titan has a thick atmosphere extending up to 800 km with a surface pressure of nearly 1.5 bars. Exchanges of angular momentum happen between the atmosphere and the surface, producing an atmospheric torque influencing Titan’s rotation. The atmosphere dynamic is mainly driven by insolation, which creates a circulation zone called a Hadley cell between hotter and cooler regions. The circulation in Titan’s troposphere is essentially dominated by one Hadley cell extending from one pole to the other. At the equinox, the circulation reverses with the formation of two Hadley cells, and the ascending air zone (also called the Intertropical Convergence Zone) moves from one pole to the other. This reversal produces a seasonal angular momentum exchange with Saturn’s semiannual period ( terrestrial days) and has been suggested as the main torque influencing Titan’s rotation (Tokano and Neubauer, 2005). The total atmospheric angular momentum (AAM) is given by (Tokano and Neubauer, 2005)
(1) 
where km is Titan’s mean radius, ms is Titan’s surface gravity, is the surface pressure, is the zonal wind speed (i.e. the horizontal wind speed in the westeast direction and relative to the surface, is positive for eastward wind), is the latitude and is the pressure.
The atmospheric torque is then
(2) 
The zonal wind speed of an air parcel in the ascending zone is approximately null due to a strong surface friction. Its angular momentum is conserved during its meridional transport from the ascending zone to the poles. Thus, the AAM is maximum at equinoxes (i.e. when the ascending zone is the farthest from Titan’s rotation axis), and minimum at solstices (i.e. when the ascending zone is the nearest to Titan’s rotation axis) (Tokano and Neubauer, 2005). This change corresponds to an exchange of angular momentum with the surface by friction of winds.
The angular momentum exchange, and so the atmospheric torque, depends essentially on the extension of the Hadley cell in latitude and altitude (Mitchell, 2009). All previous studies on Titan’s nonsynchronous rotation (Lorenz et al., 2008; Karatekin et al., 2008; Van Hoolst et al., 2009, 2013; Goldreich and Mitchell, 2010) used the atmospheric torque at Saturnian semiannual frequency from the general circulation model (GCM) of Tokano and Neubauer (2005) ( kg m s; this torque is called here TO05). In this study, the atmospheric torque from the Titan IPSL GCM (Institut PierreSimon Laplace General Circulation Model) (Lebonnois et al., 2012) is also used. This model successfully reproduces the winds and the thermal structure observed by Cassini and Huygens (Lebonnois et al., 2012; Charnay and Lebonnois, 2012). Charnay and Lebonnois (2012) show that the Hadley cell is essentially trapped in the first two kilometers. This trapping, in addition to a lower latitudinal extension of the cell, reduces the angular momentum exchange by approximately a factor 10 compared to Tokano and Neubauer (2005).
The latitudinal extension of the Hadley cell depends a great deal on the thermal inertia of the ground (i.e. the ground thermal response to solar heating). Its value for Titan is unknown, it has been estimated to be around J m s K (Tokano, 2005) and used in Tokano and Neubauer (2005), but the value may be higher. Atmospheric torques have been calculated with a thermal inertia of and J m s K, the last one corresponding to a realistic maximal value. In that case, the Saturnian semiannual frequency component of the Titan IPSL GCM torque (called here CH12) is of about kg m s and is reduced to kg m s for a model with a thermal inertia of J m s K. In this study, we focus on the first model providing the largest Saturnian semiannual frequency torque amplitude.
The atmospheric torque CH12 with a thermal inertia of J m s K and the Saturnian semiannual frequency torque TO05 are represented in Fig. 1. Many other frequencies than the Saturnian semiannual one appear in the spectrum of the atmospheric torque CH12. The most important ones are frequencies corresponding to (linked to the eccentricity of the orbit of Saturn around the Sun), and (other harmonics from the seasonal cycle) Saturnian year with amplitudes below kg m s. The other frequencies correspond to atmospheric waves that are either free (baroclinic waves, kelvin waves,…) or forced by the thermal tides and the gravitational tides caused by Saturn (Tokano and Neubauer, 2002). The free atmospheric waves are very dependent on the model, but their amplitudes appear to be small. Concerning the thermal and gravitational tides, they have a very limited impact on the tropospheric winds and therefore on the exchange of angular momentum with the surface in the model from Lebonnois et al. (2012). Other effects may influence the atmospheric torque. Mitchell (2009) suggested the methane cycle could reduce the amplitude of the torque and produce a small drift while Tokano (2012) suggested the presence of mountains could increase the AAM a little without changing the torque very much. These effects remain poorly constrained but should only have a small impact.
3 Orbital forcing
The librational motion is the rotational response of Titan to the gravitational or atmospheric torques exerted on its dynamical figure. The gravitational torque depends on the relative distance between Titan and Saturn and on the true longitude of Titan (e.g. Rambaux et al. (2011)). The orbital elements are timedependent and can be approximated by quasiperiodic Fourier series resulting from the interaction of Titan with Saturn, the Sun or Iapetus (Vienne and Duriez, 1995).
Here, we develop a linear theory of the librational motion of Titan that is approximated in small quantities, the amplitude of libration and the eccentricity. The orbital motion of Titan comes from the JPL ephemerides that contain the last observation from Cassini space mission (Giorgini et al., 1996). The ephemerides are given over the period of 07Jan1800 to 07Jan2200 corresponding to an interval of 400 years. The frequency analysis method developed by Laskar (1988, 2003) and implemented in the TRIP software (Gastineau and Laskar, 2012) is used in order to decompose the true longitude into Fourier series.
The decomposition of the true longitude is given in Table 2. The notation of Vienne and Duriez (1995) is used, where represents the mean longitude of Titan ( with the longitude of node, the mean anomaly and the argument of pericenter, see Figure 2), and the longitudes of pericenter () of Titan and Iapetus, respectively, and the longitudes of node of Titan and Iapetus, and the mean longitude of the Sun. The true longitude is given by with the true anomaly.
The largest term of the Table 2 corresponds to the first term of the development of the difference between the true anomaly and the mean anomaly as function of the eccentricity , i.e. . The magnitude is equal to 11899.3237” corresponding to an eccentricity of 0.0288 and the frequency is equal to the orbital frequency. The following term corresponds to the second term of the development in the true anomaly, i.e. and oscillates at twice the orbital frequency. The other frequencies result from the motion of Saturn around the Sun (, , and here called Saturnian annual, semiannual and terannual terms respectively, and the combination ) and interaction with Iapetus () (see Vienne and Duriez (1995)). The last term has a frequency close to the combination of () where is the mean longitude of Jupiter. The last two terms of the Table 2 with periods of and days respectively are small terms of the true longitude series with magnitudes below . These terms have small magnitudes but due to their frequency values, the librational response amplitudes are greater than the 10 meters limit which is used for the Table 4 (see Section 5).
Let us note that the frequency analysis of the orbit over years is not sufficient to identify the periods longer than several hundred years. Vienne and Duriez (1995) have identified longer periods in Titan’s orbital motion listed in the analytical ephemerides. Frequencies associated with ( years), ( years), ( years) and ( years) possess high magnitudes. Since they are the signature of time variations much larger than the 400 years of the ephemerides, these terms act like secular deviation on the true longitude. These secular terms are removed from the analysis in the determination of the linear deviation of the true longitude. However, for a longer time range, the long period terms from longitudes of nodes and pericenters have to be considered explicitly.
Freq.  Period  Magnitude  Phase  Identification 

(rad/days)  (days)  (”)  (degree)  
0.394018  15.9464  11899.3237  163.3693  
0.788036  7.9732  212.5868  32.7941  
0.394081  15.9439  56.6941  68.1211  
0.001169  5376.6331  43.7313  66.0428  
0.000584  10750.3648  37.5508  138.4821  
0.392897  15.9919  31.5673  10.8789  
0.001753  3583.9304  5.6147  250.1412  
0.009810  640.4892  1.4983  77.2905   
4 Librational model
4.1 Rigid case
The librational equations describing the variations of Titan’s rotation are obtained from the angular momentum equation applied to each layer of the satellite. In the case of rigid solid layers, the hydrostatic equilibrium shape is composed of a static tidal bulge on which the planet exerts a gravitational forcing at each instant. In addition, the misalignment of the static tidal bulges of the shell and the inner core generates an internal gravitational torque coupling the upper layers to the inner core. Following the developments of Van Hoolst et al. (2008, 2009) and Rambaux et al. (2011), the equations governing the librational motions of the rigid triaxial shell and inner core perturbed by external and internal gravitational forces are written as
(3) 
where subscripts and refer to the icy shell and to the inner core respectively, is the libration angle, the mean anomaly, is the true longitude, is the rotation angle of the satellite’s longest axis measured from the line of the ascending node and its initial value for the layer (where subscript refers to the shell or the inner core), is the mean motion. and are the amplitudes of the effective gravitational torques exerted by Saturn on the dynamical figure of each layer, while is the amplitude of the internal gravitational torque exerted by the shell on the inner core due to their misalignment. The gravitational torques amplitudes are given by and , where , and are the principal moments of inertia (defined as ) and , and are the ocean pressure effect on the solid layer expressed as increment of inertia. The internal gravitational torque amplitude also depends on the geophysical parameters of the body as described by Van Hoolst et al. (2009). In this study, the dynamical equations are linearized as function of the eccentricity and libration amplitude. The frequencies of the forcing torques acting on the two solid layers are given by the Fourier series of developed in the previous section.
The proper frequencies (or natural frequencies) of the satellite are described in terms of free frequencies of the harmonic oscillator as done by Dumberry (2011). First, the free frequencies corresponding to the system (3) free of internal gravitational torque () are
(4)  
(5) 
Then, the free mode of oscillation corresponding to the system (3) without Saturn’s gravitational torque () is expressed as
(6) 
Finally, the proper frequencies of the system (3) forced by Saturn and internal gravitational torques can be written as
(7)  
(8) 
with . The values of and are given for each selected model in Table 3.
Using these expressions of the proper frequencies, the solutions of the system (3) are written as by setting , where and are respectively the amplitude of libration and the magnitude of perturbation, is the frequency of perturbation and the phase. The icy shell libration amplitudes for an oceanic model are then given by
(9) 
The libration amplitudes for a rigid solid model are given in Appendix.
As seen in Sect. 2, the atmospheric torque exerted on Titan’s surface is acting like a forcing term with a main component at Saturnian semiannual frequency and many components with amplitudes below kg m s. The atmospheric torque is introduced into the righthand side of the first equation of the system (3) as , where is the atmospheric torque magnitude, the frequency and the phase difference with . The libration amplitude is then decomposed into a sine and a cosine term expressed as
(10)  
(11) 
where and are the amplitudes of the sine and cosine terms respectively of the libration under orbital and atmospheric influence, as previously is the perturbing frequency, its phase, and the phase difference with the orbital perturbation. The cosine term implies a difference of phase with the libration amplitude due to orbital perturbation. For frequencies of the atmospheric torque which are not present in the orbital motion, the librational motion is obtained by doing in Eq. (10).
4.2 Elastic case
The elasticity is introduced by modeling the radial deformations of the surfaces of the satellite layers and the resulting variations of the moments of inertia are computed. Each surface composing the satellite is parametrized by where is the mean radius of the equivalent sphere and the angles , are the colatitude and longitude, respectively. Then is the radial deformation induced by the centrifugal and tidal potential.
Hinderer et al. (1982) used a decomposition of into spherical harmonics. At the second order, we have
(12) 
where and quantify the radial deformation normalized to and are associated Legendre polynomials. By using the definition of , the nonzero components of the inertia tensor can be written at first order in terms of , \linenomath
(13)  
(14)  
(15)  
(16) 
Love (1911) defined at the surface the Love number which characterizes the body response to second order perturbing potential. Here, we use the following definition for the radial deformation
(17) 
where represents the surface gravity of Titan and is the perturbing potential developed at order . is the dimensionless radial function corresponding to at the surface. By using eqs. (12) and (17), a straight relation exists between the , and the response to the perturbing potential proportional to .
As shown by Giampieri (2004), the development at first order in eccentricity of the perturbing potential (centrifugal and tidal) can be decomposed into a secular and a periodic part \linenomath
(18) 
where and are the mass and mean radius of the satellite, and with the mass of the planet and the semimajor axis of the satellite. The centrifugal potential has been supposed constant with a periodic variation of the potential only due to tides.
The secular part of the potential corresponds to static tides and is responsible for the static bulge of the satellite which is already included in our rigid case (see Sect. 4.1). The periodic part, which is of order of the eccentricity , governs timedependent deformation of the surfaces.
The deformations of the satellite layers due to the tides induce a timedependent gravitational potential. Mass distribution varies periodically and the inertia tensor of the body can be decomposed into a static part (corresponding to the static potential) and a periodic part associated to deformations around the static bulge. Using the definitions of and of the potential , the periodic components of the inertia tensor are given by
(19)  
(20)  
(21)  
(22)  
(23) 
where the following definitions are used:
(24) 
and
(25) 
The total deformation factors are given by a combination of static deformations () and periodic deformations ().
The periodic components of the inertia tensor and are associated with and variations called radial tides and is associated with variations called zonal tides. The component corresponds to a periodic bulge shifted by in longitude from the principal axis of inertia (called librational tide)(Murray and Dermott, 1999; Greenberg et al., 1998).
The external gravitational torque exerted on a layer is given by where is the external gravitational potential of the deformed satellite layer given by (Jeffreys, 1976) \linenomath
(26) 
where and are respectively the top and bottom radius of the layer . The zcomponent of this torque can be written at second order in eccentricity \linenomath
(27) 
Here describes the torque amplitude exerted by the planet on the static bulge of the layer as in the rigid case (see Figure 3(a)). is the torque amplitude exerted on the radial tidal bulge of the layer with and is the contribution of the ocean pressure on the periodic figure based on equation (24).
The last two terms of equation (27) are the torque exerted by the planet on the librational bulge. It can be notice that the development in second order of for a Keplerian orbit cancels out , and in this case the radial tides terms vanishes. The signs of the remaining librational and static terms at first order in are in phase and opposed, which means that the torques are counteracting as illustrated by the gravitational forces on Fig. 3(b).
The total internal gravitational torque exerted by outer layers on the interior is given by where is the potential exerted by outer layer on an inner point given by (Jeffreys, 1976) \linenomath
(28) 
where refers to the inner core of mean radius . The zcomponent of the internal gravitational torque exerted on the inner core is then \linenomath
(29) 
where and are defined such as \linenomath
(30)  
(31) 
with the deformation factor of the static bulge. In the internal gravitational torque, three main components are identified: the interaction between the static bulges of the shell and the inner core through like in the rigid case (see Figure 3(a)), the interaction between the periodic bulge of the shell and the static bulge of the core (subscript , Figure 3(c)), and the interaction between the static bulge of the shell and the periodic bulge of the core (subscript , Figure 3(d)). Here, the attraction between periodic bulges, which is of higher order in eccentricity, is neglected. Terms in are due to radial tides while terms in are due to librational tides. All the term but one have the same sign in the internal torque, which means that, for small and positive values of , they are acting in the same direction and tend to align the shell and the inner core figures. The last term is counteracting however, it corresponds to the attraction of the static bulge of the core by the librational bulge of the shell (Fig. 3(c)).
The equations of librations are then linearized and written at second order in eccentricity:
(32) 
where the second term of the lefthand side is due to zonal tides on the layer with . The internal and external gravitational torques at first order in eccentricity are dominated by the static and librational tides and the second order is dominated by the radial tides.
4.3 Analytical resolution of the elastic case
The system (32) is composed of two differential equations of second order in and . To solve this system, we transform it into a first order equation system in by setting , , , and . The system is then decomposed into a static and a periodic part such as where z is the vector of components , is the vector of forcing terms, and are the matrices of the static and periodic coefficients, respectively. The static coefficient matrix is diagonalized by defining a vector such that where corresponds to the eigenvector matrix of . The system is then reduced to
(33) 
where is the diagonal matrix, is the transformation of and and are the transformation of decomposed into first order and higher order terms, respectively.
To solve this new system, the perturbation method used by e.g. Robutel et al. (2011) is followed and is decomposed in decreasing amplitude terms such that where subscripts and refer to first and second order solutions. We then have to solve \linenomath
(34)  
(35)  
The system is first solved and the solutions are substituted in the system to obtain the second order solutions. Then the solutions are transformed back to our libration angles also decomposed in and . The solution is not computed here since we are only interested in solutions of order 2 in eccentricity.
Since the deformation is also periodic with frequency , the solutions of the elastic case are given by
(36) 
where the third subscript corresponds to the frequency. At the orbital frequency, the amplitude is and \linenomath
(37)  
(38) 
As suggested in the previous section, the torques variation and present in the solution at first order in eccentricity are only due to librational tides components.
The solution for an elastic solid body without ocean is given in Appendix.
5 Librational response and application to interior models
5.1 Discussion of the rigid solution
The libration amplitude obtained in the previous sections for the rigid and elastic cases are described and analyzed in this section. The libration responds to the perturbations according to three regimes depending on the values of the forcing frequencies with respect to the values of the proper frequencies: a high forcing frequency regime, a low forcing frequency regime, and a resonant regime where the forcing frequencies are close to the proper frequencies.
For all regimes, the solutions are simplified with the assumption that the icy shell thickness is small. is about 10 to 20 times larger than for the six models selected in Sect. 1. So the proper frequencies can be developed at first order in and such as
(39) 
(40) 
Using these developments, it can be noticed that for low forcing frequencies , the libration amplitude of the rigid case is simply , i.e. the librational response of the body or the icy shell will follow the magnitude of the perturbation. In that case, the librational response is dominated by the external gravitational torque and the internal structure has a negligible influence.
If the atmospheric torque excites the system then the amplitude of the librational response to low frequency forcing will be different from the magnitudes of the orbital perturbations. This difference would contain information on the atmospheric torque coupled to internal structure through the sine term
(41) 
At low forcing frequency, the librational response due to the atmosphere is then inversely proportional to the external gravitational torque amplitude on the shell and dependent on the internal gravitational coupling through the constant which depends on the inner core composition.
At high forcing frequencies, the rigid solution (eq. (9)) can be simplified to
(42) 
when is large in front of and . Here, the dynamics is dominated by if . By writing as function of the equatorial flattening of the layers (see Van Hoolst et al. (2009)), and by noticing that the flattening of the ocean is almost equal to the one of the shell from our interior models (the variations are only ), is written
(43) 
where the icy shell thickness and equatorial flattening have been introduced. The free frequency and the libration amplitude are then dependent on the ratio between ocean and icy shell densities, and inversely proportional to the icy shell thickness.
Finally, if the forcing frequencies are close to the proper frequencies of the system, the librational response is dominated by the resonant behavior. The amplitude of the libration is then significantly increased due to the presence of a small divisor. If the libration amplitude becomes too large, the linearized equations used are not enough to describe the dynamics. For Titan, the proper frequencies listed in Table 3 compared to the forcing frequencies listed in Table 2 show that there is no close resonance between these frequencies.
5.2 Discussion of the elastic solution
Here, we assess the effect of the elasticity on the libration solution which is dominant only for high forcing frequencies. The elasticity induces periodic variations of the torques amplitudes that are identified in the diurnal libration amplitude through two terms denoted by and (eq. 36) corresponding respectively to the variation of the shell and the inner core inertia. The form of the elastic case diurnal solution at first order in eccentricity (eq. (36)) is similar to the rigid case solution (eq. (9)) and it can also be simplified as
(44) 
with large in front of and , and is defined as
Thus, the libration amplitude at orbital frequency will be largely reduced compared to the rigid case if is large enough in front of , i.e. if the shell librational tide contribution is large enough in front of the external gravitational torque exerted on the static bulge. To evaluate this contribution, the reduction rates of the torques amplitudes are defined as for the shell and for the inner core. Their numerical values are given in the next section.
The second order librational reponse to the forcing frequency is identical to the rigid case and the behavior is described in the previous section. Terms of frequencies also appear in the second order solution due to the modulation of the gravitational torque amplitude. Since can be small compared to , the low frequency forcing can contribute to the orbital frequency libration but these terms have small amplitudes below one meter. For the special case where is close to , the terms of frequency contribute to the libration with their amplitude being inversely proportional to the proper frequencies (see eq. 45). The solution can be written for the frequency , and
(45) 
In that case, since the values of proper frequencies and are small (see Table 3), large values of magnitude can provide a nonnegligible signature in the libration as detailed in the next section.
(rad day)  (rad day)  (rad day)  
F1    
F2    
F3      
CA10    
MC11    
FE10   
5.3 Application
Now, the librational solutions are applied to the different interior models introduced in Sect. 1. The librations amplitudes projected onto Titan’s equator are given in Table 4. As discussed previously the librations can be classified into three regimes depending on the value of the forcing frequency with respect to the values of the proper frequencies.
For the low forcing frequencies, which are lower than the proper frequencies, the libration amplitude is almost equal to the magnitude of the orbital perturbation for any interior models. As seen in Table 4, the resulting libration amplitudes are large, around meters and meters for the Saturnian semiannual and annual frequencies. The small differences in the amplitudes at Saturnian semiannual and annual frequencies are due to a second order effect in the proper frequencies (eqs. (39) and (40)) and are related to the different inertia of the models, with a larger value of libration for the solid model induced by the small value of proper frequency .
Amplitudes (m)  
Freq.  Period  F1  F2  F3  Identification 
(rad/days)  (days)  
0.394018  15.946441  62.061  86.257  50.661  
0.001169  5376.6317  552.772  552.309  560.119  
0.000584  10750.4115  470.287  470.183  471.838  
0.001753  3584.9304  72.120  72.000  74.306  
0.000063  99027.4111  47.477  50.171  31.500  
0.001121  5606.2511  26.687  28.182  17.978  
0.009810  640.4878  11.868  16.118  24.275   
0.001169  5376.6317  517.761  515.967  523.095  Tokano 05 
0.001169  5376.6317  548.164  548.107  556.305  Charnay 12 
0.394018  15.9464  319.036  384.527  52.032  Rigid case 
Amplitudes (m)  
Freq.  Period  CA10  MC11  FE10  Identification 
(rad/days)  (days)  
0.394018  15.946441  72.751  76.679  75.420  
0.001169  5376.6317  552.265  552.284  552.478  
0.000584  10750.4115  470.171  470.174  470.217  
0.001753  3584.9304  71.994  72.002  72.054  
0.000063  99027.4111  51.963  51.996  50.714  
0.001121  5606.2511  29.185  29.205  28.492  
0.009810  640.4878  16.956  17.351  16.721   
0.001169  5376.6317  515.887  515.995  516.814  Tokano 05 
0.001169  5376.6317  548.084  548.055  548.028  Charnay 12 
0.394018  15.9464  391.131  387.614  358.397  Rigid case 
In the case of rigid layers, the librational response at orbital frequency are expected to be strongly dependent on the satellite inertia (Van Hoolst et al., 2009). Differences between oceanic and solid models would then appear clearly through the variation of the amplitude by at least a factor six, as for example the values of 391 meters for CA10 and 52 meters for F3. In the rigid case, equation (43) shows that the amplitude would mainly depend on the ocean and icy shell densities ratio and on the icy shell thickness.
As described in 5.1, the total torque acting on the solid layers is modified by reduction rates and when elasticity is taken into account. To evaluate these rates, the radial function is computed by the numerical integration of the gravitoelastic equations set of Alterman et al. (1959) with the SatStress numerical code developed by Wahr et al. (2009). For the different interior models selected, the corresponding surface Love numbers are contained between and for oceanic models, and about for the solid model. In an elastic and solid Titan, only a small deformation is necessary for the shear stresses to balance the tidal force and thus a smaller value of . There are no shear stresses in an inviscid fluid layer allowing the Love number value to increase (Rappaport et al., 2008). For the oceanic models, values obtained in the shell are at least forty times larger than at the surface of the core. It means that the periodic radial deformation of the core is negligible compared to the shell as well as its radial and librational bulges. Thus the internal torque amplitude can be ignored in the forcing terms.
For the CA10 model, we find