{'fraction_non_alphanumeric': 0.03746085484553534, 'fraction_numerical': 0.014980956411341515, 'mean_word_length': 4.0901257970015505, 'pattern_counts': {'":': 0, '<': 17, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 9, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Tidal synchronization of an anelastic multi-layered body: Titan's synchronous rotation

This paper presents one analytical tidal theory for a viscoelastic multi-layered body with an arbitrary number of homogeneous layers. Starting with the static equilibrium figure, modified to include tide and differential rotation, and using the Newtonian creep approach, we find the dynamical equilibrium figure of the deformed body, which allows us to calculate the tidal potential and the forces acting on the tide generating body, as well as the rotation and orbital elements variations. In the particular case of the two-layer model, we study the tidal synchronization when the gravitational coupling and the friction in the interface between the layers is added. For high relaxation factors (low viscosity), the stationary solution of each layer is synchronous with the orbital mean motion (n) when the orbit is circular, but the spin rates increase if the orbital eccentricity increases. For low relaxation factors (high viscosity), as in planetary satellites, if friction remains low, each layer can be trapped in different spin-orbit resonances with frequencies n/2,n,3n/2,... . We apply the theory to Titan. The main results are: i) the rotational constraint does not allow us confirm or reject the existence of a subsurface ocean in Titan; and ii) the crust-atmosphere exchange of angular momentum can be neglected. Using the rotation estimate based on Cassini's observation, we limit the possible value of the shell relaxation factor, when a subsurface ocean is assumed, to 10^-9 Hz, which correspond to a shell's viscosity 10^18 Pa s, depending on the ocean's thickness and viscosity values. In the case in which the ocean does not exist, the maximum shell relaxation factor is one order of magnitude smaller and the corresponding minimum shell's viscosity is one order higher.


Introduction
Tidal torque is a key physical agent controlling the rotational and orbital evolution of systems with close-in bodies and may give important clues on the physical conditions in which these systems originated and evolved. The viscoelastic nature of a real body causes a non-instantaneous deformation, and the body continuously tries to recover the equilibrium figure corresponding to the varying gravitational finally, the conclusions are presented in Sect. 11. The paper is completed by several appendices where are given technical details of some of the topics presented in the forthcoming sections. In addition, an Online Supplement is provided with further details, not worthy of inclusion in the paper but useful for the reproduction of several developments.

Non-homogeneous Newtonian creep tide theory
Let us consider one differentiated body m of mass m T , disturbed by one mass point M of mass M orbiting at a distance r from the center of m. We assume that the body is composed of N homogeneous layers of densities d i (i = 1, . . . , N ) and angular velocities Ω i , perpendicular to the orbital plane.
The outer surface of the ith layer is ζ i ( ϕ i , θ i , t), where ζ i is the distance of the surface points to the center of gravity of m and the angles ϕ i , θ i are their longitudes and co-latitudes in a fixed reference system. At each instant, we assume that the static equilibrium figure of each layer under the action of the tidal potential and the rotation may be approximated by a triaxial ellipsoidal equilibrium surface ρ i ( ϕ i , θ i , t), whose semi-major axis is oriented towards M.
The adopted rheophysical approach is founded on the simple laẇ where γ i is the relaxation factor at the outer surface of the ith layer. This is a radial deformation rate gradient related to the viscosity through (see Appendix 1) where R i and g i are the equatorial mean radius and the gravity acceleration at the outer surface of the ith layer. η i is the viscosity of the inner layer (assumed to be larger than that of the outermost layers). Although the creep equation is valid in a reference system co-rotating with the body, we can use the coordinates in a fixed reference system. This is because only relative positions appear in the right-hand side of the creep equation. If ϕ F is the longitude of a point in one frame fixed in the body, then we have

The static equilibrium figure
The static equilibrium figure of one body composed by N homogeneous layers, under the action of the tidal potential and the non-synchronous rotation, when all layers rotate with the same angular velocity, was calculated by Folonier et al. (2015). 1 In this work, we need, beforehand, to extend these results to the case in which each layer has one different angular velocity. We assume that each layer has an ellipsoidal shape with outer semiaxes a i , b i and c i , where the axis a i is pointing towards M and c i is the axis of rotation. Then, the equatorial prolateness (i) ρ and polar oblateness (i) z of the outer surface of the ith layer can be written as where R i = √ a i b i is the outer equatorial mean radius of the ith layer, ρ is the flattening of the equivalent Jeans homogeneous spheroid and z is the flattening of the equivalent MacLaurin homogeneous spheroid in synchronous rotation: Here, G is the gravitation constant, R N is the equatorial mean radius of m and n is the mean motion of M. The Clairaut's coefficients H i and G i depend on the internal structure and are (see Appendix 2) where (E −1 ) ij are the elements of the inverse of the matrix E, whose elements are where x i = R i /R N and d i = d i /d 1 are the normalized mean equatorial radius and density, respectively, and f N = 3 1 0 d(z)z 2 dz. Finally, the static ellipsoidal surface equation of the outer boundary of the ith layer, to first order in the flattenings, can be written as (see Section A in the Online Supplement), where ϕ M is the longitude of M in the same fixed reference system used to define ϕ i .

The creep equation
Using the static equilibrium surface (8), the creep equation (1) becomeṡ where and v are the longitude of the pericenter and the true anomaly, respectively. For resolving the creep differential equation, we proceed in a similar way as FM13 and FM15. We consider the two-body Keplerian motion. The equations of the Keplerian motion of M are and v = + 2e − e 3 4 sin + 5e 2 4 sin 2 + 13e 2 12 where a, e, are the semi-major axis, the eccentricity, and the mean anomaly, respectively. The resulting equation is then an ordinary differential equation of first order with forced terms that may be written asζ where the arguments of the cosines Θ ik , Θ ik are linear functions of the time The constants C ik , C ik are where δ 0,k is the Kronecker delta (δ 0,k = 1 when k = 0 and δ 0,k = 0 when k = 0), the constant ρ is and E q,p are the Cayley functions (Cayley, 1861), defined by E q,p (e) = 1 2π 2π 0 a r 3 cos (qv + (p − q) ) d .
After integration, we obtain the forced terms The phases σ ik and σ ik are ; sin σ ik = kn where is the semi-diurnal frequency. These phases are introduced during the exact integration of the creep equation (12). If we define the angles and the equatorial and polar flattenings the solution (17) can be written as which has a simple geometric interpretation: using Eq. (8), we can identify each term of the Fourier expansion of the height δζ i , with the boundary height of one ellipsoid, with equatorial and polar flattenings (ik) ρ and (ik) z , respectively, rotated at an angle δ ik /2, with respect to the axis x.

