© ESO, 2011

1. Introduction

The Cassini spacecraft gives us the unique opportunity to obtain an accurate set of geodetic data for the icy satellites of Saturn such as the shape, the gravitational field, and the rotational data (Thomas 2010). The flybys of Mimas have provided high-resolution images of the surface in the finest detail yet seen (Roatsch et al. 2009). The Cassini spacecraft has detected temperature inhomogeneities (Howett et al. 2011), which are usually attributed to exogenic processes. The theoretical model of Mimas’ rotational state can be used to interpret the Cassini data and to better understand its interior and evolution.

As our Moon, Mimas is in synchronous rotation with its planet and shows almost the same face towards Saturn. Moreover, it is considered to have a large librational amplitude (Comstock & Bills 2003). The rotational state of a synchronous body depends on the mass distribution of the body, and therefore it is a signature of its internal structure. Here, we propose to model the rotation of Mimas considering it as a rigid body. A highly rigid interior of Mimas for most of its history is consistent with its un-relaxed shape (Thomas et al. 2007; Thomas 2010).

Because the distant spacecraft flybys of Mimas do not allow the determination of the gravity field nor the gravity harmonics, the mass of Mimas is determined from an analysis of its orbital resonances with Tethys and Methone (Jacobson et al. 2006). But its internal structure remains uncertain. The mean density of 1.15 g cm^-3 suggests that Mimas consists of a homogenous mixture of ice and rocks. The observed shape of Mimas by Voyager has been interpreted as an indication of interior mass concentration, which can be either caused by internal differentiation (Dermott & Thomas 1988) or radially variable porosity (Eluszkiewicz 1990). However, Cassini observations showed that Mimas’ shape, although a triaxial ellipsoid, deviates slightly from a hydrostatic shape and therefore an interpretation of the interior configuration from the shape is limited (Thomas et al. 2007; Thomas 2010). In the present study, we consider Mimas to be composed of two rigid layers. We consider a hydrostatic and a non-hydrostatic interior model. The interior models considering compaction of ice-silica particle mixtures (Yasui & Arakawa 2009) are expected to yield realistic principal moment of inertia A < B < C (Eluszkiewicz 1990). Because Mimas orbits close to its parent planet, the present-day diurnal tidal stresses can be significant, and we took the tidal effects into account as well.

The paper is structured as follows: we first model the internal structure of Mimas considering two different assumptions: that Mimas is in hydrostatic equilibrium, and that its ellipsoid of gravity is proportional to its ellipsoid of shape. From these two assumptions we derive 23 models of Mimas. Then we perform numerical integrations of the rotation of these “Mimases” in a full three-degree of freedom conservative models. Finally, we check the influence of the tides on the equilibrium orientation of Mimas’ long axis.

2. Internal structure

Interior structure models of planets and natural satellites are in general non-unique because one has fewer constraints than unknowns. For Mimas we have only two constraints: the mean radius R and the mean density ρ, or the mass m (Table 1). We are here concerned with the rotation, and accordingly the moment of inertia differences are the main point of interest.

Because we have only few constraints, we preferred a most simple interior structure for our models. We assumed a two-layer interior structure model with a rocky core and icy mantle. The models with a variable porosity, which is supported by recent compaction experiments (Yasui & Arakawa 2009), yield realistic moments of inertia (Eluszkiewicz 1990). The mean density of Mimas is close to the density of water ice and the interior is consistent with an icy mantle and a small rocky core; alternatively Mimas can be homogenous with a variable or constant porosity.

With Mimas’ figure departed from hydrostatic shape (Thomas et al. 2007; Thomas 2010) and without the knowledge of C₂₂ and J₂, we do not know the moment of inertia differences. As for the Moon (see e.g. Lambeck & Pullan 1980), the internal mass distribution could be a fossil shape that dates back from an earlier orbital position where the tidal heating was significant with higher orbital eccentricity and obliquity or a frozen shape following a large impact. Mimas shows a heavily cratered surface without signs of geological activity for billions of years. Its high free eccentricity gives another reason for low internal activity. This anomalously high eccentricity of Mimas can be explained by a passage through several past resonances (Meyer & Wisdom 2008). A determination of the gravity coefficients C₂₂ and J₂ =  −C₂₀ is necessary to conclude on the hydrostatic equilibrium because the figure of Mimas may not represent the real flattening or the internal mass determination, as is the case for Titan (Zebker 2009). As suggested by Johnson et al. (2006), large impact craters and heating of Mimas’ hemispheres by Saturn to different amounts may be potential sources of the considerable non-hydrostatic anomalies that could impede an accurate interpretation of the shape data.

2.1. Hydrostatic approximation

Because we do not know the gravity coefficients C₂₂ and J₂, we used a simple approach. For a two-layer interior model the core radius R_c can be determined if the densities of the rocky core ρ_c and the icy mantle ρ_s are known, i.e. ${\mathit{R}}_{\mathrm{c}}\mathrm{=}\mathit{R}{\left(\frac{\mathit{\rho }\mathrm{-}{\mathit{\rho }}_{\mathrm{s}}}{{\mathit{\rho }}_{\mathrm{c}}\mathrm{-}{\mathit{\rho }}_{\mathrm{s}}}\right)}^{\mathrm{1}\mathit{/}\mathrm{3}}\mathrm{·}$(1)The moment of inertia factor (MOI = I_p/(MR²)) is given as $\mathit{MOI}\mathrm{=}\frac{\mathrm{2}}{\mathrm{5}}\left(\frac{{\left(\mathit{\rho }\mathrm{-}{\mathit{\rho }}_{\mathrm{s}}\right)}^{\mathrm{5}\mathit{/}\mathrm{3}}}{\mathit{\rho }{\left({\mathit{\rho }}_{\mathrm{c}}\mathrm{-}{\mathit{\rho }}_{\mathrm{s}}\right)}^{\mathrm{2}\mathit{/}\mathrm{3}}}\mathrm{+}\frac{{\mathit{\rho }}_{\mathrm{s}}}{\mathit{\rho }}\right)\mathrm{·}$(2)Figure 1 shows that without additional constraints, plausible density values of the core and icy shell yield 0.3 < MOI < 0.4. The range of the MOI was used to estimate the moment of inertia of tri-axial Mimas, as shown below.

For a satellite in hydrostatic equilibrium, the MOI is related to the fluid Love number k_f, which describes the reaction of the satellite to a perturbing potential after all viscous stresses have relaxed (Munk & MacDonald 1960; Hubbard & Anderson 1978): $\mathit{MOI}\mathrm{=}\frac{\mathrm{2}}{\mathrm{3}}\left[\mathrm{1}\mathrm{-}\frac{\mathrm{2}}{\mathrm{5}}\sqrt{\frac{\mathrm{4}\mathrm{-}{\mathit{k}}_{\mathrm{f}}}{\mathrm{1}\mathrm{+}{\mathit{k}}_{\mathrm{f}}}}\right]\mathit{,}$(3)and the gravity coefficients C₂₂ and J₂ are determined from (Rappaport et al. 1997) $\begin{array}{ccc}{\mathit{C}}_{\mathrm{22}}& \mathrm{=}& \\ {\mathit{J}}_{\mathrm{2}}& \mathrm{=}& \end{array}$where , Ω being the spin velocity of Mimas, equal to its mean motion n because Mimas is in synchronous rotation. With the numerical values given in Table 1, we have q_r = 0.01854.

The differences between the three principal moments of inertia A < B < C are determined from the definitions of C₂₂ and J₂, i.e. $\begin{array}{ccc}\mathit{B}\mathrm{-}\mathit{A}& \mathrm{=}& \mathrm{4}{\mathit{C}}_{\mathrm{22}}\mathit{M}{\mathit{R}}^{\mathrm{2}}\mathit{,}\\ \mathit{C}\mathrm{-}\mathit{A}& \mathrm{=}& \mathrm{(}{\mathit{J}}_{\mathrm{2}}\mathrm{+}\mathrm{2}{\mathit{C}}_{\mathrm{22}}\mathrm{)}\mathit{M}{\mathit{R}}^{\mathrm{2}}\mathit{,}\\ \mathit{C}\mathrm{-}\mathit{B}& \mathrm{=}& \mathrm{(}{\mathit{J}}_{\mathrm{2}}\mathrm{-}\mathrm{2}{\mathit{C}}_{\mathrm{22}}\mathrm{)}\mathit{M}{\mathit{R}}^{\mathrm{2}}\mathit{.}\end{array}$The relationship between the mean moment of inertia $\mathit{I}\mathrm{=}\frac{\mathit{A}\mathrm{+}\mathit{B}\mathrm{+}\mathit{C}}{\mathrm{3}}$ and the polar moment of inertia C is $\mathit{C}\mathrm{=}\mathit{I}\mathrm{+}\frac{\mathrm{2}}{\mathrm{3}}{\mathit{J}}_{\mathrm{2}}\mathit{M}{\mathit{R}}^{\mathrm{2}}\mathit{.}$(6)We can then calculate all three moments of inertia A, B and C.

