Librations of a synchronous satellite

Interpreting the librations of a synchronous satellite – How their phase assesses Mimas’ global ocean

Benoît Noyelles NAmur Center for CompleX SYStems (NAXYS) – University of Namur – Rempart de la Vierge 8 – B-5000 Namur – Belgium

Most of the main planetary satellites of our Solar System are expected to be in synchronous rotation, the departures from the strict synchronicity being a signature of the interior. Librations have been measured for the Moon, Phobos, and some satellites of Saturn. I here revisit the theory of the longitudinal librations in considering that part of the interior is not hydrostatic, i.e. has not been shaped by the rotational and tidal deformations, but is fossil. This consideration affects the rotational behavior.

For that, I derive the tensor of inertia of the satellite in splitting these two parts, before proposing an analytical solution that I validate with numerical simulations. I apply this new theory on Mimas and Epimetheus, for which librations have been measured from Cassini data.

I show that the large measured libration amplitude of these bodies can be explained by an excess of triaxiality that would not result from the hydrostatic theory. This theory cannot explain the phase shift which has been measured in the diurnal librations of Mimas. This speaks against a solid structure for Mimas, i.e. Mimas could have a global internal ocean.

Resonances, spin-orbit – Rotational dynamics – Satellites, shapes – Celestial mechanics – Saturn, satellites

1 Introduction

The Cassini space mission has given us invaluable data on the satellites of Saturn. In particular, we now dispose of measurements of their shapes as triaxial ellipsoids (Thomas, 2010, e.g.), and the rotation of Janus, Epimetheus (Tiscareno et al., 2009), Mimas (Tajeddine et al., 2014), Enceladus (Thomas et al., 2016), and Titan (Stiles et al., 2008; Meriggiola et al., 2016) have been measured. An issue is to get clues on the interior from these informations.

Most of the natural satellites of the giant planets are expected to have reached a state of synchronous rotation, known as Cassini State. This is a dynamical equilibrium from which small departures are signatures of the internal structure. These small departures are longitudinal and latitudinal librations, the latter ones translating into an obliquity. The main part of the longitudinal librations is a periodic diurnal oscillation, named physical librations. For a rigid body, they read (Murray & Dermott, 2000, Eq. 5.123, e.g.)


being the orbital eccentricity of the satellite, its orbital frequency, and the frequency of the small proper oscillations around the equilibrium. We have:


where the quantities are diagonal elements of the tensor of inertia of the satellite, and is an eccentricity function made popular by Kaula (1966):


This function is also known as the Hansen function , and is present in Cayley (1861).

The quantity , sometimes written as , represents the equatorial ellipticity, or the triaxiality, of the distribution of mass inside the satellite. This formula assumes a rigid shape, i.e. the tensor of inertia is constant.

A seminal paper by Goldreich & Mitchell (2010) has shown that if the satellite is not strictly rigid, but viscoelastic, then a restoring torque tends to counterbalance the tidal torque of the parent planet, to lower the amplitude of the physical librations. This has motivated recent studies (Van Hoolst et al., 2013; Richard et al., 2014; Makarov et al., 2016), which revisit the theory of the librations in including a tidal parameter, or , which characterizes the amplitude of variation of the gravity field, or of the surface, at the diurnal, or orbital, frequency . In these studies, the body is assumed to be at the hydrostatic equilibrium.

The theory of the hydrostatic equilibrium tells us that the mass distribution in the body, i.e. its inertia, is shaped by its rotation and the tidal torque of its parent planet, while the inertia rules its rotation. However, in most of the studies, the inertia is mainly composed of a constant component which has no chance to be shaped by the rotation, while being assumed to correspond to the hydrostatic equilibrium.

In this paper I go further, in assuming that part of the mass distribution of these bodies is frozen, while part of it is still being shaped by the rotational and tidal deformation. For that, I first show from the measured radii that the departures from the hydrostatic equilibrium are ubiquitous in the system of Saturn (Sect. 2). Then I express the tensor of inertia of a satellite orbiting a giant planet, in considering a frozen triaxiality superimposed with elastic deformation (Sect 3). I then deduce the librational dynamics of the satellite (Sect. 4), which I apply to the specific cases of Epimetheus and Mimas, for which longitudinal librations have been measured, and which are assumed to have rigid structure. Finally, I introduce the dissipative part of the tides (Sect. 7), to investigate their influence on the measurements. The reader can refer to Tab. 6 for the notations.

2 Departures from the hydrostatic equilibrium

Usually the shape of such a body is assumed to be the signature of an hydrostatic equilibrium. This means that the mass distribution is an equilibrium between the gravity of the body, and the rotational and tidal deformations it is subjected to. For a natural satellite orbiting a giant planet, it is assumed that the tidal deformation is only due to the planet, and that its spin rate is equal to its orbital rate, i.e. it rotates synchronously. This gives the following relations (Matsuyama & Nimmo, 2009; Correia & Rodriguez, 2013; Tricarico, 2014):