The disturbing potential
The potential of the ith layer of m at a generic point M * (r * , θ * , ϕ * ) external to this layer, can be written as the potential of one spherical shell of outer and inner radii R i and R i−1 , respectively, plus the disturbing potential due to the mass excess or deficit corresponding to the outer and the inner boundary heights δζ i and δζ i−1 . It is important to note that since these excesses or deficits are very small, we may calculate the contribution of each term of the Fourier expansion separately and then sum them to obtain the total contribution. In this way, we assume that the ith layer has outer and inner boundary heights given by the kth term of the Fourier expansion. The equatorial and polar flattenings of the outer boundary, (ik) ρ and (ik) z , are given by Eq. (21), and the bulge is rotated at an angle δ ik /2 with respect to the axis x. Similarly, the inner boundary height δζ (1) i−1 , can be identified with the boundary height of one ellipsoidal surface, with equatorial and polar flattenings rotated at an angle δ i−1,k /2, with respect to the axis x.
The disturbing potential at an external point M * (r * , θ * , ϕ * ), due to the mass excess or deficit, corresponding to the kth term of the Fourier expansion of the outer and the inner boundary heights δζ i and δζ i−1 , is where C i is the axial moment of inertia of the ith layer (see Section A in the Online Supplement) and ∆(f i ) = f i − f i−1 , denotes the increment of one function f i , between the inner and the outer boundaries of this layer. Taking into account that the total disturbing potential of the ith layer, can be approximated by the sum of the contribution of each term of the Fourier expansion, we obtain

Forces and torques
To calculate the force and torque due to the ith layer of m, acting on one mass M * located in M * (r * , θ * , ϕ * ), we take the negative gradient of the potential of the ith layer at the point M * and multiply it by the mass placed in the point, that is, The corresponding torque is M i = r * × F i , or since r * = (r * , 0, 0), that is 5 Forces and torques acting on M Since we are interested in the force acting on M due to the tidal deformation of the ith layer of m, we must substitute (M * , r * , θ * , ϕ * ) by (M, r, π 2 , + v). Replacing the angles δ ik and δ ik given their definitions (Eq. 20), the forces, then are and the corresponding torques are After Fourier expansion, the torque along to the axis z (M zi = −M 2i ), can be written as Finally, replacing the coefficient C ik given by Eq. (14), the time average of the tidal torque over The above expression for the time average, which is equivalent to take into account only the terms with j = 0, is only valid if ν i is constant. This condition is satisfied, for example, by homogeneous bodies with γ n, as stars and giant gaseous planets, where the stationary rotation is ∼ 6nγe 2 /(n 2 + γ 2 ). However, the final rotation of the homogeneous rocky bodies, with γ n, as satellites and Earth-like planets, is dominated by a forced libration ∼ B 1 cos ( + φ 1 ) with the same period as the orbital motion of the system (see Chap. 3 of FM15). In this case, any time average that involves the rotation, should also take into account this oscillation. It is worth emphasizing that in this paper we calculate the time average of some quantities, as the work done by the tidal forces and the variations in semi-major axis and eccentricity, assuming that ν i is constant, which is valid only for bodies with low viscosity. The applications to Titan in this paper were done using the complete equations, where the distinction between these extreme cases is not necessary.
6 Work done by the tidal forces acting on M The time rate of the work done by the tidal forces due to the ith layer isẆ i = F i · v, where v is the relative velocity vector of the external body 2 whose components in spherical coordinates are Using the tidal force, given by the Eq. (29), the rate of the work corresponding to the ith layer is or after Fourier expansion 3 The time-average over one period is The average of the term involving δ 0,k G i z in the last term of Eq. (35), for k = 0, is (see Section C in the Online Supplement). 2 The definition of power (the time derivative of work) used in this section is the most general definition of the power done by the force couple formed by the disturbing force Fi acting on the external body and its reaction −Fi acting on the ith layer of the deformed body. It may be written as Fi · V M + (−Fi) · Vi where V M and Vi are, respectively, the velocities of the body M and of the ith layer of the body m, w.r.t. a fixed reference frame. It is equivalent to Fi(V M − Vi), that is, Fi · v (see Scheeres, 2002;Ferraz-Mello et al., 2003). 3 For the details of the calculation see the Online Supplement of FM15.

Variations in semi-major axis and eccentricity
In this section, we calculate the variation in semi-major axis and eccentricity. As in FM13 and FM15, we use the energy and angular momentum definitions. 4 If we differentiate the equation where a is the semi-major axis of the relative orbit, we obtain the equation for the variation in semimajor axis:ȧ ReplacingẆ by the Eq. (35) and summing over all layers, we obtaiṅ After the averaging over one period, we obtain In the same way, if we differentiate the total angular momentum equation where e is the eccentricity of the relative orbit, and useṅ/n = −3ȧ/2a, we obtain the equation for the variation in eccentricity eė whereL = M z is the total torque exerted by the tidal forces. The interaction torques between the layers do not affect the orbital motion, because they are action-reaction pairs (that is M ij = −M ji , ∀ i, j = 1, . . . , N and i = j), then they mutually cancel themselves. ReplacingẆ and M z by the Eqs. (31) and (35), and summing over all layers, we obtaiṅ After the time-average over one period, we obtain that the variation in eccentricity are 8 The two-layer model In the previous sections we have studied the tidal effect on one body composed of N homogeneous layers. However, in contrast with a homogeneous body, in one differentiated body we must also take into account the interaction between the different layers. In this paper we consider the gravitational coupling of the layers and the friction that occurs at each interface of two layers in contact. An important point to keep of mind is that the number of free parameters increases significantly as the number of layers increases.
In this section, we study the simplest non-homogeneous problem: one body formed by two independent rotating parts. The inner layer, or core, is denoted with the subscript c and the outer layer, or shell, is denoted with the subscript s. Despite its simplicity, the two-layer model allows us to study the main features of the stationary rotations, introducing a minimum number of free parameters.

The tidal torques
The tidal torques due to the core and the shell, along the axis z, are (see Eq. 31) where the function T i (i = c, s) is the constants T ij are and the tidal parameter T , is defined as R c , C c are the mean outer radius and moment of inertia of the core, and R s , C s are the mean outer radius and moment of inertia of the shell. The parameters H c , γ c are the Clairaut parameter and the relaxation factor at the core-shell interface and H s , γ s are the Clairaut parameter and the relaxation factor at the body's surface.