2.2. Nonhydrostatic shape

We here use the observed shape (a = 207.8 km, b = 196.7 km, c = 190.6 km, Thomas 2010) to calculate the moments of inertia of Mimas, assuming that the shape of the core is proportional to that of Mimas, i.e. we assume $\frac{{\mathit{a}}_{\mathrm{c}}}{\mathit{a}}\mathrm{=}\frac{{\mathit{b}}_{\mathrm{c}}}{\mathit{b}}\mathrm{=}\frac{{\mathit{c}}_{\mathrm{c}}}{\mathit{c}}\mathrm{=}\frac{{\mathit{R}}_{\mathrm{c}}}{\mathit{R}}\mathit{,}$(7)where a_c, b_c and c_c are the dimensions of the core, and R_c its mean radius (Eq. (1)).

A quadrature over the volume of the core and the shell gives $\begin{array}{ccc}{\mathit{C}}_{\mathrm{c}}& \mathrm{=}& \$\mathrm{core}{\mathit{\rho }}_{\mathrm{c}}\mathrm{(}{\mathit{x}}^{\mathrm{2}}\mathrm{+}{\mathit{y}}^{\mathrm{2}}\mathrm{)}\mathrm{d}\mathit{x}\hspace{0.17em}\mathrm{d}\mathit{y}\hspace{0.17em}\mathrm{d}\mathit{z}\\ & \mathrm{=}& \frac{\mathrm{4}}{\mathrm{15}}\mathit{\pi }{\mathit{a}}_{\mathrm{c}}{\mathit{b}}_{\mathrm{c}}{\mathit{c}}_{\mathrm{c}}\mathrm{(}{\mathit{a}}_{\mathrm{c}}^{\mathrm{2}}\mathrm{+}{{\mathit{b}}_{\mathrm{c}}^{\mathrm{2}}}^{\mathrm{)}}{\mathit{\rho }}_{\mathrm{c}}\\ & \mathrm{=}& \\ {\mathit{C}}_{\mathrm{s}}& \mathrm{=}& \$\mathrm{Mimas}{\mathit{\rho }}_{\mathrm{s}}\mathrm{(}{\mathit{x}}^{\mathrm{2}}\mathrm{+}{\mathit{y}}^{\mathrm{2}}\mathrm{)}\mathrm{d}\mathit{x}\hspace{0.17em}\mathrm{d}\mathit{y}\hspace{0.17em}\mathrm{d}\mathit{z}\\ & & \mathrm{-}{\$}_{\mathrm{core}}{\mathit{\rho }}_{\mathrm{s}}\mathrm{(}{\mathit{x}}^{\mathrm{2}}\mathrm{+}{\mathit{y}}^{\mathrm{2}}\mathrm{)}\mathrm{d}\mathit{x}\hspace{0.17em}\mathrm{d}\mathit{y}\hspace{0.17em}\mathrm{d}\mathit{z}\\ & \mathrm{=}& \end{array}$We then obtain C = C_c + C_s. The other moments of inertia A and B being obtained similarly, we have $\begin{array}{ccc}\mathit{A}& \mathrm{=}& \\ \mathit{B}& \mathrm{=}& \\ \mathit{C}& \mathrm{=}& \end{array}$We finally see that the ratio of the moments of inertia A/C = (b² + c²)/(a² + b²) and B/C = (a² + c²)/(a² + b²) are independent of the mean radius R_c and density ρ_c of the core, so every model of the internal structure of Mimas based on its observed shape (in neglecting the uncertainties on the radii a, b and c) will present the same rotational response.

The interior models considered in the present study are gathered in Table 2.

3. Computing the rotation of Mimas

In this section, Mimas is assumed to be a two-layer rigid body and the tidal contributions will be investigated in Sect. 5. Its rotation is highly constrained by the gravitational perturbation of Saturn, and so depends on the variations of the distance Mimas-Saturn. That is the reason why we must understand the orbital dynamics of Mimas before investigating its rotation.

3.1. The orbital dynamics of Mimas

Mimas is the smallest of the main Saturnian satellites, and also the closest to its parent planet and the rings. Discovered by Herschel in 1789, it is known since Struve (1891) to be in 2:1 mean-motion resonance with Tethys. More precisely, these two bodies are locked in an inclination-type resonance whose argument is 2λ₁ − 4λ₃ + ☊₁ + ☊₃, the subscript 1 standing for the satellite S-1 Mimas, 3 for S-3 Tethys, λ_i being the mean longitudes, and ☊_i the longitudes of the ascending nodes. This resonance tends to raise the inclinations of the satellites to  ≈ 1.5° for Mimas and  ≈ 1° for Tethys (Allan 1969), and stimulates librations of the resonant argument around 0 with an amplitude of  ≈ 95° and a period of  ≈ 70 years. The trapping of the system into this resonance can be explained in considering a non-null eccentricity for Tethys that induces secondary resonances that strongly enhance the capture probability (Champenois & Vienne 1999a,b).

It is convenient to work on a Fourier-type representation of the orbital motion of Mimas that allows us to identify every proper mode of the motion. The basic idea is that the variables describing the orbital motion of Mimas can be represented as a quasi-periodic series (and a slope for precessing angles such as the ascending node, the pericenter and the mean longitude), i.e. infinite but converging sums of trigonometric series. The arguments of these series can be expressed as an integer combination of a few proper modes of constant frequencies. The existence of these modes comes both from the KAM (Arnold 1963; Moser 1962) and the Nekhoroshev theories (Nekhoroshev 1977, 1979). The KAM theory states that for a quasi-integrable Hamiltonian system (i.e. like ℋ = ℋ₀ + ϵℋ₁ where ℋ₀ is an integrable Hamiltonian and ϵℋ₁ a small perturbation) verifying classical assumptions, the motion can be considered to be on invariant tori (i.e. with constant amplitudes and angles depending linearly on time) in action-angle coordinates. For a bigger perturbation the Nekhoroshev theory says that the invariant tori survive over a timescale that is exponentially long with respect to the invert of the amplitude of the perturbation ϵ, provided that the Hamiltonian of the system presents a property of steepness, that is, an extension of the convexity.

This representation is given by the TASS1.6 ephemerides (Vienne & Duriez 1995) where the orbital motion of Mimas can be described using the five proper modes λ, ω, φ, ζ and Φ. λ is the linear part of Mimas’ mean longitude, ω is the main oscillation mode of the librations of the resonant argument 2λ₁ − 4λ₃ + ☊₁ + ☊₃, ζ (called ρ₁ in Vienne & Duriez 1995) is the mean slope of λ₁ − 2λ₃, and φ − ζ and Φ − ζ are the mean slopes of the longitudes of the pericenter of Mimas and its ascending node, respectively. The values of the associated frequencies are gathered in Table 3.

3.2. Rotational model

As for most of the natural satellites of the solar system, Mimas is expected to follow the three Cassini Laws originally described for the Moon (Cassini 1693; Colombo 1966), i.e.

• 1.

  the Moon rotates uniformly about its polar axis with a rotational period equal to the mean sidereal period of its orbit about the Earth;

• 2.

  the inclination of the Moon’s equator to the ecliptic is a constant angle (approximately 1.5°);

• 3.

  the ascending node of the lunar orbit on the ecliptic coincides with the descending node of the lunar equator on the ecliptic. This law could also be expressed like this: the spin axis and the normals to the ecliptic and orbit plane remain coplanar.

For natural satellites, these laws can be rephrased this way: the rotation of the satellite is synchronous, its angular momentum has a nearly constant inclination on an inertial reference plane, and is located in the plane defined by the normal to the orbital plane and to the Laplace plane. The Laplace plane is the plane normal to the rotation axis of the orbital frame, i.e. it is defined with respect to the orbital precessional motion. It has the property to minimize the variations of the orbital inclinations. For satellites orbiting close to their planet, as is the case here, the equatorial plane of Saturn is so close to the Laplace plane that it can be used to describe the rotational dynamics.

Our rotational model is similar to the one already used in e.g. Noyelles et al. (2008); Noyelles (2010) for studying the rigid rotation of the Saturnian satellites Titan, Janus and Epimetheus.

We consider Mimas to be a rigid triaxial body whose matrices of inertia reads $\mathit{I}\mathrm{=}\left(\begin{array}{c}\\ \mathit{A}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathit{B}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathit{C}\end{array}\right)$(13)with A ≤ B ≤ C.