being respectively the 3 planet-facing, orbit-facing, and polar radii, and the classical Stokes coefficients, and . In this last formula, is the mass of the satellite, and its mean radius. The relation (4) holds for the shape, and (5) for the gravity field. The synchronous rotation results in a forcing of the triaxiality of the satellite, since it on average always presents the same face to its parent planet.

The Cassini mission has given us invaluable data on the shapes (Tab. 2) and gravity fields (Tab. 1) of some Saturnian satellites, and a comparison with the numbers predicted by the hydrostatic theory reveals some departure (Tab. 3), as for the Moon (Jeffreys, 1937; Lambeck & Pullan, 1980; Garrick-Bethell et al., 2014, e.g.) and Phobos (Le Maistre et al., 2013, e.g.).

These departures from the hydrostatic equilibrium may have different origins, such as a fossil shape, internal processes, or impacts…In the following, I consider them as frozen inertia.

3 The tensor of inertia

I consider a homogeneous, triaxial and synchronous satellite. Its orbital dynamics comes from the oblate two-body problem, i.e. the semimajor axis and the eccentricity are constant, and the mean longitude and longitude of the pericentre have a uniform precessional motion, the associated frequencies being and , respectively. Moreover, its orbit is assumed to lie in the equatorial plane of the planet, as a consequence the satellite has no obliquity. I also neglect the polar motion, the angular momentum of the satellite being collinear to its polar axis of inertia.

The tensor of inertia of the satellite can be decomposed as:



  • is the frozen component. It is constant and its physical origin is here not addressed,

  • is the inertia of a spherical body. We have


    which gives, for our homogeneous body:


    being the classical Kronecker symbol, the density at the distance , and a volume element.

  • is the deformation of the inertia induced by the tides, and

  • is the rotational deformation.

In this section, I neglect the dissipative part of the deformation. As a consequence, the sum represents an elastic deformation, and will be denoted , e standing for elastic. The dissipation will be considered in the Sect. 7.

3.1 The frozen inertia

The frozen part of the tensor of inertia is constant. As a symmetric tensor, it can be written under a diagonal form in an appropriate reference frame , whose axes are the principal axes of inertia, i.e.:


with . In the literature dealing with rigid rotation, these 3 moments of inertia are often written , and , respectively.

3.2 The elastic inertia

The elastic tensor results from the combined action of the tidal torque of the parent planet and the rotation of the satellite. Following (Williams et al., 2001, e.g.), we have




where , and are the coordinates of the unit vector pointing to the parent planet in the reference frame of the principal axes of inertia of the satellite. In this frame, is the rotation vector of the satellite. is the second order gravitational Love number, which characterizes the amplitude of deformation of the gravity field. is the mass of the parent planet, is the mean radius of the satellite, and is the distance between the barycentres of the planet and the satellite.

Our assumptions imply , and , this results in:


being the semimajor axis of the satellite.

The coordinates of the parent planet and can now be expressed with respect to the orbital elements of the satellite and to its orientation. Since it is assumed to have only a longitudinal motion, its orientation is given by an instantaneous rotation angle around such that and its origin corresponds to a state in which the long axis of the satellite points to the barycenter of the parent planet.

In the reference frame , we have


where stands for the true anomaly, and for the longitude of the pericenter. These quantities can be introduced in the formulae (12) to (15), before being expanded with respect to the eccentricity using the classical formulae (Duriez, 2002, Eq. 3.117), (Murray & Dermott, 2000, Eq. 2.84 & 2.85):


where are the Bessel functions of the first kind. I now introduce the argument of the synchronous spin-orbit resonance (Henrard, 2005, e.g.) and the mean anomaly to get, up to the second order in the eccentricity:


The argument of the spin-orbit resonance is usually set to , except in (Correia & Rodriguez, 2013, Eq. 37-38). Here, I keep as a variable, since it will be involved in the equation of the librations. Since the rotation is assumed to shape the inertia, and the inertia rules the rotation, then should be let free to vary. I will anyway assume that the frequency associated with is null, i.e. I assume to be constant on average, which is consistent with the resonant locking. The formulae (21) to (24) imply that the tidal response of the satellite does not depend on the frequency of the excitation. Actually the tidal parameter is frequency-dependent. To model this dependency I define the following frequencies:


and the elastic tensor of inertia becomes:


In the literature, is sometimes denoted as the fluid, or secular, Love number . It represents an indefinitely slow deformation. However, and are often assumed to be equal and denoted as . The periods of the deformations associated are respectively the orbital and half the orbital ones, which are also the diurnal and semi-diurnal periods of the satellite.