The gravitational core-shell coupling
When the principal axes of inertia of two layers are not aligned, a restoring torque appears which tends to align these axes again. This torque was calculated by several authors (e.g. Buffett, 1996 Here, we use one similar expression for this torque adapted to a body formed by two layers whose boundaries are prolate ellipsoids, whose flattenings are defined by the composition of the main elastic and anelastic tidal components. If we follow the same composition adopted in FM13, these flattenings are where 0 < λ c < 1 and 0 < λ s < 1 are the relative measurements of the actual maximum heights of the elastic tides of the core and the shell, respectively. The geodetic lags of these two ellipsoidal surfaces, when one elastic component is added are In this case, the torques, along the axis z, are where ξ = ϑ s − ϑ c is the offset of the geodetic lags of the two ellipsoidal boundaries and the constant of gravitational coupling K is (see Appendix 3 for more details). We may pay attention to the sign of these torques. If ϑ s > ϑ c , the motion of the shell is braked, while the motion of the core is accelerated. This is consistent with the signs of the above equations.

Linear drag
The model considered here also assumes that a linear friction occurs between the two contiguous layers. For the two-layer model, the torques acting on the core and the shell, along the axis z, are where the friction coefficient µ is an undetermined ad-hoc constant that comes from assuming that a linear friction occurs between two contiguous layers. When we consider that the body m has solid layers, but not rigid, we can assume that between the core and the shell exists one thin fluid boundary with viscosity η o and thickness h. If this interface is a Newtonian fluid, the Eq. (52) is the law corresponding to liquid-solid boundary for low speeds, and µ can be written as (see Appendix 4 for more details).

Rotational equations
Putting together all contributions to the torque, we obtain the rotational equations where M core z and M shell z are the z-components of the total torque acting on the core and on the shell. These torques include the reaction of the tidal torque M zi , the gravitational coupling Γ i and the friction Φ i .

Comparison with the homogeneous case
In this section, we compare some of the main features of the homogeneous creep tide theory, developed in FM15, with the non-homogeneous creep tide theory for the two-layer model developed in this article. The main difficulty lies in the number of free parameters in these approaches. In the homogeneous case, with a suitable choice of dimensionless variables, the final state of rotation depends only on the ratio n/γ and on the eccentricity e (Eq. 42 of FM15). However, even in the most simple non-homogeneous case (the two-layer model), we need to set 12 free parameters. In order to proceed, we use the typical values for Titan and also Titan's eccentricity e = 0.028 (see Tables 1-4 in Sect. 10.1), and let as free parameters, only n/γ i , e and µ.
Following FM15, we introduce the adimensional variables y i = ν i /γ and the scaled time x = /γ, where γ = 2γ c γ s (γ c + γ s ) −1 . If we consider the case in which γ c = γ s , the behavior of the evolutions of y c and y s is similar to that observed in the homogeneous case. Figure 1 shows the time evolution of y s , with initial conditions y c = 0.3, y s = 0.15 and differents values of α = log 10 (n/γ c ) = log 10 (n/γ s ). When γ i n (i.e. rocky bodies), after a transient, the solution oscillates around zero, independently of the initial conditions (left panel), and the amplitude of oscillation decreases when α decreases. In the case where α = 4, we also plot the solution with initial conditions y c = 0.3 and y s = −0.15 (dashed black line). This solution increases quickly, becoming indistinguishable from the solution with initial value y s = 0.15. When γ i ∼ n, the stationary solution becomes a super-synchronous rotation with the amplitude of oscillation tending to zero, and, finally, when γ i n, the stationary solution of y s becomes closer zero (right panel). The evolution of y c is very similar, and the friction does not have any relevant role. However, when we analyze the time evolution of ν s instead of y s , we observe that when γ i n, the solution oscillates around zero and the amplitude of the oscillation of ν s increases when α decreases (left panel of Fig. 2). When γ i n. this amplitude decreases when α decreases and ν s tends to 12ne 2 , independently of the value of α (right panel of Fig. 2).
When γ c = γ s , we can have different core and shell rotation behavior. In Fig. 3, we show the core and shell rotation (left and right, respectively) for log 10 (n/γ c ) = 2 and log 10 (n/γ s ) = 4. We also set two very different values for the friction: the frictionless case µ = 0 (black) and a very high value of friction µ = 10 28 kg km 2 s −1 (red lines), larger than the expected value in the case of Titan (µ = 10 11 − 10 13 kg km 2 s −1 ), which corresponds to a typical ocean viscosity η o = η H2O ≈ 10 −3 Pa s and a large range for the ocean thickness h (see Eq. 53). In the frictionless case, we can observe the differential rotation between the core and the shell. After a transient, both solutions oscillate around zero with very different amplitudes, depending on the value of γ of each surface. For very high friction parameter, both layers rotate with the same angular velocity. The core and the shell have the same amplitude of oscillation and phase, keeping the relative velocity equal to zero.  Finally, we study the dependence of the stationary solutions on the eccentricity. For that sake, we choose a grid of initial conditions ν c /n and ν s /n, and integrate the system (54) until the stationary solution is reached. When n/γ c = n/γ s 1, all initial conditions lead to the same equilibrium point (a super-synchronous rotation), independently of the value of the friction parameter. The value of the excess of rotation depends only on the eccentricity. In the left panels of Fig. 4, we show the family of stationary solutions, where each point corresponds to a different eccentricity value in 0 ≤ e ≤ 0.5. If the eccentricity is zero, the rotations are synchronous to the orbital motion. When the eccentricity increases, the rotations become super-synchronous, and the excess of rotation ν i /n is proportional to e 2 (right panels).
When n/γ c and n/γ s increase, that is, when the viscosities increase, the excess in the supersynchronous rotation decreases. If the eccentricity is low, the only attractor is the super-synchronous solution. When the eccentricity increases, captures in other attractors ν i n, 2n, 3n, . . . appear gradually. This behavior is the same studied by in FM15 and also in Correia et al. (2014) in the case of homogeneous bodies. Figure 5 shows the families of stationary rotation for n/γ c = n/γ s = 1, 0 ≤ e ≤ 0.5 and two values of the friction parameter: the frictionless case, with µ = 0 (top panels), and a very high friction case, with µ = 10 20 kg km 2 s −1 (bottom panels). In the frictionless case, when the eccentricity is smaller than ∼ 0.48, only the super-synchronous solution is possible. If the eccentricity is larger than 0.48, besides the super-synchronous solution, three new stationary configurations appear: The core and the shell in the 3/2 commensurability (ν c n and ν s n), the core in super-synchronous rotation and the shell in the 3/2 commensurability (ν c 0 and ν s n), and the core in the 3/2 commensurability and the shell in super-synchronous rotation (ν c n and ν s 0). Figure 6 shows in more detail these stationary solutions. The labels R pq denote the stationary families indicating the resonances ν c = pn and ν s = qn. It is important to note that the excesses in the rotations are large because the eccentricity is high. In the high friction case (bottom panels of Fig. 5), only the stationary solutions with the same commensurabilities survive because in these configurations, the relative velocity of rotation between the core and the shell is zero. If n/γ c and n/γ s continue to increase and the friction parameter is low (not necessarily zero), the core and the shell may tend to different resonances, depending on the eccentricity. If the friction increases, the attractors with higher differential rotation, begin to disappear, until eventually, as from a certain value limit of µ only survive the attractors with differential rotation zero Fig. 7.  Table 1. In this section, we assume the interior model given by Sohl et al. (2014) (hereafter reference model ), which is given in Table 2. In this model, Titan is formed by four homogeneous layers: i) an inner hydrated silicate core (inner core); ii) a layer of high-pressure ice (outer core); iii) a subsurface water-ammonia ocean and iv) a thin ice crust. For the sake of simplicity, we construct one two-layer equivalent model, where the core is a layer formed by the inner core and the high-pressure ice layer, and the shell is a layer formed by the subsurface ocean and the ice crust, but keeping some features of the four-layer model (e.g. axial moments of inertia and Clairaut numbers). In this way, we can use the rotational equations (54), retaining the main features of the realistic reference model. This simplified model is given in Table 3, and some calculated parameters of each layer are listed in Table 4. The existence of relative translational motions due to the non-coincidence of the barycenters of the several layers, as discussed by Escapa and Fukushima (2011) in the case of an icy body with an internal ocean and solid constituents, has not been taken into account.
In order to estimate the relative height of the elastic tide λ s , we assume that the difference between the observed surface flattening s with the tidal flattening s = H s ρ E 2,0 cos σ s0 ≈ H s ρ cos σ s0 (calculated) is due to the existence of an elastic component, with flattening (el) s = λ s H s ρ (see Appendix 3 for more details). If we use Eq. (48), and assume that near the synchronous rotation cos 2 σ s0 ≈ 1, we obtain For the relative heights of the elastic tide λ c , we assume λ c ≈ λ s def = λ.