The dynamical model is a three-degree of freedom one. We used the Andoyer variables, which requires a decomposition with three references frames:

• 1.

  an inertial reference frame (e₁,e₂,e₃). We used the one in which the orbital ephemerides are given, i.e. mean Saturnian equator and mean equinox for the J2000.0 epoch;

• 2.

  a frame (n₁,n₂,n₃) bound to the angular momentum of Mimas;

• 3.

  a frame (f₁,f₂,f₃) rigidly linked to Mimas.

We first used Andoyer’s variables (Andoyer 1926; Deprit 1967), which are based on two linked sets of Euler’s angles. The first set (h,K,g) locates the position of the angular momentum in the first frame (e₁,e₂,e₃), while the second one, (g,J,l), locates the body frame (f₁,f₂,f₃) in the second frame tied to the angular momentum (see Fig. 2).

The canonical set of Andoyer’s variables consists of the three angular variables l,g,h and their conjugated momenta L,G,H defined by the norm G of the angular momentum and two of its projections: $\begin{array}{c}\\ \mathit{l,}& & \mathit{L}\mathrm{=}\mathit{G}\cos \mathit{J,}\\ \mathit{g,}& & \mathit{G,}\\ \mathit{h,}& & \mathit{H}\mathrm{=}\mathit{G}\cos \mathit{K}\mathit{.}\end{array}$(14)Unfortunately, these variables present two singularities: when J = 0 (i.e., the angular momentum is colinear to f₃), l and g are undefined, and when K = 0 (i.e., when Mimas’ principal axis of inertia is perpendicular to its orbital plane), h and g are undefined. That is the reason why we used the modified Andoyer’s variables: $\begin{array}{c}\\ \mathit{p}\mathrm{=}\mathit{l}\mathrm{+}\mathit{g}\mathrm{+}\mathit{h,}& & \mathit{P}\mathrm{=}\frac{\mathit{G}}{\mathit{nC}}\mathit{,}\\ \mathit{r}\mathrm{=}\mathrm{-}\mathit{h,}& & \mathrm{ℛ}\mathrm{=}\frac{\mathit{G}\mathrm{-}\mathit{H}}{\mathit{nC}}\mathrm{=}\mathit{P}\mathrm{(}\mathrm{1}\mathrm{-}\cos \mathit{K}\mathrm{)}\mathit{,}\\ & & \mathrm{=}\mathrm{2}\mathit{P}{\sin }^{\mathrm{2}}\frac{\mathit{K}}{\mathrm{2}}\mathit{,}\\ {\mathit{\xi }}_{\mathit{q}}\mathrm{=}\sqrt{\frac{\mathrm{2}\mathit{Q}}{\mathit{nC}}}\sin \mathit{q,}& & {\mathit{\eta }}_{\mathit{q}}\mathrm{=}\sqrt{\frac{\mathrm{2}\mathit{Q}}{\mathit{nC}}}\cos \mathit{q,}\end{array}$(15)where n is the body’s mean orbital motion, q =  −l, and $\mathit{Q}\mathrm{=}\mathit{G}\mathrm{-}\mathit{L}\mathrm{=}\mathit{G}\mathrm{(}\mathrm{1}\mathrm{-}\cos \mathit{J}\mathrm{)}\mathrm{=}\mathrm{2}\mathit{G}{\sin }^{\mathrm{2}}\frac{\mathit{J}}{\mathrm{2}}$. With these new variables, the singularity on l has been dropped. Using these variables has a great mathematical interest, because they are canonical, so they simplify an analytical study of the system, as was done in the previous works mentioned above. Our study here is essentially numerical, but we keep these variables to be consistent with previous studies. We derive other output variables below, that are more relevant from a physical point of view.

In these variables, the kinetic energy $\mathit{T}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}{\omega }\mathrm{·}{G}$ of the system reads $\begin{array}{ccc}\mathit{T}& \mathrm{=}& \frac{\mathit{n}{\mathit{P}}^{\mathrm{2}}}{\mathrm{2}}\mathrm{+}\frac{\mathit{n}}{\mathrm{8}}\mathrm{[}\mathrm{4}\mathit{P}\mathrm{-}{\mathit{\xi }}_{\mathit{q}}^{\mathrm{2}}\mathrm{-}{{\mathit{\eta }}_{\mathit{q}}^{\mathrm{2}}}^{\mathrm{]}}\\ & & \mathrm{\times }\left[\frac{{\mathit{\gamma }}_{\mathrm{1}}\mathrm{+}{\mathit{\gamma }}_{\mathrm{2}}}{\mathrm{1}\mathrm{-}{\mathit{\gamma }}_{\mathrm{1}}\mathrm{-}{\mathit{\gamma }}_{\mathrm{2}}}{\mathit{\xi }}_{\mathit{q}}^{\mathrm{2}}\mathrm{+}\frac{{\mathit{\gamma }}_{\mathrm{1}}\mathrm{-}{\mathit{\gamma }}_{\mathrm{2}}}{\mathrm{1}\mathrm{-}{\mathit{\gamma }}_{\mathrm{1}}\mathrm{+}{\mathit{\gamma }}_{\mathrm{2}}}{\mathit{\eta }}_{\mathit{q}}^{\mathrm{2}}\right]\end{array}$(16)with ${\mathit{\gamma }}_{\mathrm{1}}\mathrm{=}\frac{\mathrm{2}\mathit{C}\mathrm{-}\mathit{A}\mathrm{-}\mathit{B}}{\mathrm{2}\mathit{C}}\mathrm{=}{\mathit{J}}_{\mathrm{2}}\frac{\mathit{M}{\mathit{R}}^{\mathrm{2}}}{\mathit{C}}$(17)and ${\mathit{\gamma }}_{\mathrm{2}}\mathrm{=}\frac{\mathit{B}\mathrm{-}\mathit{A}}{\mathrm{2}\mathit{C}}\mathrm{=}\mathrm{2}{\mathit{C}}_{\mathrm{22}}\frac{\mathit{M}{\mathit{R}}^{\mathrm{2}}}{\mathit{C}}\mathrm{·}$(18)In these last three formulae, ω is the instantaneous vector of rotation, M is the mass of Mimas, R its mean radius, and J₂ and C₂₂ the two classical normalized gravitational coefficients related to the oblateness and equatorial ellipticity of the considered body, respectively.

The gravitational disturbing potential caused by an oblate perturber p reads (Henrard 2005c) ${\mathit{V}}_{\mathrm{p}}\mathrm{=}{\mathit{V}}_{\mathit{p}\mathrm{1}}\mathrm{+}{\mathit{V}}_{\mathit{p}\mathrm{2}}$(19)with ${\mathit{V}}_{\mathit{p}\mathrm{1}}\mathrm{=}\mathrm{-}\frac{\mathrm{3}}{\mathrm{2}}\mathit{C}\frac{\mathrm{𝒢}{\mathit{M}}_{\mathrm{p}}}{{\mathit{d}}_{\mathrm{p}}^{\mathrm{3}}}\mathrm{[}{\mathit{\gamma }}_{\mathrm{1}}\mathrm{(}{\mathit{x}}_{\mathrm{p}}^{\mathrm{2}}\mathrm{+}{{\mathit{y}}_{\mathrm{p}}^{\mathrm{2}}}^{\mathrm{)}}\mathrm{+}{\mathit{\gamma }}_{\mathrm{2}}{\mathrm{(}{\mathit{x}}_{\mathrm{p}}^{\mathrm{2}}\mathrm{-}{{\mathit{y}}_{\mathrm{p}}^{\mathrm{2}}}^{\mathrm{)}}}^{\mathrm{]}}$(20)and $\begin{array}{ccc}{\mathit{V}}_{\mathit{p}\mathrm{2}}& \mathrm{=}& \mathrm{-}\frac{\mathrm{15}}{\mathrm{4}}\mathit{C}{\mathit{J}}_{\mathrm{2}\mathit{p}}\frac{\mathrm{𝒢}{\mathit{M}}_{\mathrm{p}}}{{\mathit{d}}_{\mathrm{p}}^{\mathrm{3}}}{\left(\frac{{\mathit{R}}_{\mathrm{p}}}{{\mathit{d}}_{\mathrm{p}}}\right)}^{\mathrm{2}}\\ & & \end{array}$(21)where is the gravitational constant, M_p the mass of the perturber, J_2p its J₂, R_p its mean radius, d_p the distance between the perturber’s and Mimas’ centers of mass, and x_p and y_p the two first components of the unit vector pointing to the center of mass of the perturber, from the center of mass of the body, in the reference frame (f₁,f₂,f₃). V_p1 expresses the perturbation caused by a pointmass perturber, while V_p2 represents the perturbation caused by its J₂, assuming that the body is in the equatorial plane of the perturber. As shown in (Henrard 2005c), it is a good approximation if the sine of the angle between Saturn’s equatorial plane and the orbit is small. In the case of Mimas, this angle (i.e. Mimas’ orbital inclination) is  ≈ 1.5° ≈ 2.6 × 10^-2 rad, so we can consider that its sine is always smaller than 3 × 10^-2. This assertion also assumes that the obliquity of Mimas is very low, which we will check in this study.