4 The librational dynamics

4.1 The librational equations

The gravitational torque of the parent planet, which acts on the satellite, is collinear to the polar axis, since we consider a planar orbit. Its non-null component reads (Williams et al., 2001, e.g.):


After expansion with respect to the orbital elements, I get, up to the second order in the eccentricity:


while is also the time-derivative of the norm of the angular momentum. This yields




The frozen component of , i.e. , should appear in this equation, but since is not affected by , then its secular part has the same behavior than its frozen part, distinguishing it is in that specific case useless. So, it is kind of absorbed in .

In equating the formula (33) with the (35), I get the differential equation ruling the longitudinal libration of the satellite, i.e.




4.2 The librational solution

The equation (36) can be simplified for an analytical resolution in neglecting the variations of with respect to the ones of , and . It becomes:




and this results in


after damping of the proper oscillations, their frequency being . Here the semi-diurnal oscillations have been dropped. The numerical application shows that their amplitude is negligible.

The Tab. 4 proposes a comparison between this study and two previous ones, by Van Hoolst et al. (2013) and Richard et al. (2014). These two studies aimed at modeling the librations of a satellite composed of a viscoelastic crust coating a global ocean, itself enshrouding a rigid core. Hence, they consider pressure and gravitational couplings between the different layers, that I do not have here. In considering only the outer shell in their studies, I managed to get equations like my Eq. (48), from (Van Hoolst et al., 2013, Eq. 45) and (Richard et al., 2014, App. A), with the help of (Richard, 2014).

In those two studies, , and stand for the mean values of the diagonal elements of the tensor of inertia , i.e. , and . Moreover, the triaxiality is supposed to be entirely due to the rotational and tidal deformations, i.e. there is no frozen component in , and the second degree in the eccentricity has been neglected. In (Van Hoolst et al., 2013), stands for , while in (Richard et al., 2014), the topographical Love number has been used.

My formula for is consistent with the one of Richard et al. (2014), while Van Hoolst et al. (2013) have a signature of the tidal deformation at the diurnal frequency. However, these two previous studies are consistent with each other for the coefficient for .

A difference between these studies and mine is that they assume the resonant argument to be fixed to in the expression of the tensor of inertia (Eq. 28 to 30), while it is a variable of the problem. If we set in the tensor of inertia, then the gravitational torque defined by the Eq. 32 becomes, after averaging over the mean anomaly :


i.e. only the frozen component remains on average. If we had only an elastic inertia, then the average torque would be


that comes from


So, a frozen component in is needed to get a non-null mean gravitational torque.

5 Modeling the tides

5.1 The Maxwell model

The classical Maxwell model (Karato, 2008, e.g.) gives a pretty good estimation of the frequency-dependency of . It depends on one parameter, the Maxwell time , being the viscosity and the rigidity. The complex Love number reads




and is the complex compliance defined as:


being the tidal frequency, and the surface gravity of the body. is the real part of , i.e.


5.2 An improvement at high frequencies

At high frequencies, the Maxwell model lacks accuracy because it does not render the fact that anelasticity dominates viscoelasticity, i.e. the response of the material is not instantaneous anymore. The Andrade model is therefore more appropriate. This is the reason why Efroimsky (2012a) proposed the so-called Andrade-Maxwell model, that corresponds to the Andrade model at high frequencies and to the Maxwell model at lower frequencies. Its complex compliance reads:


being the classical function defined as


We can see that this model depends on 3 tidal parameters, which are the Maxwell time , an Andrade time that has been introduced by Efroimsky (2012a, Eq. 78), and an Andrade parameter . The Andrade time should be equal to the Maxwell time to have a continuous transition between viscoelasticy and anelasticity, i.e. . is usually assumed to lie between 0.1 and 0.5, I will take the classical value . The resulting expression for the Andrade-Maxwell model is




The Andrade parameter is an order of magnitude of the period above which the excitation will generate the Andrade creep, responsible for anelasticity. Setting to the infinity renders the Maxwell rheology.

The Fig. 1 illustrates the elastic tides given by these two models, for Mimas and Epimetheus, in using the physical parameters given in Tab. 5.

6 Application to Epimetheus and Mimas

Diurnal librations have been measured for Epimetheus and Mimas, thanks to Cassini data (Tiscareno et al., 2009; Tajeddine et al., 2014). These two bodies are a priori assumed to be solid bodies, which legitimates the use of this model to explain their librations.

6.1 Methodology