Atmospheric influence on Titan's rotation
The seasonal variation in the mean and zonal wind speed and direction in Titan's lower troposphere causes the exchange of a substantial amount of angular momentum between the surface and the atmosphere. The variation calculated from the observed zonal wind speeds shows that the atmosphere angular momentum undergoes a periodic oscillation between 3 × 10 18 and 3 × 10 19 kg km 2 s −1 (Tokano and Neubauer, 2005, hereafter TN05) with a period equal to half Saturn's orbital period and maxima at Titan's equinoxes (when the Sun is in the satellite's equatorial plane). The angular momentum of the atmosphere may be written as L atm = L 0 + L 1 cos 2α where L 0 = 1.65 × 10 19 kg km 2 s −1 , L 1 = 1.35 × 10 19 kg km 2 s −1 and α is the Saturnian right ascension    of the Sun. The variation of the angular momentum isL atm = −2L 1 n sin 2α . If we neglect external effects (as atmospheric tides), this variation may be compensated by an equal variation in the shell's angular momentum: δL s = −L atm , which corresponds to an additional shell acceleration We must emphasize that we have considered in these calculations the moment of inertia of the ice crust C k , since the winds are acting on the crust and do not have direct action on the liquid part of the shell.
In a more recent work, Richard et al. (2014) (hereafter R14) re-calculate the amplitude of the variation of the angular momentum with the Titan IPSL GCM (Institut Pierre-Simon Laplace General Circulation models) (Lebonnois et al., 2012). They obtain L 1 = 8.20×10 17 kg km 2 s −1 , which is ∼ 16.5 times less than the TN05 value.