Usually the orbital ephemerides give us the location of the perturber in the inertial frame, therefore we have to perform five rotations to convert the coordinates from the inertial frame to (f₁,f₂,f₃). More precisely, if we name (x_i,y_i,z_i)^T the unit vector locating the perturber’s center of mass in the inertial frame, we have $\left(\begin{array}{c}\\ {\mathit{x}}_{\mathrm{p}}\\ {\mathit{y}}_{\mathrm{p}}\\ {\mathit{z}}_{\mathrm{p}}\end{array}\right)\mathrm{=}{\mathit{R}}_{\mathrm{3}}\mathrm{(}\mathrm{-}\mathit{l}\mathrm{)}{\mathit{R}}_{\mathrm{1}}\mathrm{(}\mathrm{-}\mathit{J}\mathrm{)}{\mathit{R}}_{\mathrm{3}}\mathrm{(}\mathrm{-}\mathit{g}\mathrm{)}{\mathit{R}}_{\mathrm{1}}\mathrm{(}\mathrm{-}\mathit{K}\mathrm{)}{\mathit{R}}_{\mathrm{3}}\mathrm{(}\mathrm{-}\mathit{h}\mathrm{)}\left(\begin{array}{c}\\ {\mathit{x}}_{\mathit{i}}\\ {\mathit{y}}_{\mathit{i}}\\ {\mathit{z}}_{\mathit{i}}\end{array}\right)$(22)with ${\mathit{R}}_{\mathrm{3}}\mathrm{\left(}\mathit{\phi }\mathrm{\right)}\mathrm{=}\left(\begin{array}{c}\\ \cos \mathit{\phi }& \mathrm{-}\sin \mathit{\phi }& \mathrm{0}\\ \sin \mathit{\phi }& \cos \mathit{\phi }& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{1}\end{array}\right)$(23)and ${\mathit{R}}_{\mathrm{1}}\mathrm{\left(}\mathit{\phi }\mathrm{\right)}\mathrm{=}\left(\begin{array}{c}\\ \mathrm{1}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \cos \mathit{\phi }& \mathrm{-}\sin \mathit{\phi }\\ \mathrm{0}& \sin \mathit{\phi }& \cos \mathit{\phi }\end{array}\right)\mathit{.}$(24)Finally, the total Hamiltonian of the problem reads $\begin{array}{ccc}\mathit{H}& \mathrm{=}& \frac{\mathit{n}{\mathit{P}}^{\mathrm{2}}}{\mathrm{2}}\mathrm{+}\frac{\mathit{n}}{\mathrm{8}}\mathrm{[}\mathrm{4}\mathit{P}\mathrm{-}{\mathit{\xi }}_{\mathit{q}}^{\mathrm{2}}\mathrm{-}{{\mathit{\eta }}_{\mathit{q}}^{\mathrm{2}}}^{\mathrm{]}}\\ & & \mathrm{\times }\left[\frac{{\mathit{\gamma }}_{\mathrm{1}}\mathrm{+}{\mathit{\gamma }}_{\mathrm{2}}}{\mathrm{1}\mathrm{-}{\mathit{\gamma }}_{\mathrm{1}}\mathrm{-}{\mathit{\gamma }}_{\mathrm{2}}}{\mathit{\xi }}_{\mathit{q}}^{\mathrm{2}}\mathrm{+}\frac{{\mathit{\gamma }}_{\mathrm{1}}\mathrm{-}{\mathit{\gamma }}_{\mathrm{2}}}{\mathrm{1}\mathrm{-}{\mathit{\gamma }}_{\mathrm{1}}\mathrm{+}{\mathit{\gamma }}_{\mathrm{2}}}{\mathit{\eta }}_{\mathit{q}}^{\mathrm{2}}\right]\\ & & \mathrm{-}\frac{\mathrm{3}}{\mathrm{2}\mathit{n}}\frac{\mathrm{𝒢}{\mathit{M}}_{♄}}{{\mathit{d}}_{♄}^{\mathrm{3}}}\left(\mathrm{1}\mathrm{+}\frac{\mathrm{5}}{\mathrm{2}}{\mathit{J}}_{\mathrm{2}}{\left(\frac{{\mathit{R}}_{♄}}{{\mathit{d}}_{♄}}\right)}^{\mathrm{2}}\right)\\ & & \end{array}$(25)where the index ♄ stands for Saturn. We will use this Hamiltonian for a numerical study of the rotation. An analytical study can show that the Hamiltonian (25) can be reduced to $\mathrm{\mathscr{H}}\mathrm{\left(}\mathit{u,v,w,U,V,W}\mathrm{\right)}\mathrm{=}{\mathit{\omega }}_{\mathit{u}}\mathit{U}\mathrm{+}{\mathit{\omega }}_{\mathit{v}}\mathit{V}\mathrm{+}{\mathit{\omega }}_{\mathit{w}}\mathit{W}\mathrm{+}\mathrm{𝒫}\mathrm{\left(}\mathit{u,v,w,U,V,W}\mathrm{\right)}\mathit{,}$(26)where represents a perturbation, and the three constants ω_u, ω_v and ω_w are the periods of the free oscillations around the equilibrium defined by the Cassini laws. This last Hamiltonian is obtained after several canonical transformations, the first one consisting of expressing the resonant arguments σ = p − λ + π and ρ = r + ☊ respectively associated with the 1:1 spin-orbit resonance and with the orientation of the angular momentum, λ and ☊ being the orbital variables defined above. The complete calculation is beyond the scope of this paper, the reader can find details in Henrard (2005a,b), and Noyelles et al. (2008).

3.3. A numerical study

To integrate the system numerically, we first express the coordinates of the perturber (x_♄,y_♄) with the numerical ephemerides and the rotations given in Eq. (22) in the body frame (f₁,f₂,f₃). As explained above, the ephemerides are given by the TASS1.6 ephemerides (Vienne & Duriez 1995). This way, we obtain coordinates depending of the canonical variables. Then we derive the equations coming from the Hamiltonian (25) $\begin{array}{ccc}\frac{\mathrm{d}\mathit{p}}{\mathrm{d}\mathit{t}}\mathrm{=}\frac{\mathit{\partial H}}{\mathit{\partial P}}\mathit{,}& & \frac{\mathrm{d}\mathit{P}}{\mathrm{d}\mathit{t}}\mathrm{=}\mathrm{-}\frac{\mathit{\partial H}}{\mathit{\partial p}}\mathit{,}\\ \frac{\mathrm{d}\mathit{r}}{\mathrm{d}\mathit{t}}\mathrm{=}\frac{\mathit{\partial H}}{\mathit{\partial R}}\mathit{,}& & \frac{\mathrm{d}\mathit{R}}{\mathrm{d}\mathit{t}}\mathrm{=}\mathrm{-}\frac{\mathit{\partial H}}{\mathit{\partial r}}\mathit{,}\\ \frac{\mathrm{d}{\mathit{\xi }}_{\mathit{q}}}{\mathrm{d}\mathit{t}}\mathrm{=}\frac{\mathit{\partial H}}{\mathit{\partial }{\mathit{\eta }}_{\mathit{q}}}\mathit{,}& & \end{array}$(27)We integrated over 200 years using the Adams-Bashforth-Moulton 10th order predictor-corrector integrator. The solutions consist of two parts, the forced one, directly caused by the perturbation, and the free one, which depends on the initial conditions. The initial conditions should be as close as possible to the exact equilibrium, which is assumed to be the Cassini state 1 in 1:1 spin-orbit resonance, to have low amplitudes of the free librations. For that, we used the iterative algorithm NAFFO (Noyelles et al. 2011) to remove the free librations from the initial conditions after they were identified by frequency analysis.