A numerical integration of the Eq.(36) is performed, with the 10th order Adams-Bashforth-Moulton predictor-corrector scheme. A frequency analysis is then made to decompose the variable under a quasi-periodic formulation. This way, the proper oscillations and the forced ones are clearly identified. The frequency analysis algorithm is based on Laskar’s original idea (Laskar, 1999, 2005), has been adapted for the rotational dynamics by Noyelles et al. (2008), and used many times since. The initial conditions of the numerical integration are given by the analytical solution at , i.e. the Eq.(51) and its time-derivative, while the parameters are gathered in Tab. 5.

For each of these satellites, several simulations are run, which differ by the fraction of non-hydrostaticity. This starts from the formula giving the mean value of :


and then, two approaches are considered:

  1. either the secular Love number is set to be constant. I first set it to before discussing the implications of a smaller , which could be more physically relevant. I define the parameter such that


    This approach suggests that the body is composed of the superimposition of 2 mass distributions, an elastic and a frozen ones. The elastic contribution is fully consistent with the hydrostatic theory, and the frozen contribution supplements the elastic one, as an excess of triaxiality. This excess represents a fraction of the elastic contribution,

  2. or the mean value of is considered to be constant and determined from the shape of the satellite in assuming a constant density, i.e.


    which implies


    I define the parameter such that


    and is obtained from the formula (67). This would mean that the mean inertia of the body is consistent with the hydrostatic theory, but part of this inertia is frozen, the remaining part being still elastic.

6.2 Epimetheus

For Epimetheus, a physical libration of has been measured by Tiscareno et al. (2009). Such a large number is due to the triaxiality of the satellite, which makes the frequency of the proper oscillations close to the orbital one, resulting in a large amplitude of response (Noyelles, 2010, Fig. 6).

The simulations with an excess of triaxiality result in a larger (Fig. 2) than anticipated by the measured shape, but permit to reach the measured amplitude for .

for the elastic part should be seen as an incomplete model, in the sense that no rigidity is considered. The Fig. 3 displays for smaller values of . In particular, we have for .

These numbers should be balanced by the uncertainties of a few kilometers on the radii and , and some departures from an actual triaxial shape. We can also see that the analytical formulae for (Eq. 49) and the amplitude 1 (Eq. 51) are validated by the numerical simulations.

The simulations assuming a constant triaxiality for Epimetheus (Fig. 4) still validate the analytical formulae for and the diurnal amplitude. Moreover, they do not result in a large enough amplitude to match the observed libration.

This study assumes that Epimetheus is triaxial. Actually the visual aspect of Epimetheus suggests significant departure from the ellipsoid. To the best of my knowledge, no higher order shape model has been published, this is why I assume the triaxial ellipsoid as a good enough model. Of course, the departures to the triaxiality should contribute to the librations.

6.3 Mimas

Tajeddine et al. (2014) have measured from Cassini data a surprisingly high diurnal amplitude for Mimas of arcmin, while a theoretical study by Noyelles et al. (2011) predicted half this number. That study assumed Mimas to be a rigid body made of a mixture of ice and silicates, with a radial gradient of porosity. It was modeled as two homogeneous layers, the inner one being denser than the outer one. Tajeddine et al. (2014) tested 5 plausible interiors to explain their measurements, and kept only two: Mimas had either a global ocean beneath an icy crust, or a highly elongated core of pretty heavy elements. This last explanation would be consistent with a model of formation of the Saturnian satellites proposed by Charnoz et al. (2011). In that scenario, the satellites would have been formed as droplets from the rings before migrating outward to their present location. The highly triaxial shape of the core would be the signature of a former more rapid rotation, itself meaning that the shape froze when Mimas was closer to Saturn than it is now. It translates into an excess of triaxiality.

The two figures 5 and 7 still validate the analytical formulae for and the diurnal amplitude. Moreover, they permit the measured amplitude for and , while a constant prohibits it. would imply (Fig. 6).

7 Influence of the creep

I now introduce the imaginary part of the Love number, that is often denoted as , where is a dissipation function. This imaginary part introduces a dissipation as a lag between the gravitational excitation of the parent planet and the response of the satellite. This lag depends on the tidal frequency, this is why this effect will appear as .

7.1 Equations of the problem

As shown for instance in (Williams & Boggs, 2015), the dissipative tide can be introduced in replacing in the formulae (28) to (31) by and by . This geometrically means that the parent planet raises two bulges at the surface of the satellite, i.e. a bulge which is aligned with the direction of the planet, and which is responsible for the elastic tide, and a bulge in quadrature, which is responsible for the dissipative tide (Zahn, 1966; Efroimsky, 2012b, e.g.). The elastic inertia is replaced by a viscoelastic one . This yields:


If we consider the Maxwell rheology, then the quantity is obtained from the formula (58), i.e.


It is in particular expected from any rheology that is null at the zero frequency, i.e. no dissipation occurs for a constant excitation.

The gravitational torque (Eq. 33) then becomes