The results
We fix the outer radius of the inner core R ic and the outer radius of the high-pressure ice layer R oc , the densities of the inner and outer cores d ic and d oc and the density of the crust d k , to the reference  (Tokano and Neubauer, 2005). model values in Table 2. The density of the inner core is calculated so as to verify the value of Titan's mass m T = 13.45 × 10 22 kg. Figure 8 shows the weak dependence of the parameters on the thickness The main consequence of the weak dependence of these parameters with the thickness of the subsurface ocean, is that both the effect of the tide and the gravitational coupling parameter also depend weakly on h. The strength of the acceleration of the rotation, due to the tide, is given by the product T ij T k (see Eqs. 45 and 46). While the parameter T ij only depends on the internal structure of Titan, the function T k do not depend on h. The left panel of Fig. 9 shows T ij and the gravitational coupling amplitude K i = K/C i , as function of h. We also observe that the thickness of the ocean does not have any relevant role. Then, for the tide and the gravitational coupling, the rotational evolution is driven by the ratios n/γ c , n/γ s and the orbital eccentricity e.
The right panel of Fig. 9 shows the quantity nµ i = nµ/C i as function of the thickness h, when we consider the realistic ocean viscosity η o = η H2O ≈ 10 −3 Pa s. The rotational acceleration of each layer, due to the friction, is µ i (Ω s − Ω c ). In super-synchronous rotation, the excess of rotation of each layer is of order ne 2 , then Therefore, in Titan's case, the friction term is negligible compared with the tide and the gravitational coupling terms, independently of the h value. Equations (54) and (56), allow us to calculate the velocities of rotation of the shell and of the core of Titan for a wide range of relaxation factors γ c and γ s , when different effects are considered. For that sake, we have to adopt the values of the involved parameters. We use four different values for the viscosity of the subsurface ocean: a realistic value η o = η H2O = 10 −3 Pa s, a moderate value η o = 10 0  The periodic attractor of the spin rate ν i of each layer can be approximated by the trigonometric polynomial where the constants B ij and the phases φ ij , depend on the relaxation factors. The tidal drift B i0 also depends on e 2 , while the amplitude of oscillation B ij , depends on e j . The coefficients B ij and φ ij gives rise to intricate analytical expressions, but are easy to calculate numerically (an analytical construction of these constants is presented in the Section B of the Online Supplement). Figure 11 shows one example for Titan's core and shell constants B cj and B sj , as a function of the shell relaxation factor, when the core relaxation factor is γ c = 10 −8 s −1 , and the ocean's viscosity and thickness are η o = 10 −3 Pa s and h = 178 km, respectively. We can observe that if γ s 10 −7.5 s −1 , the shell oscillates around the super-synchronous rotation. When γ s 10 −7.5 s −1 , the tidal drift B s0 tends to zero and the shell oscillates around the synchronous rotation, with a period of oscillation equal to the orbital period. Finally, if γ s 10 −8 s −1 , the amplitude of the shell rotation decreases, tending to zero when γ s decreases. On the other hand, the core oscillates around the synchronous rotation, with a period of oscillation equal to the orbital period, independently of the shell relaxation factor. In Fig. 12, fixing η o = 10 −3 Pa s and L 1 = 1.35 × 10 19 kg km 2 s −1 (TN05), we plot the resulting maximum and minimum of the final oscillation of the shell rotation Ω s − n, or, equivalently, the length-of-day variation  Figure 12 shows that if γ s < 10 −7 s −1 , the shell's rotation oscillate around the synchronous motion and the amplitude of oscillation depends on the relaxation factors and the ocean thickness. The average rotation (central orange line) is synchronous; it only becomes super-synchronous for relaxation values larger than ∼ 10 −6.5 s −1 . We also observe that when γ s < 10 −8 s −1 , independently of the values of γ c and h, the amplitude of oscillation of the shell tends to zero when the relaxation factor γ s decreases. Particularly, if γ s < 10 −9 s −1 , the amplitude of the oscillation of the excess of rotation reproduces the dispersion of the Ω s value of ±0.02 deg/yr around the synchronous value, observed as reported by Meriggiola (2012) and Meriggiolla et al. (2016). The results are not consistent with the previous drift reported by . We note that for larger values of the relaxation, e.g. 10 −8 s −1 , the large short period oscillation due to the tide would be much larger than the reported values and would introduce big dispersion in the measurements, much larger than the reported dispersion due to the difficulties in the precise localization of Titan's features. On the other hand, the effect of the atmospheric torque is completely negligible in the range of possible γ s that reproduces the observed values of the shell rotation, even for the high value of L 1 given Tokano and Neubauer (2005). When we consider the amplitude of the variation of the angular momentum given by Richard et al. (2014), the contribution to the rotation variations tends to zero. The results shown in Fig. 12 remain virtually unchanged when the ocean viscosity is increased up to a value of η o = 10 6 Pa s. But if the ocean viscosity is increased to η o = 10 9 Pa s, the transfer of angular momentum between the shell and the core induces in the shell accelerations of the same order as the rotational acceleration due to the others forces. As a consequence, the shell rotation will follow the core rotation closely (which is shown in Fig. 13). This high value of η o can be interpreted as the ocean thickness tending to zero. In this case, to obtain the dispersion of Titan's observed rotation as determined by  we should have a value of γ s smaller than the values obtained in the previous cases, where a low viscosity ocean was assumed between the shell and the core. It is worth noting yet that, in this case, the observed dispersion could also be obtained taking for γ c an extremely low value (10 −9 s −1 ) and for γ s a much larger and unexpected value (10 −5 s −1 ). It is important to note that, in any case, the rotational constraint does not allow us to estimate the value of the core relaxation factor γ c . For realistic values of the ocean viscosity (η o = 10 −3 − 10 6 Pa s), the shell relaxation factor may be such that γ s 10 −9 s −1 . The actual value will depend on the values of h and γ c and on the interpretation of the dispersion determined by Meriggiola, which may include the forced short-period oscillation of Ω s . Equivalently, using Eq. (2), the shell viscosity may be such that η s 10 18 Pa s. These values remain without significant changes if η o < 10 9 Pa s. For the case in which a subsurface ocean does not exist, the shell relaxation factor may be such that γ s 10 −10 s −1 , one order less than when an ocean is considered. Equivalently, the shell viscosity may be such that η s 10 19 Pa s. It is worth noting that in this case, when γ s 10 −7 s −1 , the rotation of the core remains stuck to the rotation of the shell even when γ c is larger, notwithstanding the larger moment of inertia of the core (Fig. 14).

Conclusion
In this article we extended the static equilibrium figure of a multi-layered body, presented in Folonier et al. (2015), to the viscous case, adapting it to allow the differential rotation of the layers. For this sake, we used the Newtonian creep tide theory, presented in Ferraz-Mello (2013) and Ferraz-Mello (2015a). Once solved the creep equations for the outer surface of each layer, we obtained the tidal equilibrium figure, and thereby we calculated the potential and the forces that act on the external mass producing the tide.
In order to apply the theory to satellites of our Solar System, we calculated the explicit expression in the particular case of one body formed by two layers. We may remember that the number of free parameters and independent variables increases quickly when the number of layers increases. The simplest version of the non-homogeneous creep tide theory (the two-layer model), allows us to obtain the main features due to the non-homogenity of the body, by introducing a minimal quantity of free parameters. In the used model, we have also calculated the tidal torque, which acts on each layer and also the possible interaction torques, such as the gravitational coupling and the friction at the interface between the contiguous layers (general development of these effects are given in Appendices 3 and 4). The friction was modeled assuming two homogeneous contiguous layers separated by one thin Newtonian fluid layer. This model of friction is particularly appropriate for differentiated satellites with one subsurface ocean, as are various satellites of our Solar System (e.g. Titan, Enceladus, and Europa).
The two-layer case was compared with the homogeneous case. For that sake, we fixed the free parameters of Titan and studied the main features of the stationary solution of this model in function of a few parameters, such as the relaxation factors γ i , the friction parameter µ and the eccentricity e. When γ c ≈ γ s , the behavior of the stationary rotations turned out to be identical to the homogeneous case. When γ c ≈ γ s n, the stationary solutions oscillate around the synchronous rotation. When γ c and γ s increase, the oscillation tends to zero. Finally, if γ c ≈ γ s n, the stationary solution is damped to super-synchronous rotation. We have also calculated the possible attractors when the eccentricity and the friction parameter µ are varied. We recovered the resonances trapping in commensurabilities Ω c ≈ Ω s ≈ 2+k 2 n (where k = 1, 2, 3, . . . ∈ N) as shown in Ferraz-Mello (2015a) and Correia et al. (2014) for the homogeneous case, and we found that if friction remains low, the non-zero differential rotation commensurabilities Ω c ∼ 2+i 2 n and Ω s ∼ 2+j 2 n, with i, j = 1, 2, 3, . . . ∈ N and i = j, are possible. When the friction increases, the resonances with higher differential rotation are destroyed. If µ continues increasing, only the resonances in which core and shell have the same rotation survive. This behavior is also observed in the non-homogeneous Darwin theory extension, when one particular ad hoc geodetic lag and one dynamical Love number for each layer are chosen (Folonier, 2016).
The two-layer model was applied to Titan, but adding to it the torques due to the exchange of angular momentum between the surface and the atmosphere, as modeled by Tokano and Neubauer The main result was that the rotational constraint does not allow us to confirm or reject the existence of a subsurface ocean on Titan. Only the maximum shell's relaxation factor γ s can be determined, or equivalently, the minimum shell's viscosity η s , because the icy crust is rotationally decoupled from the Titan's interior. When a subsurface ocean is considered, the maximum shell's relaxation factor is such that γ s 10 −9 s −1 , depending on the ocean's thickness and viscosity values considered. Equivalently, this maximum value of γ s , corresponds with a minimum shell's viscosity η s 10 18 Pa s, some orders of magnitude higher than the modeled by Mitri et al. (2014). When the non-ocean case is considered, the maximum shell's relaxation factor is such that γ s 10 −10 s −1 and the corresponding minimum shell's viscosity is η s 10 19 Pa s. For these values of γ s , the amplitude of the oscillation of the excess of rotation reproduces the dispersion of the Ω s value of ±0.02 deg/yr around the synchronous value, observed as reported by Meriggiola (2012) and Meriggiolla et al. (2016). It is important to note that in all the cases studied, the influence of the atmosphere can be neglected, since it does not affect the results in the ranges of γ c and γ s where the excess of rotation calculated is compatible with the excess of rotation observed.