The frequency analysis algorithm we used is based on Laskar’s original idea, named NAFF for numerical analysis of the fundamental frequencies (see for instance Laskar (1993) for the method, and Laskar (2005) for the convergence proofs). It aims at identifying the coefficients a_k and ω_k of a complex signal f(t) obtained numerically over a finite time span  [ − T;T]  and verifying $\mathit{f}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{\approx }\sum _{\mathit{k}\mathrm{=}\mathrm{1}}^{\mathit{n}}{\mathit{a}}_{\mathit{k}}\exp \mathrm{\left(}\mathrm{i}{\mathit{\omega }}_{\mathit{k}}\mathit{t}\mathrm{\right)}\mathit{,}$(28)where ω_k are real frequencies and a_k complex coefficients. If the signal f(t) is real, its frequency spectrum is symmetric and the complex amplitudes associated with the frequencies ω_k and  − ω_k are complex conjugates. The frequencies and amplitudes associated are found with an iterative scheme. To determine the first frequency ω₁, one searches for the maximum of the amplitude of $\mathit{\phi }\mathrm{\left(}\mathit{\omega }\mathrm{\right)}\mathrm{=}\mathrm{⟨}\mathit{f}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathit{,}\exp \mathrm{\left(}\mathrm{i}\mathit{\omega t}\mathrm{\right)}\mathrm{⟩}\mathit{,}$(29)where the scalar product  ⟨ f(t),g(t) ⟩  is defined by $\mathrm{⟨}\mathit{f}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathit{,g}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{⟩}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}\mathit{T}}{\int }_{\mathrm{-}\mathit{T}}^{\mathit{T}}\mathit{f}\mathrm{\left(}\mathit{t}\mathrm{\right)}\overline{\mathit{g}\mathrm{\left(}\mathit{t}\mathrm{\right)}}\mathit{\chi }\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{d}\mathit{t,}$(30)and where χ(t) is a weight function, i.e. a positive function with $\frac{\mathrm{1}}{\mathrm{2}\mathit{T}}{\int }_{\mathrm{-}\mathit{T}}^{\mathit{T}}\mathit{\chi }\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{d}\mathit{t}\mathrm{=}\mathrm{1.}$(31)Once the first periodic term exp(iω₁t) is found, its complex amplitude a₁ is obtained by orthogonal projection, and the process is started again on the remainder f₁(t) = f(t) − a₁exp(iω₁t). The algorithm stops when two detected frequencies are too close to each other, which alters their determinations, or when the number of detected terms reaches a maximum set by the user. This algorithm is very efficient, except when two frequencies are too close to each other. In that case, the algorithm is not confident in its accuracy and stops. When the difference between two frequencies is larger than twice the frequency associated with the length of the total time interval, the determination of each fundamental frequency is not perturbed by the other ones. Although the iterative method suggested by Champenois (1998) allows one to reduce this distance, some difficulties remain when the frequencies are too close to each other.

3.4. Outputs

To deliver theories of rotation that can be easily compared with observations, we chose to express our results in the following variables:

• longitudinal librations;

• latitudinal librations;

• orbital obliquity ϵ (the orientation of the angular momentum of Mimas with respect to the normal to the instantaneous orbital plane);

• motion of the rotation axis about the pole axis.

There are at least two ways to define the longitudinal librations. We can for instance consider the librations about the exact synchronous rotation, i.e. p −  ⟨ n ⟩ t, usually called physical librations. In this case, we used for  ⟨ n ⟩  the frequency associated with the proper mode λ, i.e. Mimas’ mean longitude. Another way to consider the longitudinal librations is to work on the librations about the Mimas-Saturn direction. We will call these librations tidal librations because they represent the misalignment of the tidal bulge of the satellite. The difference between these two librations is known as optical librations, only caused by the orbital motion of Mimas around Saturn. The reader can find graphical descriptions of these librations in Murray & Dermott (1999, Fig. 5.16).