Appendix 1: Relaxation factor
Let us consider the equilibrium surface ρ i (φ, θ) between two adjacent homogeneous layers of the body m whose densities are d i (inner) and d i+1 (outer). We consider that at a given instant, the actual surface between the two layers ζ i (φ, θ) does not coincide with the equilibrium surface (Fig. 15). In some parts, the separation surface is above the equilibrium surface (as in region I) and in other parts it is below the equilibrium surface (as in region II). Let us now consider one small element of the equilibrium surface in region I. The pressure in the base of this element is positive because the weight of the column above the element is larger than its weight in the equilibrium configuration. Note that the column is now partly occupied by the fluid with density d i and d i > d i+1 . The pressure surplus is given by where ∆w = (d i −d i+1 )g is the difference of the specific weight of the two columns in the neighborhood of the separation surface, and h is the distance of the element of the equilibrium surface to the actual separation surface. g is the local acceleration of gravity. Fig. 15 Interface between two adjacent homogeneous layers of m whose densities are di (inner) and di+1 (outer). ζi(φ, θ) and ρi(φ, θ) are the actual and the equilibrium surfaces, respectively, of the outer boundary of the ith layer. Region I (resp. II) is where the actual surface is above (resp. below) the equilibrium surface. FI (resp. FII ) is the force acting on one small element of the equilibrium surface in the region I (resp. region II) due to the pressure surplus (resp. pressure deficit).
The radial flow in the considered element is ruled by the Navier-Stokes equation: where F ext is the external force per unit volume, u is the radial velocity and η i is the viscosity of the layer i (assuming η i > η i+1 ). We notice that ∆ is operating on a vector, contrary to the usual ∆. Actually, in this pseudo-vectorial notation, the formula refers to the components of u and means the vector formed by the operation of the classical ∆ on the three components of the vector u. We assume that the flow is orthogonal to the equilibrium surface. We remind that, by the definition of the equilibrium surface, the tangential component of the resultant forces 5 acting on the fluid vanishes at the equilibrium surface. Since the equilibrium surface is an almost spherical ellipsoid, we may consider in a first approximation that the motion of the fluid in that region is a radial flow.
If we consider that F ext = 0 (no other external forces are acting on the fluid) and restricting u to its radial component u r , there follows Hence, The general solution of this equation is where C 1 and C 2 are integration constants. The task of interpreting and determining its integration constants becomes easier if the solution is linearized in the neighborhood of r = ρ i (i.e. h = 0): Hence, u r (ρ i ) = 0, that is, there is no pressure surplus (or deficit) when the actual separation surface coincides with the equilibrium and the linear approximation of the solution is obtained when we assume u r (ρ i ) = 0. Therefore, Hence, u r (ρ i ) = ρ i ∆w/2η i , and the linear approximation corresponding to the Newtonian creep of the fluid is where In the region II, the calculation is similar; however, instead of a pressure surplus we have a pressure deficit because the equilibrium assumes one fluid with density d i below the equilibrium surface, which is now occupied by fluid of density d i+1 < d i . The equations are the same as above. We note that in the new equations, the adopted viscosity continues being η i since we assumed it larger than η i+1 . The relaxation of the surface to the equilibrium is governed by the larger of the viscosities of the two layers.
In the homogeneous case we have one layer body (N = 1). If we consider d N +1 = 0 (neglected the density of the atmosphere), we recover the expression of the relaxation factor given by Ferraz-Mello (2013; 2015a) where w = d N g is the specific weight and ρ N ≈ R N .

Appendix 2: Equilibrium ellipsoidal figures
In this appendix we calculate the equatorial and the polar flattenings of the equilibrium ellipsoidal figures for differentiated non-homogeneous bodies in non-synchronous rotation when each layer has a different angular velociy. For this, we extend the results obtained by Folonier et al. (2015), where the rigid rotation hypothesis has been assumed. Let us consider one body m of mass m T and one mass point M of mass M orbiting at a distance r from the center of m. We assume that the body is composed of N homogeneous layers of density d i (i = 1, . . . , N ) and angular velocity Ω i = Ω ik , perpendicular to the orbital plane (Fig. 16). We also assume that each layer has an ellipsoidal shape with outer semiaxes a i , b i and c i ; the axis a i is pointing towards M while c i is along the axis of rotation.
The equatorial prolateness and polar oblateness of the ith ellipsoidal surface, respectively, are where R i = √ a i b i is the outer equatorial mean radius of the ith layer.
where G is the gravitation constant and If we assume that where M is the flattening of the equivalent MacLaurin homogeneous spheroid in synchronous rotation and J is the flattening of the equivalent Jeans homogeneous spheroids: the Eq. (70) can be written as where x i = R i /R N is the normalized mean equatorial radius and the coefficients α ij , β ij and γ i are Then, the Clairaut's coefficients H i and G i are where (E −1 ) ij are the elements of the inverse of the matrix E, whose elements are where d i = d i /d 1 is the normalized density of the ith layer and f N = 3 1 0 d(z)z 2 dz.