The latitudinal librations are the north-south librations of the large axis of the considered body in the saturnocentric reference frame that follows the orbital motion of the body. They are analogous to the tidal librations that are the east-west librations. To obtain the tidal longitudinal librations and the latitudinal librations, we first express the unit vector f₁ (i.e. the direction of Mimas’ long axis) in the inertial frame (e₁,e₂,e₃). From (Eq. (22)) and the definitions of the Andoyer modified variables (Eq. (15)), we obtain: $\begin{array}{ccc}{{f}}_{\mathbf{1}}& \mathrm{=}& \mathrm{\left(}\cos \mathit{r}\mathrm{\right(}\cos \mathrm{(}\mathit{p}\mathrm{+}\mathit{r}\mathrm{-}\mathit{l}\mathrm{)}\cos \mathit{l}\mathrm{-}\sin \mathrm{(}\mathit{p}\mathrm{+}\mathit{r}\mathrm{-}\mathit{l}\mathrm{)}\cos \mathit{J}\sin \mathit{l}\mathrm{)}\\ & & \mathrm{+}\sin \mathit{r}\mathrm{\left(}\cos \mathit{K}\mathrm{\right(}\sin \mathrm{(}\mathit{p}\mathrm{+}\mathit{r}\mathrm{-}\mathit{l}\mathrm{)}\cos \mathit{l}\\ & & \mathrm{+}\cos \mathrm{(}\mathit{p}\mathrm{+}\mathit{r}\mathrm{-}\mathit{l}\mathrm{)}\cos \mathit{J}\sin \mathit{l}\mathrm{)}\mathrm{-}\sin \mathit{K}\sin \mathit{J}\sin \mathit{l}\mathrm{)}\mathrm{)}{{e}}_{\mathbf{1}}\\ & & \mathrm{+}\mathrm{(}\mathrm{-}\sin \mathit{r}\mathrm{(}\cos \mathrm{(}\mathit{p}\mathrm{+}\mathit{r}\mathrm{-}\mathit{l}\mathrm{)}\cos \mathit{l}\mathrm{-}\sin \mathrm{(}\mathit{p}\mathrm{+}\mathit{r}\mathrm{-}\mathit{l}\mathrm{)}\cos \mathit{J}\sin \mathit{l}\mathrm{)}\\ & & \mathrm{+}\cos \mathit{r}\mathrm{\left(}\cos \mathit{K}\mathrm{\right(}\sin \mathrm{(}\mathit{p}\mathrm{+}\mathit{r}\mathrm{-}\mathit{l}\mathrm{)}\cos \mathit{l}\\ & & \mathrm{+}\cos \mathrm{(}\mathit{p}\mathrm{+}\mathit{r}\mathrm{-}\mathit{l}\mathrm{)}\cos \mathit{J}\sin \mathit{l}\mathrm{)}\mathrm{-}\sin \mathit{K}\sin \mathit{J}\sin \mathit{l}\mathrm{)}\mathrm{)}{{e}}_{\mathbf{2}}\\ & & \mathrm{+}\mathrm{\left(}\sin \mathit{K}\mathrm{\right(}\sin \mathrm{(}\mathit{p}\mathrm{+}\mathit{r}\mathrm{-}\mathit{l}\mathrm{)}\cos \mathit{l}\mathrm{+}\cos \mathrm{(}\mathit{p}\mathrm{+}\mathit{r}\mathrm{-}\mathit{l}\mathrm{)}\cos \mathit{J}\sin \mathit{l}\mathrm{)}\\ & & \end{array}$(32)The tidal longitudinal librations ψ and the latitudinal ones η are found by $\mathit{\psi }\mathrm{=}{t}\mathrm{·}{{f}}_{\mathbf{1}}$(33)and $\mathit{\eta }\mathrm{=}{n}\mathrm{·}{{f}}_{\mathbf{1}}\mathit{,}$(34)where n is the unit vector normal to the orbit plane, and t the tangent to the trajectory. We compute these last two vectors by ${n}\mathrm{=}\frac{{x}\mathrm{\times }{v}}{\mathrm{\left|}\mathrm{\right|}{x}\mathrm{\times }{v}\mathrm{\left|}\mathrm{\right|}}$(35)and ${t}\mathrm{=}\frac{{n}\mathrm{\times }{x}}{\mathrm{\left|}\mathrm{\right|}{n}\mathrm{\times }{x}\mathrm{\left|}\mathrm{\right|}}\mathit{,}$(36)where x is the position vector of the body, and v its velocity.

Finally, the motion of the rotation axis about the pole is derived from the wobble J, it is given by the two variables Q₁ and Q₂ defined as ${\mathit{Q}}_{\mathrm{1}}\mathrm{=}\sin \mathit{J}\sin \mathit{l}\left(\mathrm{1}\mathrm{+}\frac{{\mathit{J}}_{\mathrm{2}}\mathrm{+}\mathrm{2}{\mathit{C}}_{\mathrm{22}}}{\mathit{C}}\right)$(37)and ${\mathit{Q}}_{\mathrm{2}}\mathrm{=}\sin \mathit{J}\cos \mathit{l}\left(\mathrm{1}\mathrm{+}\frac{{\mathit{J}}_{\mathrm{2}}\mathrm{-}\mathrm{2}{\mathit{C}}_{\mathrm{22}}}{\mathit{C}}\right)\mathit{,}$(38)they are the first two components of the unit vector pointing at the instantaneous north pole of Mimas’ rotation axis in the body frame of Mimas. These quantities are finally multiplied by the polar radius of the satellite, i.e. 190.6 km (Thomas 2010) to obtain a deviation in meters.

4. Results

We here present the outputs of our numerical study of the rotation of Mimas. We first give the example of a non-hydrostatic model of Mimas based on its observed shape, then we compare the results with the rotational response of the first 22 models of Table 2, obtained considering Mimas to be in hydrostatic equilibrium.

4.1. Non-hydrostatic Mimas based on its shape

As already mentioned, this case is unique, because changes in the size of the core do not affect the ratios of the moments of inertia A/C, and B/C, and the coefficients γ₁ and γ₂ (Eqs. (17) and (18)). Consequently, there is a unique rotational behaviour of Mimas for any homogenous or two-layer model using this specific model based on the observed shape.

The free librations around the equilibrium are assumed to be damped, it is in any case important to know their frequencies ω_u, ω_v and ω_w (or periods T_u, T_v and T_w) because they characterize the way the system reacts to external sinusoidal excitations, which are here caused by the variations of the distance between the Sun and Mimas.

The frequencies of the free librations are listed in Table 4. The proper mode u roughly represents the free longitudinal librations, v the free librations of the obliquity, and w the wobble, i.e. the free polar motion of Mimas. These frequencies were deduced from the frequency analysis of the modified Andoyer variables (cf. Eq. (15)).

The proper modes involved in the Fourier representations of the librations of Mimas are the forced modes caused by the orbital motion of Mimas around Saturn (cf. Table 3) and the free ones (Table 4). The arguments of the sinusoidal components of the quasi-periodic decompositions of the rotation variables are integer combinations of these proper modes. If we consider that the free librations are damped, the solutions should be only composed of the forced modes.

The forced librations of Mimas modelled from its observed shape are given in Tables 5 to 8. These tables give the solutions in the form of periodic time series in cosines. Evidently, the main difference between the physical and the tidal librations is in the presence in the physical librations of a long-period term ( ≈ 70 years) with a high amplitude ( ≈ 43°, i.e.  ≈ 86° peak-to-peak) caused by the librations of the argument of the orbital resonance between Mimas and Tethys. As explained in Rambaux et al. (2010, 2011), the amplitude of the long-period librations are equal to the magnitude of the orbital perturbations because at long period the body is oriented towards the central planet. As a consequence, by analysing the tidal librations, the long period librations vanish. There is also a considerable difference in the amplitude given for the tidal and physical longitudinal librations. As explained above, this difference is caused by optical librations, with amplitude 2e ≈ 3.8 × 10^-2 rad  ≈ 2.2°.

The latitudinal librations of Mimas (Table 7) are significantly smaller ( ≈ 2 arcmin vs. 2.5° for the tidal longitudinal librations), and so could hardly be used in the framework of observations of the rotation of Mimas (except if there are free oscillations due to a recent unexpected excitation). The mean obliquity of Mimas (Table 8) is of the same order of magnitude.

4.2. For a hydrostatic Mimas

We performed the same numerical study for the 22 hydrostatic configurations of Mimas given in Table 2. The results are gathered in Table 9. In this table, the amplitudes of the tidal and physical longitudinal librations indicated are related to the mode λ − φ + ζ (period: 0.944898 d), while the latitudinal ones are related to the mode λ + φ − ζ (period: 0.939962 d). The main physical reason for these librations is the variation of the distance Mimas-Saturn during an orbital period.

To make the results more readable, we present them graphically in Fig. 3. The plots present a clear dependency of the amplitudes of librations on the densities of the core and the shell. We can notice in particular that the longitudinal (i.e. tidal and physical) librations have a larger amplitude when the density of the core is lower, which is because a concentration of the mass in the core lowers the moments of inertia of the body, and so tends to limit its amplitude of response to solicitations. Finally we can see that the dependency on ρ_c is small for ρ_s = 1100 kg/m³, it is because in this case, ρ_s is close to the mean density of Mimas (i.e. 1150.03 kg/m³), as a consequence the core is small and Mimas is close to be homogeneous.

4.3. A small polar motion

We present the forced polar motion of Mimas (i.e. after removal of the free wobble) as a sum of a trigonometric series (Table 10). We can see that this motion is expected to be small, the highest amplitude being  ≈ 15 m. The sum of all these amplitudes can reach 40 m, so we can consider these 40 m as the upper bound of the polar motion. An analysis of the polar motions for the different hydrostatic Mimas does not exhibit significant differences.

5. Tidal dissipation

This section is dedicated to study the influence of the tidal torque on the rotational motion of Mimas. We introduce the tidal torque in a Lagrangian formalism and follow the approach of Williams et al. (2001) that was also recently used in Rambaux et al. (2010) and Robutel et al. (2011). The starting equation is the angular momentum equation $\frac{\mathrm{d}{G}}{\mathrm{d}\mathit{t}}\mathrm{+}{\omega }\mathrm{\wedge }{G}\mathrm{=}{T}\mathit{,}$(39)where ω is the angular velocity vector, the angular momentum G = Iω with I the tensor of inertia, and T is the external gravitational torque expressed as ${T}\mathrm{=}\frac{\mathrm{3}\mathrm{𝒢}{\mathit{M}}_{♄}}{{\mathit{r}}^{\mathrm{3}}}{u}\mathrm{\wedge }\mathit{I}{u}\mathit{,}$(40)where u is the cosine director of Saturn in the reference frame tied to Mimas, and M_♄ its mass.

The dissipation is caused by the tidal and centrifugal potentials that deform the satellite. In this case, the tensor of inertia I becomes a constant plus a time-variable part resulting from the deformation. The time-variable part does not react instantaneously and therefore presents a time delay δt characteristics of the rheological properties of the body (see Sect. 2).

In addition, the dynamical equation of the rotational motion Eq. (39) may be linearized by using the synchronous spin-orbit resonance of the body implying that ω₁,ω₂ ≪ ω₃ ~ n and u₂,u₃ ≪ u₁ ~ 1 where u₁,u₂,u₃ are the coordinates of the cosine director along the principal axis of inertia of Mimas.

By using these approximations and focusing on the libration in longitude, the main tidal torque is expressed as (Williams et al. 2001) ${\mathit{T}}_{\mathit{t}}\mathrm{=}\mathrm{-}{\mathit{k}}_{\mathrm{2}}{\mathit{R}}^{\mathrm{5}}\frac{\mathrm{3}\mathrm{𝒢}{\mathit{M}}_{♄}^{\mathrm{2}}}{{\mathit{a}}^{\mathrm{6}}}\mathrm{(}{\mathit{U}}_{\mathrm{11}}{\mathit{U}}_{\mathrm{12}}^{\mathrm{\ast }}\mathrm{-}{\mathit{U}}_{\mathrm{12}}{\mathit{U}}_{\mathrm{11}}^{\mathrm{\ast }}\mathrm{)}\mathit{,}$(41)where ${\mathit{U}}_{\mathit{ij}}\mathrm{=}{{}^{\mathrm{\right(}}\frac{\mathit{a}}{\mathit{r}}^{\mathrm{\left)}}}^{\mathrm{3}}{\mathit{u}}_{\mathit{i}}{\mathit{u}}_{\mathit{j}}$ and the star indicates the time-delay part.

Then, we used the same approach as in Rambaux et al. (2010). We introduce the rotation angle ϕ similar to the sum of the Andoyer angles l + g because the polar motion is small, as shown in Fig. 2, where J and θ (the nutation angle) are small. The libration angle γ is defined as ϕ = M + γ representing the physical libration in longitude of the body. The cosine director u₂ ~ s − γ is of the order of the difference between the orbital variation s and the physical libration γ. We note that u₂ corresponds to the tidal libration ψ defined in Sect. 3.4 and their amplitude is small as shown in Table 5. The quantity s, the orbital variation, is defined as the difference between the true and the mean longitude of Mimas and represents the oscillation of the orbital longitude of Mimas that may be expressed in Fourier series as $\mathit{s}\mathrm{=}\sum _{\mathit{i}}{\mathit{H}}_{\mathit{i}}\sin \mathrm{(}{\mathit{\omega }}_{\mathit{i}}\mathit{t}\mathrm{+}{\mathit{\alpha }}_{\mathit{i}}\mathrm{)}\mathit{.}$(42)Then, by developing ${\mathit{u}}_{\mathrm{1}}^{\mathrm{\ast }}$ and ${\mathit{u}}_{\mathrm{2}}^{\mathrm{\ast }}$ in Taylor series for small δt, the dynamical equation becomes $\mathit{C}\begin{array}{}\mathrm{¨}\\ \mathit{\gamma }\end{array}\mathrm{+}\frac{\mathrm{3}}{\mathrm{2}}\mathrm{(}\mathit{B}\mathrm{-}\mathit{A}\mathrm{)}\frac{\mathrm{𝒢}\mathit{m}}{{\mathit{r}}^{\mathrm{3}}}\sin \mathrm{2}\mathrm{(}\mathit{\gamma }\mathrm{-}\mathit{s}\mathrm{)}\mathrm{=}\mathrm{-}{\mathit{k}}_{\mathrm{2}}{\mathit{R}}^{\mathrm{5}}\frac{\mathrm{3}\mathrm{𝒢}{\mathit{M}}_{♄}^{\mathrm{2}}}{{\mathit{a}}^{\mathrm{6}}}\mathit{\delta t}\mathrm{(}\mathit{\gamma ̇}\mathrm{-}\mathit{ṡ}\mathrm{)}\mathit{.}$(43)As shown in the previous section, the quantity γ − s is always small (see Table 5) allowing us to simplify the sine function by its angle. In addition, the eccentricity of Mimas is low and so a/r is equal to 1 at first order in eccentricity. Finally, we obtain a forced dissipative harmonic oscillator with the frequency ${\mathit{\omega }}_{\mathrm{0}}\mathrm{=}\mathit{n}\sqrt{\mathrm{3}\mathrm{(}\mathit{B}\mathrm{-}\mathit{A}\mathrm{)}\mathit{/}\mathit{C}}$ and the dissipative rate λ expressed as $\mathrm{2}\mathit{\lambda }\mathrm{=}\frac{\mathrm{3}{\mathit{k}}_{\mathrm{2}}{\mathit{R}}^{\mathrm{3}}}{\mathit{C}}\frac{{\mathit{n}}^{\mathrm{4}}}{\mathrm{𝒢}\mathit{m}}\mathit{\delta t,}$(44)m being the mass of the satellite and k₂ is the Love number of Mimas.

In the conservative case, the amplitude of terms associated with the long period is almost equal to the magnitude of the oscillation s. The solution may be expressed as $\mathit{\gamma }\mathrm{=}{\mathit{A}}_{\mathit{d}}\sin \mathrm{(}{\mathit{\omega }}_{\mathit{d}}\mathit{t}\mathrm{+}{\mathit{\phi }}_{\mathit{d}}\mathrm{)}{\mathrm{e}}^{\mathrm{-}\mathit{\lambda t}}\mathrm{+}\sum _{\mathit{i}}{\mathit{x}}_{\mathit{i}}\cos \mathrm{(}{\mathit{\omega }}_{\mathit{i}}\mathit{t}\mathrm{+}{\mathit{\alpha }}_{\mathit{i}}\mathrm{)}\mathrm{+}{\mathit{y}}_{\mathit{i}}\sin \mathrm{(}{\mathit{\omega }}_{\mathit{i}}\mathit{t}\mathrm{+}{\mathit{\alpha }}_{\mathit{i}}\mathrm{)}\mathit{,}$(45)where A_d and φ_d are constants of integration. The first term decays with time scale 1/λ and its resonant frequency is ${\mathit{\omega }}_{\mathit{d}}\mathrm{=}\sqrt{{\mathit{\omega }}_{\mathrm{0}}^{\mathrm{2}}\mathrm{-}{\mathit{\lambda }}^{\mathrm{2}}}$. The periodic term of the particular solution γ is composed of the in-phase y_i and out-of-phase x_i terms ${\mathit{y}}_{\mathit{i}}\mathrm{=}{\mathit{H}}_{\mathit{i}}\frac{\mathrm{(}{\mathit{\omega }}_{\mathrm{0}}^{\mathrm{2}}\mathrm{-}{\mathit{\omega }}_{\mathit{i}}^{\mathrm{2}}\mathrm{)}{\mathit{\omega }}_{\mathrm{0}}^{\mathrm{2}}\mathrm{+}\mathrm{4}{\mathit{\lambda }}^{\mathrm{2}}{\mathit{\omega }}_{\mathit{i}}^{\mathrm{2}}}{\mathrm{(}{\mathit{\omega }}_{\mathrm{0}}^{\mathrm{2}}\mathrm{-}{\mathit{\omega }}_{\mathit{i}}^{\mathrm{2}}{\mathrm{)}}^{\mathrm{2}}\mathrm{+}\mathrm{4}{\mathit{\lambda }}^{\mathrm{2}}{\mathit{\omega }}_{\mathit{i}}^{\mathrm{2}}}\mathit{,}{\mathit{x}}_{\mathit{i}}\mathrm{=}{\mathit{H}}_{\mathit{i}}\frac{\mathrm{-}\mathrm{2}\mathit{\lambda }{\mathit{\omega }}_{\mathit{i}}^{\mathrm{3}}}{\mathrm{(}{\mathit{\omega }}_{\mathrm{0}}^{\mathrm{2}}\mathrm{-}{\mathit{\omega }}_{\mathit{i}}^{\mathrm{2}}{\mathrm{)}}^{\mathrm{2}}\mathrm{+}\mathrm{4}{\mathit{\lambda }}^{\mathrm{2}}{\mathit{\omega }}_{\mathit{i}}^{\mathrm{2}}}\mathrm{·}$(46)At first order, the expression of x_i may be simplified as ${\mathit{x}}_{\mathit{i}}\mathrm{=}\mathrm{-}\mathrm{0.9054}\frac{{\mathit{k}}_{\mathrm{2}}}{\mathit{Q}}{\mathit{H}}_{\mathit{i}}\frac{{\mathit{\omega }}_{\mathit{i}}^{\mathrm{3}}}{\mathrm{(}{\mathit{\omega }}_{\mathrm{0}}^{\mathrm{2}}\mathrm{-}{\mathit{\omega }}_{\mathit{i}}^{\mathrm{2}}{\mathrm{)}}^{\mathrm{2}}}$(47)expressed in radians and by introducing the dissipation factor as δt = (nQ)^-1. For short-period librations at 0.944898 days the x_i is 1.32 mas with $\frac{{\mathit{k}}_{\mathrm{2}}}{\mathit{Q}}\mathrm{=}{\mathrm{10}}^{-6}$ (this is the value used by Meyer & Wisdom 2008) and the resulting displacement at the surface of the satellite at the periaster passage is also negligible 0.0013 m. The damping time 1/λ is about 6000 years. If we consider $\frac{{\mathit{k}}_{\mathrm{2}}}{\mathit{Q}}$ to be 100 times bigger, i.e. closer to the expected value of Enceladus, we have a displacement at the surface of  ≈ 0.13 m and a damping time of  ≈ 60 years. For librations at long period x_i is definitely negligible because ω_i is very small.

6. Discussion

One of the aims of this theoretical study is to prepare the interpretation of potential observations of the rotation of Mimas. After a restricted analytical approach to validate the numerical results, we discuss the possibility to observe the rotation of Mimas and in particular to distinguish the different interior models. Then we focus on the non-hydrostatic contributions.

6.1. Analytical approach

We here compare with classical analytical formulae for the main term of the physical and tidal longitudinal librations and the mean obliquity, for which deriving these amplitudes accurately is quite straightforward.

6.1.1. Longitudinal librations

An analytical study of the longitudinal librations of a satellite in 1:1 spin-orbit resonance and on a Keplerian orbit can be found for instance in Murray & Dermott (1999). Let us call ψ the amplitude of the main term (i.e. associated with the mode λ + φ − ζ) of the tidal librations, and γ for the physical ones. We have from Murray & Dermott (1999)$\begin{array}{ccc}\mathit{\psi }& \mathrm{=}& \\ \mathit{\gamma }& \mathrm{=}& \end{array}$and $\begin{array}{ccc}{\left(\frac{{\mathit{\omega }}_{\mathit{u}}}{\mathit{n}}\right)}^{\mathrm{2}}& \mathrm{=}& \mathrm{3}\frac{\mathit{B}\mathrm{-}\mathit{A}}{\mathit{C}}\left(\mathrm{1}\mathrm{-}\mathrm{5}{\mathit{e}}^{\mathrm{2}}\mathrm{+}\frac{\mathrm{13}}{\mathrm{16}}{\mathit{e}}^{\mathrm{4}}\right)\\ & \mathrm{=}& \end{array}$(50)e being the eccentricity of Mimas. We can see that this amplitude is bigger when the ratio is closer to unity, or when 12C₂₂ is closer to C/(mR²). We can see from Table 2 that C₂₂ is of the order 5 × 10^-3 while C ≈ 0.4mR², i.e. C₂₂/C/(mR²) ≈ 1/80. Thus, the ratio $\mathrm{12}\frac{{\mathit{C}}_{\mathrm{22}}}{\mathit{C}\mathit{/}\mathrm{\left(}\mathit{m}{\mathit{R}}^{\mathrm{2}}\mathrm{\right)}}$ is closer to unity when C₂₂ is bigger, which is the case for the lowest values of ρ_c. Figure 3 confirms this trend, while the Table 11 settles the validity of the analytical formulae (48) and (49).

6.1.2. Mean obliquity

We here use the analytical study of Ward & Hamilton (2004; see Noyelles 2010, for an application to natural satellites in spin-orbit resonances) for the location of the Cassini states. Mimas is expected to be locked at the Cassini state 1, i.e. the most stable one, characterized by $\mathit{ϵ}\mathrm{=}\mathrm{-}\frac{\sin \mathit{I}}{\frac{\mathrm{3}\mathit{n}}{\mathrm{2}\dot{☊}}\frac{{\mathit{J}}_{\mathrm{2}}\mathrm{+}\mathrm{2}{\mathit{C}}_{\mathrm{22}}}{\mathit{C}\mathit{/}\mathrm{\left(}\mathit{m}{\mathit{R}}^{\mathrm{2}}\mathrm{\right)}}\mathrm{+}\cos \mathit{I}}\mathit{,}$(51)ϵ being the mean obliquity of Mimas, $\dot{☊}$ the precessional rate of its orbital ascending node, and I its inclination on the Laplace plane, here assumed to be the equator of Saturn at J2000.

From the definition the orbital proper modes of Mimas, we can approximate $\dot{☊}$ by rad/y, this yields a regressional period of 360.1657 days. In assuming J₂ ≈ 2 × 10^-2, C₂₂ ≈ 6 × 10^-3 and C ≈ 0.4mR² from Table 2, we have $\frac{\mathrm{3}\mathit{n}}{\mathrm{2}\dot{☊}}\frac{{\mathit{J}}_{\mathrm{2}}\mathrm{+}\mathrm{2}{\mathit{C}}_{\mathrm{22}}}{\mathit{C}\mathit{/}\mathrm{\left(}\mathit{m}{\mathit{R}}^{\mathrm{2}}\mathrm{\right)}}\mathrm{\approx }\mathrm{-}\mathrm{37.26}$ while sinI is very small and cosI close to unity (the mean inclination of Mimas I being of the order of 1.5° = 2.6 × 10^-2 rad). This means that higher values of the quantity J₂ + 2C₂₂ will yield a lower obliquity. Once more, these values are reached for the lowest estimations of ρ_c, Fig. 3 confirming this tendency. The validity of the analytical formula (51) is checked in Table 11.

The remaining question for analytical validation is: how does one evaluate the mean eccentricity and inclination required in the analytical formulae, i.e. how does one average them? These formulae were derived in assuming a Keplerian orbit, while the orbit of Mimas is perturbed by the oblateness of Saturn and the mutual interactions with the other satellites, inducing an orbital resonance with Tethys. As a consequence, its eccentricity and inclination are far from constant.

We have from Vienne & Duriez (1995)$\begin{array}{ccc}\mathit{z}\mathrm{\left(}\mathit{t}\mathrm{\right)}& \mathrm{=}& \mathit{e}\mathrm{\left(}\mathit{t}\mathrm{\right)}\exp \left(\mathrm{i}\mathit{\varpi }\mathrm{\left(}\mathit{t}\mathrm{\right)}\right)\\ & \mathrm{=}& \mathrm{1.59817}\mathrm{\times }{\mathrm{10}}^{-2}\exp \left(\mathrm{i}\left(\mathrm{6.38121472}\mathit{t}\mathrm{+}{\mathrm{356.521}}^{\mathrm{◦}}\right)\right)\\ & & \\ & & \mathrm{+}\mathrm{7.1114}\mathrm{\times }{\mathrm{10}}^{-3}\exp \left(\mathrm{i}\left(\mathrm{6.47026010}\mathit{t}\mathrm{+}{\mathrm{35.846}}^{\mathrm{◦}}\right)\right)\mathrm{+}\mathit{...}\mathit{,}\\ \mathit{\zeta }\mathrm{\left(}\mathit{t}\mathrm{\right)}& \mathrm{=}& \sin \left(\frac{\mathit{I}\mathrm{\left(}\mathit{t}\mathrm{\right)}}{\mathrm{2}}\right)\left(\mathrm{i}☊\mathrm{\left(}\mathit{t}\mathrm{\right)}\right)\\ & \mathrm{=}& \mathrm{1.18896}\mathrm{\times }{\mathrm{10}}^{-2}\exp \left(\mathrm{i}\left(\mathrm{-}\mathrm{6.37188169}\mathit{t}\mathrm{+}{\mathrm{234.213}}^{\mathrm{◦}}\right)\right)\\ & & \\ & & \mathrm{+}\mathrm{5.3017}\mathrm{\times }{\mathrm{10}}^{-3}\exp \left(\mathrm{i}\left(\mathrm{-}\mathrm{6.28283631}\mathit{t}\mathrm{+}{\mathrm{273.538}}^{\mathrm{◦}}\right)\right)\mathrm{+}\mathit{...}\mathit{,}\end{array}$the frequencies are given in rad/year, and the time origin J1980. Evidently, the mean eccentricity should be at least  ≈ 1.6 × 10^-2, probably higher (same for the mean inclination, which should be at least  ≈ 1.4°). In Table 11 we use e = 1.92 × 10^-2 and I = 1.68°, this arbitrary choice minimizes the relative errors and is consistent with the TASS1.6 theory.

6.2. Observational possibilities

It would be challenging to constrain the orientation and interior structure of Mimas using its rotation. The first expected result is the confirmation that Mimas is in the Cassini state 1 with the 1:1 spin-orbit resonance. Another challenge would be to detect the longitudinal librations, which have been actually observed for the Moon (Koziel 1967), the Martian satellite Phobos (Burns 1972), and the Saturnian satellite Epimetheus (Tiscareno et al. 2009). To estimate the required accuracy of the observations, we convert the rotation outputs into kilometres (Table 12).

As expected, the longitudinal librations are significantly bigger (a few kilometres) than the mean obliquity and the latitudinal librations (with an amplitude smaller than 200 m). The amplitude of the librations given are related to the quasi-periodic decompositions, so the peak-to-peak amplitudes are twice as big. The reader should keep in mind that the physical and tidal librations are two expressions of the same quantity, so are not independent. We can consider that the detection of the longitudinal librations would require an accuracy of about 1 km, while using them to invert the internal structure of Mimas would require an accuracy at least ten times better.

6.3. Non-hydrostatic contributions

The study of the non-hydrostatic Mimas based on the shape model does not exhibit a significant possibility to distinguish a non-hydrostatic Mimas from a hydrostatic one from observations. This is not surprising considering Mimas’ nearly hydrostatic global shape. But a non-hydrostatic Mimas could result in an offset between the ellipsoid of shape and the ellipsoid of inertia, as investigated for Janus by Robutel et al. (2011), for which an offset in longitude and in latitude has actually been detected (Tiscareno et al. 2009). Consequently, detection of non-hydrostatic contributions from observation of Mimas’ orientation should not a priori be excluded.

7. Conclusion

We have presented a theoretical study of the rotation of Mimas considering the three degrees of freedom of the rigid rotation, and different possible interior models, assuming Mimas to be in hydrostatic equilibrium, or not. Moreover, we considered a complete orbital motion, and also investigated the influence of tides on the rotation of Mimas.

We estimated the physical longitudinal librations to have an amplitude of about 0.5°, i.e. nearly 2 km, the exact value depending on the internal structure of Mimas. For a hydrostatic Mimas, a dense core lowers this amplitude. Non-hydrostatic contributions were shown to be small, as expected from Mimas’ shape in near hydrostatic equilibrium. Moreover, we expect an obliquity between 2 and 3 arcmin, while the polar motion can be neglected. The tidal deviation of Mimas’ long axis should be negligible as well, while this is the most inner main Saturnian satellite.

The Cassini spacecraft has already completed its initial four-year mission and the first extended mission, with a limited number of Mimas flybys. Its orbit close to Saturn makes Mimas a difficult target for Cassini observations. Since September 2010 Cassini is in a second extended mission called the Cassini Solstice Mission during (and especially at the end) of which Cassini will likely have additional Mimas observations. We hope that future observations of Mimas will allow us to constrain its rotation and to obtain clues on its internal structure and orientation.

Acknowledgments

Numerical simulations were made on the local computing ressources (Cluster URBM-SYSDYN) at the University of Namur. This work has been supported by EMERGENCE-UPMC grant (contract number: EME0911).

References

Appendix A: Notations used in the paper

All Tables

All Figures