Appendix 3: Gravitational coupling
When the principal axes of inertia of two layers (of one body composed by N homogeneous layers) are not aligned, a restoring gravitational torque, which tends to align these axes appears. The torque acting on the inner jth layer due to the outer ith layer (not necessarily contiguous) is where d j , m j are the density and the mass in the jth layer and δU i is the disturbing potential of the ith layer at an external point (Fig. 17). The limits of the integral in Eq. (78), ζ j and ζ j−1 , are the real outer and inner boundaries of the jth layer, respectively. In our model we have to consider the actual flattening of the surfaces, which is the composition of the main elastic and anelastic tidal components (see Sec. 10 of Ferraz-Mello, 2013). The addition of the two components is virtually equivalent to the use the Maxwell viscoelastic model ab initio as done by Correia et al. (2014) (Ferraz-Mello, 2015b).
Assuming that the elastic and the anelastic components have ellipsoidal surfaces (not aligned), the resulting surface can be approximated by a prolate ellipsoid with equatorial flattening and rotated by an angle ϑ with respect to M. For the sake of simplicity, we also assume that the relative motion of the outer body M is circular. Then, neglecting the axial term does not contribute to the calculation of the gravitational coupling, the height of the outer surface of the jth layer with respect to the one sphere of radius R j , in polar coordinates, rotated by an angle ϑ j with respect to M and to first order in the flattenings (see Fig. 18), is where 0 < λ j < 1 is a relative measurement of the maximum height of the elastic tides of the outer boundary of the jth layer. The angle ϑ j is often called the geodetic lag of the surface.  If we open the trigonometric functions, by identification of the terms with same trigonometric arguments, the resulting equatorial flattening of the outer boundary of the jth layer is and the geodetic lag is The height of the inner boundary of the jth layer, taking into account the composition of the main elastic and anelastic tides has an identical expression: where 0 < λ j−1 < 1 is the relative measurement of the maximum height of the elastic tides of the inner boundary of the jth layer. Then, the resulting equatorial flattening is and the geodetic lag is In the same way, we assume that the ellipsoidal shape of this layer is also given by the composition of the main elastic and anelastic tidal components. Then, the inner and outer equatorial flattenings, respectively, are and the corresponding geodetic lags are where 0 < λ i , λ i−1 < 1 are the relative measurements of the maximum heights of the elastic tides of the outer and inner boundaries of the ith layer. Using the expression of the disturbing potential, given in Section A in the Online Supplement, and neglecting the axial term, we obtain where ∆(f i ) = f i − f i−1 , denotes the increment of one function f i between the inner and the outer boundaries of this layer. Then, the vectorial product in Eq. (78) is Using the polar unitary vectors in Cartesian coordinates θ = cos θ cos ϕ x + cos θ sin ϕ y − sin θ z ϕ = − sin ϕ x + cos ϕ y, and the approximation of ln ζ j /ζ j−1 to first order in the flattenings then, we may perform the integrals of Eq. (78) and obtain the torque acting on the inner jth layer due to the outer ith layer where As the torque acting on the outer ith layer, due to the inner jth layer, is the reaction then, the total gravitational coupling, acting on the jth layer can be written as If we consider the two-layer model, the torque acting on the core and the shell, are where the gravitational coupling parameter K is The equatorial flattenings are and the geodetic lags are The parameters 0 < λ c , λ s < 1 are relative measurements of the heights of the elastic tides of the outer surfaces of the core and the shell, respectively. The trigonometric functions in (96)-(97) are frequency functions. Using Eq. (18), an elementary calculation shows that It is important to note that some works, as Karatekin et al. (2008), use a different gravitational coupling parameter K. When applied to a two-layer model, their K differs from Eq. (95) by a multiplicative factor (1 − d s /d c ). The reason for this difference is simple: while here we calculate the torque due to the mutual gravitational attraction of two layers, through Eq. (91), they calculate the gravitational coupling between regions that involve various layers simultaneously (see Fig. 2 and Eq. 16 of Van Hoolst et al., 2008).
The element of area is R i dθ × R i sin θdφ. The integral of dΦ i,i+1 over the sphere is easy to calculate giving If we compare with the law used to introduce the friction, we obtain This is the law corresponding to a liquid-solid boundary for low speeds. The torque, along the axis z, acting on the inner (i + 1)th layer due to the outer ith layer is Then, the total torque, due to the friction, acting on the ith layer is the sum of the torque due to the outer (i + 1)th layer plus the the torque due to the inner (i − 1)th layer In the two-layer model, the torque acting on the core due to the shell and the torque acting on the shell due to the core are, respectively where η o and h are the viscosity and the thickness, respectively, of the core-shell boundary and Online Supplement A Shape and gravitational potential of one ellipsoid and one ellipsoidal layer

A.1 Homogeneous ellipsoid
Let us consider a homogeneous triaxial ellipsoid with density d, semi axes a > b > c, equatorial mean radius R = √ ab and equatorial and polar flattenings are Then, the semi axes of this ellipsoid, to first order in the flattenings, can be written as Let us consider the equation of surface of this homogeneous triaxial ellipsoid, in a reference system where the semi axes a and c are aligned to the coordinates axes x and z, respectively: If we use the semi axes (A.2), the spherical coordinates x = ρ sin θ cos ϕ; y = ρ sin θ sin ϕ; z = ρ cos θ, (A.4) and expand to first order in the flattenings, we obtain The mass of this ellipsoids is The principal moments of inertia are and its differences are The corresponding gravitational potential of this homogeneous triaxial ellipsoid, at an external point r * , is principal axes of inertia do not coincide with the axes of symmetry of both boundaries. In the sequence we focus in the particular case in which the inner and the outer boundaries are aligned.
In this case, we can use the same scheme used to calculate the mass of the layer. The principal moments of inertia of the layer, can be written as the subtraction of the principal moments of inertia of two homogeneous ellipsoids of same density di: the principal moments of inertia of one homogeneous ellipsoid of mass m i and the same surface as the outer boundary of the layer, less the principal moments of inertia of the homogeneous ellipsoid of mass m i and the same surface as the inner boundary of the layer. Using the semi axes (A.13) and the masses (A. 19), the principal moments of inertia can be approximated to first order in the flattenings as and its differences are where ∆(fi) = fi − fi−1, denotes the increment of one function fi, between the inner and the outer boundaries of this layer. Using the same scheme used to calculate the mass and the principal moments of inertia, the corresponding gravitational potential of this homogeneous triaxial layer at an external point r * is If we consider the static equilibrium figure, the flattenings can be written as where H k and G k are the Clairaut numbers. Then, the difference of the principal moments of inertia can be approximated to first order in the flattenings as where the parameters Li and L i are The coefficients Li and L i play a role equivalent to the coefficients Hi and Gi for the quantities Ci − Ai, Ci − Bi and Bi − Ai. In this case, the moments of inertia Bi − Ai (resp. Ci − Bi) of the ith layer can be written as the homogeneous moments multiplied by the coefficients Li (resp. L i ), characteristics of this layer. The difference between Li and L i comes from the fact that the body has a differential rotation. If we assume a rigid rotation, then L i = Li(Ω/n) 2 .
The corresponding gravitational potential of this homogeneous triaxial layer at an external point r * is Although we do not calculate the principal moments of inertia when the inner and the outer boundaries are not aligned, it is possible to calculate easily the gravitational potential with the same scheme used to calculate the mass of the layer and the principal moments of inertia. The potential of the layer, can be written as the subtraction of the potential of two homogeneous ellipsoids of same density di: the potential of one homogeneous ellipsoid of mass m i and the same surface as the outer boundary of the layer, given by the Eq. (A.16), less the potential of the homogeneous ellipsoid of mass m i and the same surface as the inner boundary of the layer, given by the Eq. (A.17).
The corresponding gravitational potential is
come from the periodic terms with |j| = 1. The harmonic terms with |j| = 2, do not contribute to the stationary rotation at order e 2 .
It is worth emphasizing that in the absence of friction and gravitational coupling, that is, K = µ = 0, the coefficient κij = δi,j. Then, the non-periodic tidal drift of the ith layer has the same expression that the excess rotation in the case of a homogeneous body, with γi instead of γ (B.28) In the case n/γ1 1, n/γ2 1, an elementary calculation shows that each coefficient κij becomes independent of the friction parameter µ, depending only on the internal structure and on the relaxation factors γ1 and γ2, with fij tending to In the case n/γ1 1, n/γ2 1, each coefficient κij becomes independent of T , K and µ, depending only on the internal structure and on the relaxation factors γ1 and γ2, tending to and the stationary solution tends to synchronous rotation. The periodic terms have amplitudes Bi1 and Bi2, given by the Eq. (B.20). The coefficients cij and sij gives rise to intricate analytical expressions, but are easy to calculate numerically. Fig. 20 shows one example for the Titan's core and the shell constants B1j and B2j, respectively, in function of the shell relaxation factor γ2 (see Table 1 -4). We use that the core relaxation factor is γ1 = 10 −8 s −1 , and fix the ocean's viscosity and thickness values to ηo = 10 −3 Pa s and h = 178 km, respectively. We also plot the non-periodic Ni and periodic Pi terms, separately, and the total tidal drift Bi0 = Ni + Pi. We can observe that if γ2 10 −7.5 s −1 , the shell oscillates around the super-synchronous rotation. When γ2 10 −7.5 s −1 , the tidal drift B20 becomes negative and tends to zero, that is, the shell oscillates around the synchronous rotation, with a period of oscillation equal to the orbital period. The negative sign of the tidal drift B20, is due to the contribution of the periodic terms P2, which becomes negative and |P2| N2. Finally, if γ2 10 −8 s −1 , the amplitude of the shell rotation decreases, tending to zero when γ2 decreases. On the other hand, the core oscillates around the synchronous rotation, with a period of oscillation equal to the orbital period, independently of the shell relaxation factor. This behavior is confirmed by the numerical simulations of non-approximate system (see Sec. 10).
In Fig. 21, we show the comparison of the Titan's shell rotation in the complete non-linear system given by Eq. (B.1) and in the approximate analytical solution given by Eq. (B.19), for some values of the core's relaxation factor γ1 and ocean thickness h. The dashed red lines show the maximum and minimum values of Ω2 − n given by the approximate solution, taking into account only the first harmonic (|j| ≤ 1), while the solid black lines show the maximum and minimum values of Ω2 − n when the complete non-linear system is integrated (using |j| ≤ 7). The approximate solution is in excellent agreement with numerical integration of the equations.

B.2 Atmospheric influence
When we consider the effect of the atmosphere, the rotational system becomeṡ  We assume that the particular solution y1 = C1 cos 2α + S1 sin 2α y2 = C2 cos 2α + S2 sin 2α , (B.33) can be added to (B.6) to obtain the general solutions of the complete equation. Cj and Sj are undetermined coefficients to be obtained by substitution of the parts of the solution into Eq. (B.31) and identification. The derivative of (B.33) isẏ 1 = −2n C1 sin 2α + 2n S1 cos 2α ẏ2 = −2n C2 sin 2α + 2n S2 cos 2α . (B.34) The tidal function can be approximated by  In Fig. 22, we show the same comparison of the Titan's shell rotation in the complete non-linear system and in the approximate analytical solution of the above section. The approximate solution, also is in excellent agreement with numerical integration. It is important to note that the fact that the approximate solution of the non-linear system (B.31) can be expressed as the sum of solutions (B.6) and (B.33), it means that this system has a behavior quasi-linear, at least for the Titan's problem.   where the constants In low-γ approximation (γi n), we can neglet the terms k = 0. If we consider only the terms j = 0, the system becomesẏ (C.5) In the same way in Ferraz-Mello (2015), each solution of this system tends to zero. The role of the terms j = 0 are periodic fluctuations which are the harmonics of the orbital period are that added to the solution. If we consider the terms j = 0, we have that yi 1, and the rotational system iṡ Kij sin (jx), (C. 6) where Kij = (T * ii − T * ii−1 )E2,0E2,j. The solution of this differential equation is Kij j cos (jx), (C. 7) or, in term of the angular velocity, we obtain In high-γ approximation (γi n), we can neglet P ik , then, the system can be written aṡ  which is identical to the system C.5, with k∈Z E 2 2,k instead of E 2 2,0 . Therefore, each solution of this system tends to zero. As in low-γ approximation, the role of the terms j = 0 are periodic fluctuations which are the harmonics of the orbital period are added to the solution. If we consider the terms j = 0, we have that yi 1, and the rotational system isẏ i = − j∈Z j =0 K ij sin (jx), (C. 13) where K ij = (T * ii − T * ii−1 ) k∈Z E 2,k E 2,k+j . Using the solution of the low-γ approximation, then, the angular velocity is γiK ij 2j cos (jnt), (C.14) and the integral (C.1) is zero.