# Andrade rheology in time-domain. Application to Enceladus' dissipation of energy due to forced libration

Yeva Gevorgyan, Gwenaël Boué, Clodoaldo Ragazzo, Lucas Ruiz, Alexandre

C.M. Correia

► To cite this version:

Yeva Gevorgyan, Gwenaël Boué, Clodoaldo Ragazzo, Lucas Ruiz, Alexandre C.M. Correia. Andrade rheology in time-domain. Application to Enceladus’ dissipation of energy due to forced libration. Icarus, 2020, 343, pp.113610. 10.1016/j.icarus.2019.113610. hal-03259889

HAL Id: hal-03259889

https://hal.science/hal-03259889v1

Submitted on 26 Sep 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Andrade rheology in time-domain. Application to Enceladus’ dissipation of energy due to forced libration

Yeva Gevorgyan$^{a,b}$, Gwenaël Boué$^a$, Clodoaldo Ragazzo$^b$, Lucas S. Ruiz$^c$, Alexandre C. M. Correia$^d$

$^a$ASD/IMCCE, CNRS-UMR8028, Observatoire de Paris, PSL University, Sorbonne Université, 77 Avenue Denfert-Rochereau, 75014 Paris, France

$^b$Instituto de Matemática e Estatística, Universidade de São Paulo, 05508-090 São Paulo, SP, Brazil

$^c$Instituto de Matemática e Computação, Universidade Federal de Itajubá, 37500-903 Itajubá, MG, Brazil

$^d$CFisUC, Department of Physics, University of Coimbra, 3004-516 Coimbra, Portugal

## Abstract

The main purpose of this work is to present a time-domain implementation of the Andrade rheology, instead of the traditional expansion in terms of a Fourier series of the tidal potential. This approach can be used in any fully three dimensional numerical simulation of the dynamics of a system of many deformable bodies. In particular, it allows large eccentricities, large mutual inclinations, and it is not limited to quasi-periodic perturbations. It can take into account an extended class of perturbations, such as chaotic motions, transient events, and resonant librations.

The results are presented by means of a concrete application: the analysis of the libration of Enceladus. This is done by means of both analytic formulas in the frequency domain and direct numerical simulations. We do not a priori assume that Enceladus has a triaxial shape, the eventual triaxiality is a consequence of the satellite motion and its rheology. As a result we obtain an analytic formula for the amplitude of libration that incorporates a new correction due to the rheology.

Our results provide an estimation of the amplitude of libration of the core of Enceladus as 0.6% of that of the shell. They also reproduce the observed 10 GW of tidal heat generated by Enceladus with a value of $0.17 \times 10^{14}$ Pa·s for the global effective viscosity under both Maxwell and Andrade rheology.

**Keywords:** Enceladus, Rheology: Andrade, Maxwell, Tides, Librations

## 1. Introduction

Andrade (1910) observed that metallic wires subject to a constant stress tend to stretch as a function of time $t$, according to a law in $t^{1/3}$. He attributed this phenomenon to a rearrangement and/or a rotation of small parts of a crystalline structure switching between successive equilibrium states. It has been shown to accurately model the transient deformation of many metals, but also of polycrystalline ice (Glen, 1955; McCarthy and Castillo-Rogez, 2013, and references therein). A theoretical derivation of the 1/3 power has been proposed by Mott (1953). At high pressure and temperature conditions, rock and rock-forming minerals also behave like metals (Griggs et al., 1960). Laboratory experiments on rock periodically stressed at frequencies spanning $10^{-4}$ to $10^1$ rad/s reproduce an Andrade like creep rheology with a power $\alpha$ varying as a function of the composition between 0.33 and 0.59 (e.g. Goetze, 1971; Goetze and Brace, 1972).

The application of Andrade's rheology in problems of convection in the Earth mantle has been discussed by Stacey (1963). However this law can only be valid over a finite range of frequencies otherwise it would prevent thermal convection and continental drift (Jeffreys, 1972). Geophysical data stemming from body waves, surface waves, free oscillations, and Chandler wobble, covering a large frequency domain also suggest that the mantle dissipates energy

<footnote>Email addresses: yeva@ime.usp.br (Yeva Gevorgyan), gwenael.boue@obspm.fr (Gwenaël Boué), ragazzo@ime.usp.br (Clodoaldo Ragazzo), lucasruiz@unifei.edu.br (Lucas S. Ruiz), acor@uc.pt (Alexandre C. M. Correia)</footnote>
according to $Q^{-1} \propto \omega^{-\alpha}$, where Q is the quality factor, $\omega$ the excitation frequency and $\alpha$ a coefficient between 1/5 and 1/3, in agreement with Andrade's creep model (Anderson and Minster, 1979). This behaviour is observed in the lower mantle for $\omega$ in the range $10^{-8} - 10^{-1}$ rad/s and in the upper mantle for $\omega$ in $10^{-3} - 10^2$ rad/s. For lower and greater frequencies, the expected dissipation rates are $Q^{-1} \propto \omega$ and $Q^{-1} \propto \omega^{-1}$, respectively (ibid.). Anderson and Minster (1979) attribute the plateau feature at intermediary frequencies as the consequence of a continuum spectrum of retardation times, which turns out to be a physical motivation for the construction of the approximate Andrade model presented in Section 5. The width of the frequency band where $Q^{-1}$ evolves as $\omega^{-\alpha}$ is related to the range of dislocation density and grain size involved in grain boundary processes and defect motions (Karato and Spetzler, 1990). More recent geophysical data suggest for the Earth an exponent $\alpha$ of the order of 0.15 (Petit and Luzum, 2010) and an upper frequency of the plateau $Q^{-1} \propto \omega^{-\alpha}$ in the Earth mantle around $10^{-3}$ rad/s (Lau and Faul, 2019).

The Andrade creep model has been introduced in the astronomical community by Efroimsky and Lainey (2007) and Efroimsky (2012b,a). In comparison with more standard rheologies, it can significantly affect the probability of capture in spin-orbit resonance and/or amplify the amount of dissipated energy (Makarov and Efroimsky, 2013; Leconte et al., 2015; Makarov, 2015; Makarov et al., 2016; Walterová and Běhounková, 2017; Renaud and Henning, 2018). The Andrade rheology has been applied in studies of specific celestial bodies such as the Moon (Nimmo et al., 2012; Williams et al., 2014; Williams and Boggs, 2015), Mercury (Padovan et al., 2014; Noyelles et al., 2014; Knibbe and van Westrenen, 2017), Enceladus (Rambaux et al., 2010; Shoji et al., 2013; Běhounková et al., 2013, 2015; Souček et al., 2019), Iapetus (Castillo-Rogez et al., 2011), Io (Bierson and Nimmo, 2016), binary asteroids (Efroimsky, 2015), GJ581 d (Makarov et al., 2012), and Proxima Century b (Ribas et al., 2016).

The main objective of this paper is to present finite dimensional sets of ordinary differential equations that approximate the Andrade rheology with an arbitrary degree of accuracy. Many of the visco-elastic models used in geophysics, Kelvin-Voigt, Maxwell, Burgers, etc, can easily be written in the time domain allowing for applications in cases of large eccentricities or when perturbations are not quasi-periodic (Ferraz-Mello, 2013; Correia et al., 2014; Ferraz-Mello, 2015a,b; Boué et al., 2016; Folonier et al., 2018). In the frequency domain the only difference of the Andrade model from the others is the appearance of a fractional power of the frequency in the complex compliance function. While in the frequency domain this difference is innocuous, in time domain it becomes quite challenging. Indeed, while integral power laws of the frequency are related to ordinary derivatives, the fractional power laws are related to fractional derivatives, which are integro-differential operators. For this reason, as far as we know, equations of motion that describe simultaneously the rotation, the position, and tide deformation under Andrade rheology were never considered in the time domain. Nevertheless, as suggested by Anderson and Minster (1979), Andrade's fractional power law can be interpreted as a continuum of Voigt elements in series, which can itself be approximated by a discrete series of Voigt elements (e.g., Birk and Song, 2010; Ben Jazia et al., 2013). The resulting model, which is known as the extended Burgers model, allows to describe any rheology in the time domain.

To validate the method we analyse the energy dissipation of Enceladus due to the libration which naturally appears in the simulations. In this framework this effect can indeed be independently estimated analytically. This goal allowed us to integrate, for the first time, the equations of motion of a celestial body, Enceladus, using the Andrade rheology to describe tidal deformations.

Enceladus is one of the medium sized moons of Saturn with a radius of 252.1 km, mostly covered by ice. Currently it is synchronised with its host planet. It is orbiting the planet in a slightly eccentric orbit ($e = 0.0045$) due to it being in 2:1 mean motion resonance with another satellite of Saturn, Dione. During the flyby of July 14th, 2005, the Cassini Composite Infrared Spectrometer (CIRS) found plume activity around four large fractures near the south pole of Enceladus (see Porco et al. (2006)). Currently, it is believed that the plum activity on Enceladus is due to tidal heating. No appreciable radiogenic heating is believed to be involved in this process. The latest estimate for the energy flux from the fractures is around 10 GW (Kamata and Nimmo, 2017). According to existing scenarios of satellite formation Canup and Ward (2006); Charnoz et al. (2011); Salmon and Canup (2017); Asphaug and Reufer (2013) radioactive decay no longer has contribution to the generation of the heat on Enceladus (for details see Efroimsky (2018b)).

The Andrade rheology has some parameters that must be fit from observations. In Efroimsky (2018b), using the Fourier expansion to compute the tidal energy dissipation, with additional terms accounting for the forced longitudinal libration, the value found for the mean viscosity was $\eta \approx 0.24 \times 10^{14}$ Pa·s, which is expected for an ice shell. A range $\eta \approx (0.6 - 1.9) \times 10^{14}$ Pa·s was estimated in Folonier et al. (2018). Here we use $\eta = 0.17 \times 10^{14}$ Pa·s both in the direct integrations of the unaverage equations of motion and in the analytical formulae.

The paper is organised as follows.
In Section 2 we first construct non-dissipative equations of motion for the Enceladus-Saturn system. Motion on a fixed orbit is introduced to include the effect of the Enceladus-Dione resonance. In Section 3 we present dissipation in the standard way in terms of the Love number $k_2$. Then in Section 4 the dissipation is introduced into the equations of motion by means of the Association Principle given in Ragazzo and Ruiz (2017). This is the key point in the implementation of the Andrade rheology in time-domain. Both Maxwell and Andrade models are used to describe the rheology of Enceladus. In this way we can evaluate the difference between distinct rheological models on the dissipation of energy. In Section 5 we build up on previous results about finite-dimensional approximations to fractional derivative operators and present a finite dimensional approximation of the Andrade rheology by means of an extended Burgers model. Section 6 is dedicated to forced librations. We present a new analytical expression for the libration magnitude and another for the average energy dissipation rate in order to validate the simulations. In Section 7 we discuss the choice of parameters and initial conditions used in the integrations and present our results. The dissipation rates obtained for different rheologies are compared. In Section 8 we analyse the difference between the angle of libration obtained from the model and that observed. This difference between the two angles and an interior model for Enceladus is used to estimate the libration angle of Enceladus core.

In the Conclusion 9 we summarise the main results in the paper.

## 2. Dynamical model

Consider Enceladus as a deformable body orbiting a point-like Saturn. Let $κ$ be an inertial reference frame at the centre of mass of the system, which is at rest. Let $q'^1 \in κ$ and $q'^2 \in κ$ be the positions, and $m_1$ and $m_2$ the masses of Enceladus and Saturn, respectively. $I_○$ is one third of the trace of the inertia tensor of Enceladus. We suppose that the deformations of the satellite are incompressible which implies that $I_○$ is constant.

In order to describe the rotation of Enceladus, let $K$ be a body frame of Enceladus with origin at its centre of mass. $K$ is defined as being a Tisserand frame, namely with respect to $K$ the angular momentum of Enceladus is null (Munk and MacDonald, 1961). Let $Y: K \to κ$ be the orientation matrix associated with $K$ and $Ω = Y^{-1}\dot{Y} \cdot K \to K$ be the angular velocity operator. Notice that $Y$, and therefore $K$, will be determined by the integration of the equations of motion, they are not known a priori.

The deformations of Enceladus imply that the moment of inertia operator $I: K \to K$ is not constant in time. It can be written as $I = I_○(I - B)$, where $I$ is the identity matrix and $B$ is a symmetric traceless matrix that encodes all the deformation of Enceladus.

Let $q = Y^{-1}(q'^2 - q'^1)$ be the relative position of the planet in the reference frame $K$. The equations for the orbital motion and the spin of the extended body, in the body frame, are (for details see Ragazzo and Ruiz (2017)),

$$ 
\begin{aligned}
\frac{\dot{\mathbf{L}}}{I_○} &= \dot{\Omega} + \dot{\mathbf{B}}\Omega + \Omega\dot{\mathbf{B}} + \mathbf{B}\dot{\Omega} + \Omega\dot{\mathbf{B}} \\
&= [\mathbf{B}, \Omega^2] + \frac{3Gm_2}{|q|^5} [q \otimes q, \mathbf{B}] \\
\dot{q} &= -\Omega q + v \\
\dot{v} &= -\Omega v + G(m_1 + m_2) \left\{ -\frac{1}{|q|^3} q \right. \\
&\qquad \left. - \frac{15}{2} \frac{I_○}{m_1} \frac{1}{|q|^7} (q \cdot \mathbf{B}q)q + 3 \frac{I_○}{m_1} \frac{1}{|q|^5} \mathbf{B}q \right\}
\end{aligned}
\quad (1) $$

where $G = 6.6743 \times 10^{-11} \text{m}^3 \cdot \text{kg}^{-1} \cdot \text{s}^{-2}$ is the gravitational constant, $q$ and $v$ are respectively the position and the velocity vectors, $[\mathbf{B}, \Omega] = \mathbf{B}\Omega - \Omega\mathbf{B}$ denotes the commutator of matrices, and $q \otimes q$ denotes the matrix with entries $(q \otimes q)_{ij} = q_i q_j$. The total angular momentum is given by

$$ \mathbf{L} = I_○(\Omega + \mathbf{B}\Omega + \Omega\mathbf{B}). \quad (2) $$

To keep Enceladus on an eccentric orbit and to attain the forced oscillations regime we would have to add Dione to the model. Here we alternatively fix the eccentricity and the semi-major axis at the present values and remove tidal evolution from the orbits.
It is more convenient to integrate the equations in the inertial reference frame

$$q' = \mathbf{Y}q, \quad v' = \mathbf{Y}v, \quad \omega = \mathbf{Y}\Omega\mathbf{Y}^T, \quad \mathbf{b} = \mathbf{Y}\mathbf{B}\mathbf{Y}^T, \qquad (3)$$

where $^T$ denotes the transpose operator. Notice that for any matrix $\mathbf{x} = \mathbf{Y}\mathbf{X}\mathbf{Y}^T$, its time derivative satisfies

$$\dot{\mathbf{X}} = \mathbf{Y}^T \dot{\mathbf{x}} \mathbf{Y} + \dot{\mathbf{Y}}^T \mathbf{x} \mathbf{Y} + \mathbf{Y}^T \mathbf{x} \dot{\mathbf{Y}} = \mathbf{Y}^T (\dot{\mathbf{x}} + \mathbf{Y} \dot{\mathbf{Y}}^T \mathbf{x} + \mathbf{x} \dot{\mathbf{Y}} \mathbf{Y}^T) \mathbf{Y} = \mathbf{Y}^T (\mathbf{x} - [\omega, \mathbf{x}]) \mathbf{Y}. \qquad (4)$$

In the inertial frame the equations of motion become

$$\begin{aligned} \dot{\omega} + \dot{\mathbf{b}}\omega + \omega\dot{\mathbf{b}} + \mathbf{b}\dot{\omega} + \dot{\omega}\mathbf{b} &= \frac{3Gm_2}{|q'|^5} [q' \otimes q', \mathbf{b}] \\ \dot{q}' &= v' \\ \dot{v}' &= -G(m_1 + m_2) \frac{1}{|q'|^3} q'. \end{aligned} \qquad (5)$$

Equations for $\mathbf{b}$, which are given in the next two sections, are still necessary to close the system.

### 3. Tidal equations in the frequency domain

Let $V'(z,t)$ be the additional quadrupolar potential at time $t$ associated with the deformation of Enceladus at a point $z \in \mathbf{K}$ of its surface, i.e. such that $|z| = R$, where $R$ is the Enceladus' mean radius, and induced by the perturbing quadrupolar potential $W(z,t)$. By definition of the Love number $k_2$, and according to the equivalence principle (Efroimsky, 2018a, Section 6.3),

$$\hat{V}'(z, \omega) = k_2(\omega) \hat{W}(z, \omega), \qquad (6)$$

where $\hat{\cdot}$ denotes the Fourier operator. The matrix $\mathbf{B}(t)$ is proportional to the quadrupolar matrix of Enceladus (Ragazzo and Ruiz, 2017, equation (43)).

$$V'(z,t) = -\frac{3}{2} \frac{\mathrm{GI}_{\circ}}{|z|^5} z^T \mathbf{B}(t) z. \qquad (7)$$

The perturbing quadrupolar potential can be decomposed into a centrifugal term and a tidal term and can be written as

$$W(z,t) = -\frac{1}{2} z^T \mathbf{F}(t) z, \qquad (8)$$

where

$$\mathbf{F}(t) = -\Omega^2 + \frac{1}{3} \mathrm{Tr}(\Omega^2) \mathbb{I} + \frac{3Gm_2}{|q|^5} \left( q \otimes q - \frac{|q|^2}{3} \mathbb{I} \right). \qquad (9)$$

Therefore, equation (6) can be rewritten as

$$3 \frac{\mathrm{GI}_{\circ}}{R^5} \hat{\mathbf{B}}(\omega) = k_2(\omega) \hat{\mathbf{F}}(\omega). \qquad (10)$$

This equation, written in the frequency domain, is not suitable for numerical integrations. In the following section we provide equivalent equations in the time domain for two particular rheologies. The method being very general, it can be applied to arbitrary rheologies.

### 4. The rheologies of Andrade and Maxwell

Linear visco-elastic rheologies are usually represented by spring-dashpot systems. Adding to this system a spring of elastic constant $\gamma$, which represents gravity, and a mass $\mu$, which represents the inertia due to deformations, we can construct a one dimensional oscillator (see Ragazzo and Ruiz (2017) for details). The oscillators associated with the Maxwell and with the Andrade rheologies are shown in Figures 1 and 2, respectively. The total displacement of the
<figure><img src="image_2.png" /><figcaption>Figure 1: Maxwell oscillator. The external force $f(t)$ splits into the force $\lambda(t)$ that acts upon the Maxwell array plus the force $f(t) - \lambda(t)$ that acts upon the $\gamma$ spring.</figcaption></figure>

oscillator is denoted by $x$. The equations for the deformation matrix **B** are introduced by means of the **Association Principle (AP)** derived by Ragazzo and Ruiz (2017)<fnref n="1" />:

> “The differential equation for **B** in the body reference frame **K** is equal to the differential equation for the viscoelastic oscillator after replacing $x$ by **B**”.

### 4.1. Maxwell rheology

The oscillator associated with the Maxwell rheology is given in Figure 1. The equations of motion of the oscillator are

$$ 
\begin{aligned}
\dot{x} &= u \\
\mu \dot{u} &= -\gamma x - \lambda + f \\
\alpha_0^{-1} \dot{\lambda} &= -\eta_0^{-1} \lambda + u,
\end{aligned} 
\quad (11) $$

where $\mu$, $\gamma$, $\alpha_0$ and $\eta_0$ represent, respectively: inertia of deformation, gravity, elasticity and viscosity.

Implementing the AP we obtain the following equations of motion for **B**

$$ 
\begin{aligned}
\dot{\mathbf{B}} &= \mathbf{U} \\
\mu \dot{\mathbf{U}} &= -\gamma \mathbf{B} - \mathbf{\Lambda} + \mathbf{F} \\
J_0 \dot{\mathbf{\Lambda}} &= -\eta_0^{-1} \mathbf{\Lambda} + \mathbf{U},
\end{aligned} 
\quad (12) $$

where $\mathbf{F}(t)$ is the tidal forcing (9), and $J_0 = 1/\alpha_0$. With $\mu = 0$ (see Correia et al. (2018)<fnref n="2" />) the equations of motion (12) can be written

$$ \dot{\mathbf{B}} = \frac{J_0}{1 + \gamma J_0} \dot{\mathbf{F}} + \frac{\eta_0^{-1}}{1 + \gamma J_0} (\mathbf{F} - \gamma \mathbf{B}), \quad (13) $$

Then, applying the change of variables (3) together with the formula (4), we get the expression of the equation of motion in the inertial reference frame, namely,

$$ \dot{\mathbf{b}} = [\omega, \mathbf{b}] + \frac{J_0}{1 + \gamma J_0} (\dot{\mathbf{f}} - [\omega, \mathbf{f}]) + \frac{\eta_0^{-1}}{1 + \gamma J_0} (\mathbf{f} - \gamma \mathbf{b}), \quad (14) $$

where

$$ \mathbf{f} = -\omega^2 + \frac{1}{3} \mathrm{Tr}(\omega^2)\mathbb{I} + \frac{3Gm_2}{|q'|^5} \left( q' \otimes q' - \frac{|q'|^2}{3} \mathbb{I} \right). \quad (15) $$

### 4.2. Andrade rheology

The oscillator associated with the Andrade rheology is given in Figure 2. As will be shown in Section 5 the Andrade element can be approximated by an extended Burgers oscillator with $n$ Voigt elements given in Figure 3,

<footnote n="1">¹The equivalence of the AP and the Correspondence Principle (Efroimsky, 2012b) is addressed in Section 4 of Correia et al. (2018). The main difference of the Association Principle and the Correspondence Principle is that the first is formulated in the time domain, while the second is formulated in the frequency domain.</footnote>

<footnote n="2">²In Correia et al. (2018) it was shown that deformation inertia is negligible if $\omega \ll \omega_0$, where $\omega_0 = \sqrt{\gamma/\mu}$ is the natural frequency of oscillation of the system when damping is neglected and $\omega$ is the angular frequency of a harmonic tidal force.</footnote>
<figure><img src="image_3.png" /><figcaption>Figure 2: The Andrade oscillator formulated from the Andrade model, $\alpha$ represents the so called spring-pot element. The external force $f(t)$ splits into the force $\lambda(t)$ that acts upon the Andrade array plus the force $f(t) - \lambda(t)$ that acts upon the $\gamma$ spring.</figcaption></figure>

where $x = x_0 + \tilde{x}_0 + x_1 + \dots + x_n$ is the total displacement of the oscillator. The equations of motion of the oscillator are

$$
\begin{align}
\dot{x} &= u \nonumber \\
\mu \dot{u} &= -\gamma x - \lambda + f \nonumber \\
J_0 \dot{\lambda} &= -\left(\sum_{r=0}^{n} \eta_r^{-1}\right) \lambda + u + \sum_{r=1}^{n} \tau_r^{-1} x_r \tag{16} \\
\dot{\tau}_1 &= -\tau_1^{-1} x_1 + \eta_1^{-1} \lambda \nonumber \\
\dot{\tau}_2 &= -\tau_2^{-1} x_2 + \eta_2^{-1} \lambda \nonumber \\
\dots \nonumber \\
\dot{x}_n &= -\tau_n^{-1} x_n + \eta_n^{-1} \lambda, \nonumber
\end{align}
$$

where $\eta_r$ is the viscosity of the r-th dashpot, $\tau_r = \eta_r/\alpha_r$ is the relaxation time of the r-th Voigt element, $r = 1, \dots, n$, and we used the relation

$$ \dot{x}_0 + \dot{\tilde{x}}_0 = \dot{x} - (\dot{x}_1 + \dots + \dot{x}_n) = \frac{\dot{\lambda}}{\alpha_0} + \frac{\lambda}{\eta_0}. $$

By means of the AP we find the equations for $\mathbf{B}$, $\mathbf{\Lambda}$, and $\mathbf{B}_r$, $r = 1, \dots, n$ ($\mathbf{\Lambda}$ and $\mathbf{B}_r$ are symmetric traceless matrices):

$$
\begin{align}
\dot{\mathbf{B}} &= \mathbf{U} \nonumber \\
\mu \dot{\mathbf{U}} &= -\gamma \mathbf{B} - \mathbf{\Lambda} + \mathbf{F} \nonumber \\
J_0 \dot{\mathbf{\Lambda}} &= -(\sum_{r=0}^n \eta_r^{-1}) \mathbf{\Lambda} + \mathbf{U} + \sum_{r=1}^n \tau_r^{-1} \mathbf{B}_r \nonumber \\
\dot{\mathbf{B}}_1 &= -\tau_1^{-1} \mathbf{B}_1 + \eta_1^{-1} \mathbf{\Lambda} \nonumber \\
\dot{\mathbf{B}}_2 &= -\tau_2^{-1} \mathbf{B}_2 + \eta_2^{-1} \mathbf{\Lambda} \nonumber \\
\dots \nonumber \\
\dot{\mathbf{B}}_n &= -\tau_n^{-1} \mathbf{B}_n + \eta_n^{-1} \mathbf{\Lambda}, \tag{17}
\end{align}
$$
<figure><img src="image_4.png" /><figcaption>Figure 3: Extended Burgers oscillator associated with the Andrade model.</figcaption></figure>

where **F** is the tidal forcing in equation (9). In the inertial reference frame and with $\mu = 0$ these equations become

$$ 
\begin{aligned}
\dot{\mathbf{b}} &= [\omega, \mathbf{b}] + \frac{J_0}{1 + \gamma J_0} (\dot{\mathbf{f}} - [\omega, \mathbf{f}]) \\
&\quad + \frac{\sum_{r=0}^{n} \eta_r^{-1} (\mathbf{f} - \gamma \mathbf{b})}{1 + \gamma J_0} - \frac{\sum_{r=1}^{n} \tau_r^{-1} \mathbf{b}_r}{1 + \gamma J_0} \\
\dot{\mathbf{b}}_1 &= [\omega, \mathbf{b}_1] - \tau_1^{-1} \mathbf{b}_1 + \eta_1^{-1} (\mathbf{f} - \gamma \mathbf{b}) \\
\dot{\mathbf{b}}_2 &= [\omega, \mathbf{b}_2] - \tau_2^{-1} \mathbf{b}_2 + \eta_2^{-1} (\mathbf{f} - \gamma \mathbf{b}) \\
\dots \\
\dot{\mathbf{b}}_n &= [\omega, \mathbf{b}_n] - \tau_n^{-1} \mathbf{b}_n + \eta_n^{-1} (\mathbf{f} - \gamma \mathbf{b}),
\end{aligned}
\quad (18) $$

with **f** from equation (15).

## 5. An approximation to the Andrade rheology

In this Section we show that the spring-pot element in the spring-dashpot representation of the Andrade rheology, Fig. 2, can be replaced by a continuum of Voigt elements in series and further approximated by a finite set of *n* Voigt elements as in Fig. 3. In the particular example of this paper we show that *n* can be taken equal to six. The same procedure can be done in any other situation to determine a minimum number of Voigt elements *n* that can be used in the approximation.

### 5.1. The Andrade model

The Andrade (1910) model is represented by three mechanical elements combined in series: a spring, a dashpot and a spring-pot (fig. 2). The creep function of the model is

$$ J(t) = \left[ J_0 + \frac{t}{\eta_0} + At^{\alpha} \right] H(t), \quad (19) $$

where *H* is the Heaviside function, and $A > 0$ and $0 < \alpha < 1$ are physical parameters. In order to avoid a fractional dimension of time, Efroimsky (2012b) replaces the parameter *A* by $\tau_A$, the so called “Andrade time” (with dimension of time):

$$ A = J_0 \tau_A^{-\alpha} = \alpha_0^{-1} \tau_A^{-\alpha}. \quad (20) $$

In terms of the Andrade and Maxwell $\tau_M = \eta_0/\alpha_0$ times the creep function becomes

$$ J(t) = J_0 \left[ 1 + \frac{t}{\tau_M} + \left( \frac{t}{\tau_A} \right)^{\alpha} \right] H(t). \quad (21) $$

This creep function induces the following complex compliance in the frequency domain

$$ \hat{J}(\omega) = J_0 \left[ 1 + (i\omega\tau_A)^{-\alpha} \Gamma(1+\alpha) - i(\omega\tau_M)^{-1} \right]. \quad (22) $$
## 5.2. The representation of the Andrade model by an extended Burgers model

The creep function of a finite array of n Voigt elements is (Bland, 2016)

$$J(t) = \sum_{r=1}^{n} J_r [1 - \exp(-t/\tau_r)] H(t), \qquad (23)$$

where $J_r$ and $\tau_r = \eta_r J_r$ are the compliance and the relaxation time of the r$^{th}$ Voigt element, respectively. The complex compliance of this finite array is

$$i\omega\hat{J}(\omega) = \sum_{r=1}^{n} \frac{J_r}{1 + i\omega\tau_r}. \qquad (24)$$

The continuum limit of this array is well defined as $n \to \infty$ provided that $\sum_{r=1}^{n} J_r$ converges. In the limit $n \to \infty$, the creep function becomes

$$J(t) = H(t) \int_{0}^{\infty} j(\tau) [1 - \exp(-t/\tau)] d\tau. \qquad (25)$$

In order to obtain the continuous distribution of compliances associated with the Andrade model it is necessary to solve the equation

$$At^{\alpha} = \int_{0}^{\infty} j(\tau) [1 - \exp(-t/\tau)] d\tau \qquad (26)$$

for $j(\tau)$. The solution to the equation is

$$j(\tau) = \frac{A\alpha}{\Gamma(1-\alpha)} \tau^{-1+\alpha}, \quad 0 < \alpha < 1. \qquad (27)$$

Using the last two equations the Fourier transform of $At^\alpha H(t)$ can be written as

$$\int_{-\infty}^{\infty} At^{\alpha} H(t) e^{-i\omega t} dt = \frac{A\alpha}{\Gamma(1-\alpha)} \int_{0}^{\infty} \frac{\tau^{-1+\alpha}}{i\omega - \tau\omega^2} d\tau. \qquad (28)$$

Then the creep function becomes

$$J(t) = \left[ J_0 + \frac{t}{\eta_0} + \frac{A\alpha}{\Gamma(1-\alpha)} \int_0^\infty \tau^{-1+\alpha} [1 - \exp(-\frac{t}{\tau})] d\tau \right] H(t) \qquad (29)$$

and the complex compliance becomes

$$\hat{J}(\omega) = J_0 + (i\eta_0\omega)^{-1} + \frac{A\alpha}{\Gamma(1-\alpha)} \int_0^\infty \frac{\tau^{-1+\alpha}}{1+i\tau\omega} d\tau. \qquad (30)$$

The main idea in the discretisation of the Andrade model is to approximate the integral in the last term of equation (30) by a finite sum, equation (24), that represents the complex compliance of an extended Burgers model, namely

$$\frac{A\alpha}{\Gamma(1-\alpha)} \int_{0}^{\infty} \frac{\tau^{-1+\alpha}}{1+i\tau\omega} d\tau \approx \sum_{r=1}^{n} \frac{J_r}{1+i\omega\tau_r}. \qquad (31)$$

One way to do this is described in Birk and Song (2010). The compliance $J_r$ and the relaxation time $\tau_r$ of the Voigt elements in equation (31) for $r = 1, \dots, n$ are given by

$$J_r = \frac{8A\alpha}{\Gamma(1-\alpha)(1-\bar{q}_r)^4} \frac{\bar{\lambda}_r}{\tau_r} = \frac{(1+\bar{q}_r)^4}{(1-\bar{q}_r)^4}, \qquad (32)$$

where $\bar{\lambda}_r$ and $\bar{q}_r$ denote the weights and abscissas, respectively, of the n-point Gauss-Jacobi quadrature rule, with weighting function $(1 - \bar{q})^{1-2\alpha}(1 + \bar{q})^{-1+4\alpha}$.
<figure><img src="image_5.png" /><figcaption>Figure 4: Approximation of the Andrade rheology by an extended Burgers model with 20 Voigt elements. The phase lag (top) and the compliance imaginary part (bottom) are represented as a function of the excitation frequency $\omega$ in units of the orbital frequency $n$. The black vertical line marks the value of Enceladus mean rotation rate.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 5: Same as Fig. 4 but for an extended Burgers model made of 6 Voigt elements only.</figcaption></figure>
### 5.2.1. The number of Voigt elements

In the following, we fix the parameters of the Andrade and Maxwell models as $\alpha = 0.33$, $J_0 = 1/\alpha_0 = 10^{10}$ Pa, $\eta_0 = 0.17 \times 10^{14}$ Pa·s, and $\tau_M = \tau_A = 1700$s (see next section). In this section these values and the dominant tidal forcing frequency upon Enceladus are used as prototypes. The goal is to show how to estimate a number of Voigt elements in an extended Burgers model that guarantees a good approximation for the Andrade rheology in the frequency range of the tidal forcing.

The complex compliance of an extended Burgers model that approximates the Andrade model is

$$ \hat{J}(\omega) = J_0 + (i\eta_0\omega)^{-1} + \sum_{i=1}^{n} J_i \left( \frac{1}{1 + i\omega\tau_i} \right), \qquad (33) $$

where $J_i$ and $\tau_i$ are given in equation (32). The phase lag associated with this complex compliance is

$$ \tan \delta(\omega) = - \frac{\Im[\hat{f}(\omega)]}{\Re[\hat{f}(\omega)]} = \frac{\frac{1}{\tau_m \omega} + \sum_{i=1}^{n} \frac{J_i}{J_0 (1 + (\omega \tau_i)^2)} \omega \tau_i}{1 + \sum_{i=1}^{n} \frac{J_i}{J_0 (1 + (\omega \tau_i)^2)} \omega^2}. \qquad (34) $$

In Figure 4, we show the phase lag and the imaginary part of the complex compliance deduced from the Andrade model on the one hand and from the extended Burgers model with 20 Voigt elements on the other as a function of the frequency. In Figure 5 we show the same comparisons replacing the extended Burgers model with 20 Voigt elements by another with 6 elements. As can be seen from these figures, the use of an extended Burgers rheology with 6 Voigt elements is sufficient to approximate the Andrade rheology within the range $0.1 \le \omega/n \le 2000$. In this case, this is enough because the main fundamental frequency of tidal forcing of Enceladus by Saturn actually is $n = 5.308 \times 10^{-5}$ rad·s$^{-1}$ (the orbital frequency of Enceladus) therefore all harmonics of order $k$ with $1 \le k \le 2000$ are within the region where the approximation is good.

## 6. Expected forced librations

Forced librations may be an important source of energy dissipation in Enceladus (see Efroimsky (2018a,b)). In this section we derive analytical expressions for the amplitude and dissipation power of forced librations under the hypotheses that the amplitude is small and the spin axis is perpendicular to the orbital plane. This calculation will allow us to validate the numerical results presented in Section 7.

### 6.1. Amplitude of libration

The satellite frame $\mathbf{K}$ and the inertial frame $\mathbf{\kappa}$, which are defined by the orthonormal vectors ($E_1, E_2, E_3$) and ($e_1, e_2, e_3$), respectively, are chosen to have $E_3$ and $e_3$ aligned with the satellite spin axis. The angle from $e_1$ to $E_1$ is denoted by $\theta$. If $M$ is the mean anomaly of the satellite, then $\theta = M + \zeta$, where $\zeta$ is the libration angle. The angular velocity is $\dot{\theta} = \Omega_3 = n + \dot{\zeta}$, where $n$ is the mean rotation rate that is equal to the mean motion$^3$.

We now focus on the first equation in (1), where $\mathbf{B}$ is determined by the generic equation (10) that with the notation $K_2(\omega) = \frac{R^5}{3GI_0} k_2(\omega)$ becomes

$$ \hat{\mathbf{B}}(\omega) = K_2(\omega) \hat{\mathbf{F}}(\omega). \qquad (35) $$

As it was shown by Correia et al. (2018) (Section 4, in particular equation (31)), the Association Principle in the frequency domain implies

$$ K_2(\omega) = \frac{1}{\gamma + \hat{J}^{-1}(\omega) - \mu\omega^2}, \qquad (36) $$

<footnote>$^3$At low viscosity the rotation becomes asynchronous with an average angular velocity given by $\dot{\theta} \approx n(1 + 6e^2)$ (e.g., Folonier et al., 2018). Here this excess of rotation speed is neglected because calculations are limit to first order in eccentricity.</footnote>
an expression that holds for arbitrary rheology. If the inertia of deformation is neglected ($\mu = 0$), as assumed in this paper, then<fnref n="4" />

$$K_2(\omega) = \frac{1}{\gamma + \hat{f}^{-1}(\omega)}. \qquad (37)$$

For instance, for a Maxwell rheology the complex compliance is $\hat{f}(\omega) = \alpha_0^{-1} + (i\omega\eta_0)^{-1}$ and equation (37) implies

$$K_2(\omega) = \frac{1}{\gamma} \frac{1 + i\omega\tau_M}{1 + i\omega(\tau_M + \tau_F)}, \qquad (38)$$

where

$$\tau_M = \eta_0 / \alpha_0 \quad \text{and} \quad \tau_F = \eta_0 / \gamma$$

are the Maxwell time and the fluid time (only due to gravity and viscosity), respectively.

The forcing term $\mathbf{F}$ is given in equation (9). It contains a term that represents the centrifugal stress, $-\Omega^2 + \frac{1}{3}\text{Tr}(\Omega^2)\mathbb{I}$, and another that represents the tidal stress. The centrifugal part is given by a diagonal matrix with equal elements in the first and second lines. This term has no effect on the libration and will be neglected. So,

$$\mathbf{F}(t) = \frac{3Gm_2}{|q|^5} \left( q \otimes q - \frac{|q|^2}{3} \mathbb{I} \right), \qquad (39)$$

which is the tidal stress. The position of Saturn in $\mathbf{K}$ is given by

$$q = r \begin{pmatrix} \cos(f - \theta) \\ \sin(f - \theta) \\ 0 \end{pmatrix}, \qquad (40)$$

where $f$ is the true anomaly of Enceladus and $r$ is the distance from Enceladus to Saturn. The first equation in (1) can be simplified. Since the motion is planar, $B_{13} = B_{23} = 0$ and $[\mathbf{B}, \Omega^2] = 0$. Moreover, $[\mathbb{I}, \mathbf{B}] = 0$ implies that the last term in right hand side of equation the first equation in (1) can be written as $[\mathbf{F}, \mathbf{B}]$. So, the first equation in (1) becomes

$$\begin{gathered} \dot{\Omega} + \dot{\mathbf{B}}\Omega + \Omega\dot{\mathbf{B}} + \mathbf{B}\dot{\Omega} + \dot{\Omega}\mathbf{B} = [\mathbf{F}, \mathbf{B}] \quad \text{or} \\ \dot{\Omega}_3 = \frac{(B_{11} - B_{22})F_{12} + B_{12}(F_{22} - F_{11}) + \Omega_3\dot{\mathbf{B}}_{33}}{1 - B_{33}}. \end{gathered} \qquad (41)$$

Equations (35), (39), (40), and (41) determine the librations of Enceladus.

If the eccentricity of the orbit $e$ is zero then there is a trivial solution $\zeta(t) = 0$. If $e > 0$ is small then there exists a nontrivial solution $\zeta$ that has an amplitude proportional to $e$ and a period equal to the orbital period. In order to determine this solution up to order $e$ we use the following well known relations (Yoder, 1995) (terms of order $e^2$ will be neglected in the following equations):

$$\frac{a}{r} = 1 + e \cos M, \quad f - M = 2e \sin M, \qquad (42)$$

where $a$ is the semi-major axis. Using these relations, equations (39) and (40), and $f - \theta = f - M - (\theta - M) = 2e \sin M - \zeta$ we can write the non trivial components of $\mathbf{F}$ as

$$\begin{align} F_{11} &= c_1(2/3 + 2e \cos M) \nonumber \\ F_{12} &= c_1(2e \sin M - \zeta) = F_{21} \nonumber \\ F_{22} &= -c_1(1/3 + e \cos M) = F_{33} \nonumber \\ c_1 &= \frac{3Gm_2}{a^3}. \tag{43} \end{align}$$

$$k_2(n) = \frac{3}{2} \frac{1}{1 + \frac{57}{8\pi} \frac{1}{G\rho^2 R^2} \frac{1}{J}}$$

used by Efroimsky (2015) (Eq.36). Detailed comparison of two notations can be found in (Correia et al., 2018, section 4).

<footnote n="4">⁴This expression is equivalent to $K_2(n) = \frac{R^5}{3Gd_0} k_2(n)$, where</footnote>
Notice that **F** is the sum of a constant term plus an oscillating term with period $2\pi/n$. So, **B** must be of the same form. Solving equation (35) we obtain that the constant elements of **B** are

$$ \bar{B}_{11} = K_2(0) c_1 \frac{2}{3}, \quad \bar{B}_{22} = -K_2(0) c_1 \frac{1}{3} = \bar{B}_{33}. \qquad (44) $$

It is convenient to write the oscillating part of **B** and other quantities in complex notation

$$ \begin{aligned} \bar{B}_{jk} &= \hat{B}_{jk} \exp(iM) + \text{c.c.}, \\ \tilde{F}_{jk} &= \hat{F}_{jk} \exp(iM) + \text{c.c.}, \\ \tilde{\zeta} &= \hat{\zeta} \exp(iM) + \text{c.c.}, \end{aligned} \qquad (45) $$

where c.c. stands for complex conjugate. With these definitions, equation (35) and $\dot{M} = n$ imply

$$ \hat{B}_{jk} = K_2(n) \hat{F}_{jk}. \qquad (46) $$

Equation (43) shows that $\hat{F}_{12}$ is of the order of $e$ and so is $\hat{B}_{12}$, due to equation (46). Therefore up to the order of $e$ equation (41) becomes

$$ \begin{aligned} -n^2 \hat{\zeta}(1 - \bar{B}_{33}) &= in^2 \hat{B}_{33} + (\bar{B}_{11} - \bar{B}_{22}) \hat{F}_{12} + \hat{B}_{12}(\bar{F}_{22} - \bar{F}_{11}) \implies \\ -n^2 \hat{\zeta}(1 + \frac{1}{3}c_1 K_2(0)) &= in^2 K_2(n) \hat{F}_{33} + c_1 [K_2(0) \hat{F}_{12} - \hat{B}_{12}] \implies \\ -n^2 \hat{\zeta}(1 + \frac{1}{3}c_1 K_2(0)) &= in^2 K_2(n) \hat{F}_{33} + c_1 [K_2(0) - K_2(n)] \hat{F}_{12}. \end{aligned} \qquad (47) $$

Equation (43) implies that $\hat{F}_{12} = -c_1(ie + \hat{\zeta})$ and $\hat{F}_{33} = -c_1e/2$. This and equation (47) imply that the amplitude of libration $\beta$ is given by

$$ \beta = 2|\hat{\zeta}|, \quad \hat{\zeta} = \frac{\frac{1}{2}ic_1K_2(n) - i(c_1/n)^2[K_2(n) - K_2(0)]}{1 + \frac{1}{3}c_1K_2(0) + (c_1/n)^2[K_2(n) - K_2(0)]} e. \qquad (48) $$

Notice that the factor two accounts for the two exponential terms in the expression of $\tilde{\zeta}$ (45) with frequencies $n$ and $-n$, respectively.

## 6.2. Amplitude of forced libration: comparison with other results in the literature

In this paper we do not a priori assume that Enceladus has a triaxial shape (the same is done in Folonier et al. (2018)). The eventual triaxiality is a consequence of the satellite motion and its rheology. So, before comparing our results with others in the literature, it is necessary to analyse the mean moments of inertia of Enceladus given in equation (44). The following results hold for any rheology such that $J^{-1}(0) = 0$ and whenever the inertia of deformation is neglected, $\mu = 0$. The first hypothesis is valid for the Maxwell and Andrade rheologies and for any other for which the static shape of the body is determined solely by gravity.

If $\bar{A} < \bar{B} < \bar{C}$ denote the three principal moments of inertia of Enceladus then the definition of **B** implies

$$ \frac{\bar{B} - \bar{A}}{\bar{C}} = \frac{\bar{B}_{11} - \bar{B}_{22}}{1 - \bar{B}_{33}} \approx \bar{B}_{11} - \bar{B}_{22}, \qquad (49) $$

where we used that $B_{kj} \ll 1$ for all $k, j$. The gravitational modulus $\gamma$ is related to the zero secular Love number by (Ragazzo and Ruiz, 2017, equation (14))

$$ \frac{1}{\gamma} = \frac{R^5}{3GI_0} k_2(0) = K_2(0). \qquad (50) $$

The combination of equations (49), (50), and (44) gives

$$ \frac{\bar{B} - \bar{A}}{\bar{C}} \approx \frac{c_1}{\gamma}. \qquad (51) $$
An algebraic manipulation shows that the term in the numerator of equation (48) can be written as

$$ \frac{1}{3}c_1 K_2(0) + (c_1/n)^2 [K_2(n) - K_2(0)] = \frac{1}{3}\frac{c_1}{\gamma} + 3\frac{c_1}{\gamma}\frac{m_2}{m_1+m_2}\left[\frac{k_2(n)}{k_2(0)} - 1\right]. \quad (52) $$

In Correia et al. (2018) (Section 5 Proposition 1), under the hypothesis $\mu = 0$, it is shown that $|k_2(n)/k_2(0)| \le 1$ for any $n \ge 0$. This implies that the quantity in equation (52) is of the order of $c_1/\gamma$ that is much smaller than one due to equation (51). Therefore equation (48) implies

$$ \hat{\xi} \approx \left( \frac{1}{2} i c_1 K_2(n) - i (c_1/n)^2 [K_2(n) - K_2(0)] \right) e \quad (53) $$

and, using equations (51), (52) and $K_2(n) = \frac{k_2(n)}{k_2(0)} \frac{1}{\gamma}$

$$ \hat{\xi} \approx i \frac{1}{2} e \frac{\bar{B} - \bar{A}}{\bar{C}} \frac{k_2(n)}{k_2(0)} + i 3 e \frac{\bar{B} - \bar{A}}{\bar{C}} \frac{m_2}{m_1 + m_2} \left[ 1 - \frac{k_2(n)}{k_2(0)} \right]. \quad (54) $$

In the limit $m_1 \ll m_2$, the amplitude of libration $\beta = 2|\hat{\xi}|$ of a deformable body in hydrostatic equilibrium is thus

$$ \beta \approx 6e \frac{\bar{B} - \bar{A}}{\bar{C}} \left| 1 - \frac{5 k_2(n)}{6 k_2(0)} \right|, \quad (55) $$

which coincides with the result obtained by Van Hoolst et al. (2013). This expression is compared with other results in literature by Noyelles (2017).

The classical formula for the amplitude of forced libration of a rigid triaxial body under the same circumstances is (Danby (1962) also Efraimsky (2018a) equation (13) or Thomas et al. (2016) equation (1))

$$ \beta = 2|\hat{\xi}| \approx 6e \frac{\bar{B} - \bar{A}}{\bar{C}}. \quad (56) $$

So our formula (55) shows that the amplitude of libration of a deformable body under hydrostatic equilibrium must be smaller than that of a rigid body with the same principal moments of inertia $\bar{A}, \bar{B}, \bar{C}$, by a factor $|1 - 5k_2(n)/6k_2(0)|$. The amplitudes are the same only if the complex number $k_2(n)$ is zero.

The case $k_2(n) = k_2(0)$, which holds for a body made of a perfect fluid under self-gravity (no dissipation of energy) and no inertia of deformation ($\mu = 0$), has an interesting interpretation that reveals the essence of a “Tisserand frame”. In this case the satellite do have a tidal deformation, which is easily seen from equation (46), with a principal axis aligned with the position axis of the planet. No torque acts upon the satellite and the satellite angular momentum remains constant. Nevertheless, the amplitude of the tidal deformation depends on the distance from the satellite to the planet. As this distance varies the moment of inertia of the satellite also varies and conservation of angular momentum requires a variation of angular velocity. Equation (55) shows that the corresponding amplitude of libration of the Tisserand frame $K$ is 1/6 of the amplitude of libration of a rigid body.

### 6.3. Energy dissipation rate

The energy dissipation rate averaged over an orbital period $T = 2\pi/n$ is given by (see Ragazzo and Ruiz (2017) appendix 2.1):

$$ \frac{\Delta E}{T} = \frac{I_{\circ}}{2} \frac{1}{T} \int_{0}^{T} \mathrm{Tr}(\mathbf{F}\dot{\mathbf{B}})dt. $$

Then equations (45) and (46) imply

$$ \frac{\Delta E}{T} = -n I_{\circ} \mathfrak{S} [K_2(n)] \sum_{jk} |\hat{F}_{jk}|^2. \quad (57) $$

Notice that the imaginary part of the Love number is negative, therefore the energy dissipated $\Delta E$ is positive. From equation (43) we obtain $\hat{F}_{11} = c_1 e$, $\hat{F}_{22} = \hat{F}_{33} = -c_1 e/2$, $\hat{F}_{13} = \hat{F}_{31} = \hat{F}_{23} = \hat{F}_{32} = 0$, and $\hat{F}_{12} = \hat{F}_{21} = \hat{F}_{12} = -c_1 (ie + \hat{\xi})$.
Substituting these expressions into equation (57) and replacing $c_1$ and $K_2(n)$ by their expressions, we arrive at the following formula:

$$ \frac{\Delta E}{T} = -\Im[k_2(n)] \frac{(nR)^5}{G} \left( \frac{m_2}{m_1+m_2} \right)^2 \times \left[ \frac{21}{2}e^2 + 12e\Im(\hat{\zeta}) + 6|\hat{\zeta}|^2 \right]. \quad (58) $$

### 6.4. Energy dissipation: comparison with other results in the literature

Again the comparison of our results with others in the literature is not straightforward because we do not assume the body has an a priori triaxial shape. So, in order to do the comparison we first suppose that the body librates as if it were rigid with principal moments of inertia $\bar{A} < \bar{B} < \bar{C}$. In this case the complex amplitude of libration $\hat{\zeta}$ is given by equation (54) with $k_2(n)/k_2(0) = 0$, namely $\hat{\zeta} = i\beta/2$, where $\beta$ is the classical amplitude of libration in equation (56). Substituting this value of $\hat{\zeta}$ into equation (58) we obtain

$$ \frac{\Delta E}{T} = -\Im[k_2(n)] \frac{(nR)^5}{G} \left( \frac{m_2}{m_1+m_2} \right)^2 \times \left[ \frac{21}{2}e^2 + 6e\beta + \frac{3}{2}\beta^2 \right], \quad (59) $$

which corresponds to the expression in Efroimsky (2018a) (equation (96)). This expression with $\beta = 0$ is that in Segatz et al. (1988) (equation (13)).

## 7. Integration of the equations of motion

### 7.1. Numerical settings

In this section we specify all the constants and initial conditions necessary to integrate the equations of motion. The numerical integrations of the system of first order ODEs are carried out by a Runge-Kutta method of order 8 with embedded error estimator of order 7 due to Dormand & Prince, with stepsize control (Hairer et al., 1993).

#### 7.1.1. The choice of parameters

The present-time physical and orbital parameters for Enceladus and Saturn are given in Table 1. The constant $\gamma$ is calculated using (Ragazzo, 2018)

$$ \gamma = \frac{4 G m_1}{5 R_1^3}, \quad I_\circ = \frac{2}{5} m_1 R_1^2, $$

where $m_1$ is the mass and $I_\circ$ is the mean moment of inertia of the satellite.

The physical dimensions of $\eta_0$ and $\alpha_0 = 1/J_0$ in equations (12) are 1/s and 1/s$^2$, respectively. As discussed in Correia et al. (2018), if these constants are divided by

$$ \frac{152\pi R}{15 m_1} = 7.43107 \times 10^{-14} \, \text{m} \cdot \text{kg}^{-1}, \quad (60) $$

then the resulting parameters have the usual dimensions ($\eta_0$[Pa·s] and $\alpha_0$[Pa]) and represent the microscopic viscosity and rigidity of an equivalent homogeneous body (possibly with radius different from $R$) with the same macroscopic properties (in this work at no time Enceladus is assumed to be homogeneous). Therefore these re-scaled constants can be directly compared to those used by Efroimsky (2018b), namely $\alpha_0 = 10^{10}$ Pa and $\eta_0 = 0.17 \times 10^{14}$ Pa·s. This gives the Maxwell time $\tau_M = \eta_0/\alpha_0 = 1700$ s. For the Andrade rheology we choose $\alpha = 0.33$ and $\tau_A = \tau_M = 1700$ s (Shoji et al., 2013). As discussed in Section 5, for this choice of parameters the number of Voigt elements in the extended Burgers model can be taken as $n = 6$. The rheological parameters associated with this extended Burgers model are given in Table 2. They were determined using equation (32) (see Birk and Song (2010) for details).
<table><caption>Table 1: The present-time physical and orbital parameters of Enceladus</caption><thead><tr><th>Parameter</th><th>Notation</th><th>Units</th><th>Values</th></tr></thead><tbody><tr><td>Mass of the host (Saturn)</td><td>m<sub>2</sub></td><td>kg</td><td>5.683 &times; 10<sup>26</sup></td></tr><tr><td>Mass</td><td>m<sub>1</sub></td><td>kg</td><td>1.08 &times; 10<sup>20</sup></td></tr><tr><td>Mean radius</td><td>R</td><td>m</td><td>252.1 &times; 10<sup>3</sup></td></tr><tr><td>Eccentricity</td><td>e</td><td></td><td>0.0045</td></tr><tr><td>Semi-major axis</td><td>a</td><td>m</td><td>2.38 &times; 10<sup>8</sup></td></tr><tr><td>Mean rotation rate</td><td>n</td><td>rad &middot; s<sup>-1</sup></td><td>5.308 &times; 10<sup>-5</sup></td></tr><tr><td>Orbital period</td><td></td><td>s</td><td>1.184 &times; 10<sup>5</sup></td></tr><tr><td>Mean moment of inertia</td><td>I<sub>o</sub></td><td>kg &middot; m<sup>2</sup></td><td>2.2994 &times; 10<sup>30</sup></td></tr><tr><td>Rigidity constant</td><td>&gamma;</td><td>s<sup>-2</sup></td><td>4.69504 &times; 10<sup>-7</sup></td></tr></tbody></table>

The mass of Saturn is taken from NASA Planetary Fact Sheets. The moment of inertia and the rigidity constant are calculated. All the other parameters are taken from Table 1 in Efroimsky (2018b).

<table><caption>Table 2: Rheological parameters for the extended Burgers model.</caption><thead><tr><th>r</th><th>1/τ<sub>r</sub>(yr<sup>-1</sup>)</th><th>1/η<sub>r</sub>(yr)</th></tr></thead><tbody><tr><td>0</td><td>18550.6</td><td>2.510112 × 10<sup>-8</sup></td></tr><tr><td>1</td><td>2378870</td><td>1.95734 × 10<sup>-7</sup></td></tr><tr><td>2</td><td>669017</td><td>8.89215 × 10<sup>-8</sup></td></tr><tr><td>3</td><td>170566</td><td>3.88775 × 10<sup>-8</sup></td></tr><tr><td>4</td><td>37423.7</td><td>1.5879 × 10<sup>-8</sup></td></tr><tr><td>5</td><td>6514.34</td><td>5.79894 × 10<sup>-9</sup></td></tr><tr><td>6</td><td>782.589</td><td>1.7528 × 10<sup>-9</sup></td></tr></tbody></table>

We use quantities in Table 3 to change from MKS units to Year.

### 7.1.2. Initial conditions for a Maxwell rheology

The libration of Enceladus, under Maxwell rheology and moving on a fixed slightly eccentric orbit, is determined by equations (5) and (14). In order to integrate these equations we use the following set of initial conditions:

$$ 
\begin{aligned}
q'(0) &= \begin{pmatrix} 0.001598 \\ 0 \\ 0 \end{pmatrix} \text{au}, v'(0) = \begin{pmatrix} 0 \\ 2.649 \\ 0 \end{pmatrix} \text{au} \cdot \text{yr}^{-1}, \\
\omega(0) &= \begin{pmatrix} 0 & -1672.44 & 0 \\ 1672.44 & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix} \text{yr}^{-1}, \\
\mathbf{b}(0) &= \begin{pmatrix} 0.01398 & 0 & 0 \\ 0 & -0.00399 & 0 \\ 0 & 0 & -0.00998 \end{pmatrix}.
\end{aligned}
\qquad (61) $$

### 7.1.3. Initial condition for an Andrade rheology

The libration of Enceladus under Andrade rheology is determined by equations (5) and (18). The parameters of the rheology are those in Section 7.1.1 and the initial conditions are those in equation (61). The Andrade model also requires the initialisation of the internal degrees of freedom of the rheology, namely $\mathbf{b}_1, \mathbf{b}_2, \dots$. At time $t = 0$ they are all set to zero.
<table><caption>Table 3: Changing from MKS units to Astronomical unit, Mass of Sun and Year</caption><thead><tr><th>Parameter</th><th>Notation</th><th>Units</th><th>Values</th></tr></thead><tbody><tr><td>Astronomical unit</td><td>au</td><td>m</td><td>149597871000</td></tr><tr><td>Year</td><td>yr</td><td>s</td><td>31536000</td></tr><tr><td>Mass of Sun</td><td>M<sub>Sun</sub></td><td>kg</td><td>1.9885 × 10<sup>30</sup></td></tr></tbody></table>

<figure><img src="image_7.png" /><figcaption>Figure 6: Libration in longitude of Enceladus obtained by numerical integration of the equations of motion (5) and (18) associated with the extended Burgers model with 6 Voigt elements. The evolution is plotted in terms of angular acceleration ($\dot{\omega}/n^2$) versus normalised spin rate ($\omega - n/n$) after the decay of the transient state. The same figure obtained in the Maxwell case from the integration of Eqs. (5) and (14) is not represented because the difference between the two models is smaller than the thickness of the curve.</figcaption></figure>

## 7.2. Results

### 7.2.1. Forced libration

The rotation angle $\theta$ does not appear explicitly in the equations of motion in the inertial reference frame, (5), (14), and (18), though $\dot{\theta} = \Omega_3$ does. The main component of the forced libration $\zeta = \theta - M$ is oscillating at the orbital frequency $n$ (see Section 6) and once the free libration is totally damped the amplitude of forced libration is given by

$$ \beta = \sqrt{\left(\frac{\Omega_3}{n} - 1\right)^2 + \left(\frac{\dot{\Omega}_3}{n^2}\right)^2}. \qquad (62) $$

The two quantities in the square root, obtained by numerical integration of the equations of motion, are plotted against each other in Figure 6. The circular trajectory described in this set of coordinates justifies the expression (62). The associated radius $\beta$ is plotted in Figure 7 together with the analytical approximations derived in Section 6. As expected from Section 5.2.1, at the orbital frequency of Enceladus the amplitudes of libration are very close for both Maxwell and Andrade rheologies. We obtained a value of 0.000508 rad which is about one fourth of the observed amplitude. This apparent contradiction is discussed in Section 8. The results are in good agreement with the analytical approximation (55). Nevertheless, the numerical outputs show a small modulation of this amplitude not captured by the analytical approximation and a slight offset of 0.4% in both the Maxwell and the Andrade cases. But for Enceladus these differences have negligible consequences on the averaged dissipation.

### 7.2.2. Dissipation rate

With the numerical integrations, we have access to the instantaneous dissipation rate within the system given by the sum of the dissipation rates in all dash-pots of the oscillator representing the rheology. Results are plotted in
<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /><figcaption>Figure 7: Amplitude of libration in longitude after the decay of the transient state. The orange curve is obtained by numerical integration of the equations of motion (5) and (14) in the Maxwell case (top) or (5) and (18) in the Andrade case (bottom). The blue horizontal lines show the mean values of the integrations. The red horizontal lines represent the analytical solutions (55).</figcaption></figure>
<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /><figcaption>Figure 8: Dissipation rate (in Watts) as a function of time in orbital period once the steady state is reached. The orange curve is obtained by numerical integration of the equations of motion (5) and (14) in the Maxwell case (top) or (5) and (18) in the Andrade case (bottom). The blue horizontal lines show the mean values of integration. The red horizontal lines represent the analytical expression (58).</figcaption></figure>

Figure 8 over 10 orbital periods after the damping of the free libration. They display small oscillations with amplitude of 2 GW at twice the orbital frequency around a mean value equal to 9.80 GW and 10.50 GW in the Maxwell and in the extended Burgers model, respectively. The closeness of these two values is due to the fact that the imaginary part of the Love number is about the same in both rheologies at the orbital frequency (see Figure 9).

For comparison, the analytical dissipation rates given by equation (58) are respectively 9.81 GW and 10.51 GW. These values, represented as horizontal blue lines in Figure 8, are in very good agreement with the averaged numerical values.

Given the low amplitude of the forced libration in our model compared to the observed one (see previous section), their contribution to the total dissipation rate is only of the order of 7% instead of the 30% previously reported based on the observed libration (Efroimsky, 2018b; Folonier et al., 2018).

8. Observed vs Model libration: should they be the same?

The observed magnitude of libration of Enceladus is 0.0021 rad (see Thomas et al. (2016)). The value computed in this paper, with rheological parameters from the literature, is 0.000508 rad. ⁵Should these figures be the same? The

<footnote n="5">⁵In Thomas et al. (2016) (Table 1) the amplitude of libration $\beta$ of Enceladus was computed for four different interior models: - one homogeneous ellipsoidal rigid body in hydrostatic equilibrium, - two different two-layers rigid bodies in hydrostatic equilibrium, and - an external ice shell over a subsurface ocean over a solid core. The value of $\beta$ found for the first three models was approximately 0.033°. The value of $\beta$ found in our paper, for a deformable Enceladus with two different rheologies and no a priori assumed deformation, is approximately 0.03°. The agreement between our result and that in Thomas et al. (2016) is very good and the difference may be due to the fact that they used the rigid body formula as that in</footnote>
<figure><img src="image_12.png" /><figcaption>Figure 9: Imaginary part of the Love number as a function of $\omega$. The black vertical line marks the value of mean rotation rate of Enceladus.</figcaption></figure>

answer is no, because they measure different physical entities. In this paper the Enceladus' frame $K$ is a Tisserand frame, namely, it is an orthogonal moving frame with the origin at the centre of mass of Enceladus and with an angular velocity with respect to the inertial frame such that the angular momentum of the body with respect to $K$ is null (Munk and MacDonald, 1961). While in this paper we computed the libration of $K$, the observed magnitude of libration of Enceladus is the libration of its external shell that moves over a global subsurface sea (Thomas et al., 2016).

In order to relate the libration of $K$ to that of the shell we will assume, as done in Thomas et al. (2016), an interior model for Enceladus made of three layers. The external layer is a spherical shell 23 km thick with density 850 kg/m³. The second layer is a spherical ocean 28.1 km thick with density 1000 kg/m³. The third layer is a spherical core of density 2300 kg/m³. The mean radius of Enceladus is 252.1 km. Let $\dot{\theta}_1$, $\dot{\theta}_2$, and $\dot{\theta}_3$ denote the angular velocity (in the Tisserand sense) of each layer with respect to the inertial frame. If $I_1$, $I_2$, and $I_3$ denote the moment of inertia of each layer then $I_1 + I_2 + I_3 = I_0$ and the definition of Tisserand's frame implies

$$I_1(\dot{\theta}_1 - \dot{\theta}) + I_2(\dot{\theta}_2 - \dot{\theta}) + I_3(\dot{\theta}_3 - \dot{\theta}) = 0. \qquad (63)$$

Supposing that the liquid ocean does not librate, $\dot{\theta}_2 = n$, and that the angle of libration of the shell and core are denoted as $\zeta_1$ and $\zeta_3$, with $\dot{\theta}_1 = n + \zeta_1$ and $\dot{\theta}_3 = n + \zeta_3$, then equation (63) implies

$$I_1 \zeta_1 + I_3 \zeta_3 = I_0 \zeta. \qquad (64)$$

Assuming that all librations are in phase then this equation can be integrated to give

$$I_1 \zeta_1 + I_3 \zeta_3 = I_0 \zeta. \qquad (65)$$

The substitution of the values of $I_1$, $I_3$, $I_0$, $\zeta_1 = 0.0021$, and $\zeta = 0.000508$ into equation (65) gives an estimate $\zeta_3 = 0.006 \zeta_1$. So, assuming the interior model of Enceladus given in Thomas et al. (2016) then the results in this paper imply that the amplitude of the libration of the core is approximately 0.6% of that of the shell.

## 9. Conclusions

In this paper, we present a method to simulate the evolution in the time domain of an extended body subject to tides with arbitrary rheology. In particular, Andrade's creep model is known to be well suited for modelling tidal deformation of icy satellites. But it involves a fractional power law in the frequency domain preventing a straightforward transposition into a set of ordinary differential equations in the time domain. To circumvent this issue, Andrade's

<footnote>equation (56) while we used (48) that takes into account the rheology of the satellite. The disagreement of all these values of $\beta$ with the observed amplitude of libration of Enceladus' surface, 0.120°, was the main evidence found in Thomas et al. (2016) for the existence of a subsurface ocean. Within the context of our model we present another, but equivalent, explanation for this difference, see text.</footnote>
compliance function is approximated by that of an extended Burgers model containing an arbitrary number of Voigt elements.

The procedure is applied to Saturn's satellite Enceladus for which the forced libration in longitude has been shown to be responsible for a significant amount of energy dissipation. For this particular body, we consider both the Maxwell and the Andrade visco-elastic models. Tidal effects on the orbit are switched off in order to maintain a constant eccentricity which, in the special case of Enceladus, is driven by the 2:1 mean motion resonance with Dione.

In this particular setting, and for arbitrary rheology represented by a Love number $k_2$, we show that the amplitude of libration $\beta$ is given by

$$ \beta \approx 6e^{\frac{\bar{B}-\bar{A}}{\bar{C}}} \frac{m_2}{m_1+m_2} \left| 1 - \frac{5k_2(n)}{6k_2(0)} \right| $$

This equation predicts values reduced by a factor $|1 - 5k_2(n)/6k_2(0)|$ with respect to the classical formula obtained in the case of a strictly rigid body with the same permanent triaxiality. We also recover the expression of the average rate of tidal dissipation of a synchronous body valid for any rheology, namely (Efroimsky, 2018a)

$$ \frac{\Delta E}{T} = -\Im[k_2(n)] \frac{(nR)^5}{G} \left( \frac{m_2}{m_1+m_2} \right)^2 \times \left[ \frac{21}{2} e^2 + 6e\beta + \frac{3}{2}\beta^2 \right] $$

Thanks to these two formulae, we have been able to validate the numerical procedure. The analytical estimates of the libration amplitude (0.0005 rad) and of the dissipation rate (10 GW) differ by less than 0.4% from the averaged values obtained by numerical integration. In our model, the libration in longitude is that of the Tisserand frame. The comparison with the observed libration of the shell allowed us to estimate that the core librates with an amplitude 0.6% times smaller than the shell.

In addition to the averaged libration amplitude and dissipation rate, the numerical integration also displays oscillations not captured by the first order analytical calculation. Moreover, the procedure can be applied in more complex situations where a Fourier series is not feasible, e.g., in systems with large eccentricities and/or with non quasiperiodic perturbations as in a chaotic system with close encounters.

## Acknowledgments

We would like to thank Michael Efroimsky and another anonymous referee for careful reading and many suggestions improving the manuscript. CR is partially supported by FAPESP grant 2016/25053-8. YG is partially supported by FAPESP grants 2015/26253-8 and 2018/02905-4. AC is supported from CFISUC strategic project (UID/FIS/04564/2019), ENGAGE SKA (POCI-01-0145-FEDER-022217), and PHOBOS (POCI-01-0145-FEDER-029932), funded by COMPETE 2020 and FCT, Portugal.

## References

D. L. Anderson and J. B. Minster. The frequency dependence of Q in the Earth and implications for mantle rheology and Chandler wobble. *Geophysical Journal*, 58:431–440, Aug 1979. doi: 10.1111/j.1365-246X.1979.tb01033.x.

E. N. Andrade. On the viscous flow in metals, and allied phenomena. *Proceedings of the Royal Society of London*, 84:1–12, 1910.

Erik Asphaug and Andreas Reufer. Late origin of the Saturn system. *Icarus*, 223(1):544–565, Mar 2013. doi: 10.1016/j.icarus.2012.12.009.

Abderrahmin Ben Jazia, Bruno Lombard, and Cédric Bellis. Wave propagation in a fractional viscoelastic Andrade medium: diffusive approximation and numerical modeling. *arXiv e-prints*, art. arXiv:1312.4820, Dec 2013.

C. J. Bierson and F. Nimmo. A test for Io's magma ocean: Modeling tidal dissipation with a partially molten mantle. *Journal of Geophysical Research (Planets)*, 121(11):2211–2224, Nov 2016. doi: 10.1002/2016JE005005.

Carolin Birk and Chongmin Song. An improved non-classical method for the solution of fractional differential equations. *Computational Mechanics*, 46(5):721–734, Oct 2010. doi: 10.1007/s00466-010-0510-4.

D. R. Bland. *The Theory of Linear Viscoelasticity*. Dover Books on Physics. Dover Publications, 2016. ISBN 9780486816388. URL https://books.google.com.br/books?id=YL8zDQAAQBAJ.

Gwenaël Boué, Alexandre C. M. Correia, and Jacques Laskar. Complete spin and orbital evolution of close-in bodies using a Maxwell viscoelastic rheology. *Celestial Mechanics and Dynamical Astronomy*, 126(1-3):31–60, Nov 2016. doi: 10.1007/s10569-016-9708-x.
Marie Běhounková, Gabriel Tobie, Gaël Choblet, and Ondřej Čadek. Impact of tidal heating on the onset of convection in Enceladus's ice shell. *Icarus*, 226(1):898-904, Sep 2013. doi: 10.1016/j.icarus.2013.06.033.

Marie Běhounková, Gabriel Tobie, Ondřej Čadek, Gaël Choblet, Carolyn Porco, and Francis Nimmo. Timing of water plume eruptions on Enceladus explained by interior viscosity structure. *Nature Geoscience*, 8:601–604, Aug 2015. doi: 10.1038/ngeo2475.

Robin M. Canup and William R. Ward. A common mass scaling for satellite systems of gaseous planets. *Nature*, 441(7095):834–839, Jun 2006. doi: 10.1038/nature04860.

Julie C. Castillo-Rogez, Michael Efroimsky, and Valéry Lainey. The tidal history of Iapetus: Spin dynamics in the light of a refined dissipation model. *Journal of Geophysical Research (Planets)*, 116(E9):E09008, Sep 2011. doi: 10.1029/2010JE003664.

Sébastien Charnoz, Aurélien Crida, Julie C. Castillo-Rogez, Valery Lainey, Luke Dones, Özgür Karatekin, Gabriel Tobie, Stephane Mathis, Christophe Le Poncin-Lafitte, and Julien Salmon. Accretion of Saturn's mid-sized moons during the viscous spreading of young massive rings: Solving the paradox of silicate-poor rings versus silicate-rich moons. *Icarus*, 216(2):535–550, Dec 2011. doi: 10.1016/j.icarus.2011.09.017.

A. C. M. Correia, C. Ragazzo, and L. S. Ruiz. The effects of deformation inertia (kinetic energy) in the orbital and spin evolution of close-in bodies. *Celestial Mechanics and Dynamical Astronomy*, 130(8):51, Aug 2018. doi: 10.1007/s10569-018-9847-3.

Alexandre C. M. Correia, Gwenaël Boué, Jacques Laskar, and Adrián Rodríguez. Deformation and tidal evolution of close-in planets and satellites using a Maxwell viscoelastic rheology. *Astronomy & Astrophysics*, 571:A50, Nov 2014. doi: 10.1051/0004-6361/201424211.

J. M. A. Danby. *Fundamentals of Celestial Mechanics*. Macmillan, 1962.

Michael Efroimsky. Tidal Dissipation Compared to Seismic Dissipation: In Small Bodies, Earths, and Super-Earths. *Astrophysical Journal*, 746 (2):150, Feb 2012a. doi: 10.1088/0004-637X/746/2/150.

Michael Efroimsky. Bodily tides near spin-orbit resonances. *Celestial Mechanics and Dynamical Astronomy*, 112(3):283–330, Mar 2012b. doi: 10.1007/s10569-011-9397-4.

Michael Efroimsky. Tidal Evolution of Asteroidal Binaries. Ruled by Viscosity. *Ignorant of Rigidity*. *Astronomical Journal*, 150(4):98, Oct 2015. doi: 10.1088/0004-6256/150/4/98.

Michael Efroimsky. Dissipation in a tidally perturbed body librating in longitude. *Icarus*, 306:328–354, May 2018a. doi: 10.1016/j.icarus.2017.10.020.

Michael Efroimsky. Tidal viscosity of Enceladus. *Icarus*, 300:223–226, Jan 2018b. doi: 10.1016/j.icarus.2017.09.013.

Michael Efroimsky and Valéry Lainey. Physics of bodily tides in terrestrial planets and the appropriate scales of dynamical evolution. *Journal of Geophysical Research (Planets)*, 112(E12):E12003, Dec 2007. doi: 10.1029/2007JE002908.

Sylvio Ferraz-Mello. Tidal synchronization of close-in satellites and exoplanets. A rheophysical approach. *Celestial Mechanics and Dynamical Astronomy*, 116(2):109–140, Jun 2013. doi: 10.1007/s10569-013-9482-y.

Sylvio Ferraz-Mello. The small and large lags of the elastic and anelastic tides. The virtual identity of two rheophysical theories. *Astronomy & Astrophysics*, 579:A97, Jul 2015a. doi: 10.1051/0004-6361/201525900.

Sylvio Ferraz-Mello. Tidal synchronization of close-in satellites and exoplanets: II. Spin dynamics and extension to Mercury and exoplanet host stars. *Celestial Mechanics and Dynamical Astronomy*, 122(4):359–389, Aug 2015b. doi: 10.1007/s10569-015-9624-5.

H. A. Folonier, S. Ferraz-Mello, and E. Andrade-Ines. Tidal synchronization of close-in satellites and exoplanets. III. Tidal dissipation revisited and application to Enceladus. *Celestial Mechanics and Dynamical Astronomy*, 130(12):78, Dec 2018. doi: 10.1007/s10569-018-9872-2.

J. W. Glen. The creep of polycrystalline ice. *Proceedings of the Royal Society of London*, 228(1175):519–38, 1955.

C. Goetze and W. F. Brace. Laboratory observations of high-temperature rheology of rocks. *Tectonophysics*, 13(1):583–600, Apr 1972. doi: 10.1016/0040-1951(72)90039-X.

Christopher Goetze. High temperature rheology of westerly granite. *Journal of Geophysical Research*, 76(5):1223–1230, Feb 1971. doi: 10.1029/JB076i005p01223.

D. T. Griggs, F. J. Turner, and H. C. Heard. Chapter 4: Deformation of Rocks at 500Å° to 800Å° C. In *Rock Deformation (A Symposium)*. Geological Society of America, 03 1960. ISBN 9780813710792. doi: 10.1130/MEM79-p39. URL https://doi.org/10.1130/MEM79-p39.

E. Hairer, S. P. Nørsett, and Wanner G. *Solving Ordinary Differential Equations I*. Springer, 1993.

H. Jeffreys. Creep in the earth and planets* *read by e.r. lapwood. In A.R. RITSEMA, editor, *The Upper Mantle*, volume 4 of *Developments in Geotectonics*, pages 569 – 581. Elsevier, 1972. doi: https://doi.org/10.1016/B978-0-444-41015-3.50035-5. URL http://www.sciencedirect.com/science/article/pii/B9780444410153500355.

Shunichi Kamata and Francis Nimmo. Interior thermal state of Enceladus inferred from the viscoelastic state of the ice shell. *Icarus*, 284:387–393, Mar 2017. doi: 10.1016/j.icarus.2016.11.034.

S. Karato and H. A. Spetzler. Defect microdynamics in minerals and solid state mechanisms of seismic wave attenuation and velocity dispersion in the mantle. *Reviews of Geophysics*, 28(4):399–421, Nov 1990. doi: 10.1029/RG028i004p00399.

Jurriën Sebastiaan Knibbe and Wim van Westrenen. On Mercury's past rotation, in light of its large craters. *Icarus*, 281:1–18, Jan 2017. doi: 10.1016/j.icarus.2016.08.036.

Harriet C. P. Lau and Ulrich H. Faul. Anelasticity from seismic to tidal timescales: Theory and observations. *Earth and Planetary Science Letters*, 508:18–29, Feb 2019. doi: 10.1016/j.epsl.2018.12.009.

Jérémy Leconte, Hanbo Wu, Kristen Menou, and Norman Murray. Asynchronous rotation of Earth-mass planets in the habitable zone of lower-mass stars. *Science*, 347(6222):632–635, Feb 2015. doi: 10.1126/science.1258686.

Valeri V. Makarov. Equilibrium Rotation of Semiliquid Exoplanets and Satellites. *Astrophysical Journal*, 810(1):12, Sep 2015. doi: 10.1088/0004-637X/810/1/12.

Valeri V. Makarov and Michael Efroimsky. No Pseudosynchronous Rotation for Terrestrial Planets and Moons. *Astrophysical Journal*, 764(1):27, Feb 2013. doi: 10.1088/0004-637X/764/1/27.

Valeri V. Makarov, Ciprian Berghea, and Michael Efroimsky. Dynamical Evolution and Spin-Orbit Resonances of Potentially Habitable Exoplanets: The Case of GJ 581d. *Astrophysical Journal*, 761(2):83, Dec 2012. doi: 10.1088/0004-637X/761/2/83.

Valeri V. Makarov, Julien Frouard, and Bryan Dorland. Forced libration of tidally synchronized planets and moons. *Monthly Notices of the Royal Astronomical Society*, 456(1):665–671, Feb 2016. doi: 10.1093/mnras/stv2735.

Christine McCarthy and Julie C. Castillo-Rogez. Planetary Ices Attenuation Properties. In Murthy S. Gudipati and Julie Castillo-Rogez, editors,
*Astrophysics and Space Science Library*, volume 356, page 183, Jan 2013. doi: 10.1007/978-1-4614-3076-6_7.

N.F. Mott. Lxxviii. a theory of work-hardening of metals ii: Flow without slip-lines, recovery and creep. *The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science*, 44(354):742–765, 1953. doi: 10.1080/14786440708521052.

W. H. Munk and G. J. F. MacDonald. *The rotation of the Earth*. Cambridge University Press, New York, 1961.

F. Nimmo, U. H. Faul, and E. J. Garnero. Dissipation at tidal and seismic frequencies in a melt-free Moon. *Journal of Geophysical Research (Planets)*, 117(E9):E09005, Sep 2012. doi: 10.1029/2012JE004160.

Benoît Noyelles. Interpreting the librations of a synchronous satellite - How their phase assesses Mimas' global ocean. *Icarus*, 282:276–289, Jan 2017. doi: 10.1016/j.icarus.2016.10.001.

Benoît Noyelles, Julien Frouard, Valeri V. Makarov, and Michael Efroimsky. Spin-orbit evolution of Mercury revisited. *Icarus*, 241:26–44, Oct 2014. doi: 10.1016/j.icarus.2014.05.045.

Sebastiano Padovan, Jean-Luc Margot, Steven A. Hauck, William B. Moore, and Sean C. Solomon. The tides of Mercury and possible implications for its interior structure. *Journal of Geophysical Research (Planets)*, 119(4):850–866, Apr 2014. doi: 10.1002/2013JE004459.

Gérard Petit and Brian Luzum. IERS Conventions (2010). IERS Technical Note No 36, Frankfurt am Main: Verlag des Bundesamts für Kartographie und Geodäsie, 2010. 179 pp., ISBN 3-89888-989-6.

C. C. Porco, P. Helfenstein, P. C. Thomas, A. P. Ingersoll, J. Wisdom, R. West, G. Neukum, T. Denk, R. Wagner, and T. Roatsch. Cassini Observes the Active South Pole of Enceladus. *Science*, 311(5766):1393–1401, Mar 2006. doi: 10.1126/science.1123013.

C. Ragazzo and L. S. Ruiz. Viscoelastic tides: models for use in Celestial Mechanics. *Celestial Mechanics and Dynamical Astronomy*, 128(1): 19–59, May 2017. doi: 10.1007/s10569-016-9741-9.

Clodoaldo Ragazzo. The theory of figures of Clairaut with focus on the gravitational rigidity modulus: inequalities and an improvement in the Darwin-Radau equation. *arXiv preprint arXiv:1811.07759*, 2018.

Nicolas Rambaux, Julie C. Castillo-Rogez, James G. Williams, and Özgür Karatekin. Librational response of Enceladus. *Geophysical Research Letters*, 37(4):L04202, Feb 2010. doi: 10.1029/2009GL041465.

Joe P. Renaud and Wade G. Henning. Increased Tidal Dissipation Using Advanced Rheological Models: Implications for Io and Tidally Active Exoplanets. *Astrophysical Journal*, 857(2):98, Apr 2018. doi: 10.3847/1538-4357/aab784.

Ignasi Ribas, Emeline Bolmont, Franck Selsis, Ansgar Reiners, Jérémy Leconte, Sean N. Raymond, Scott G. Engle, Edward F. Guinan, Julien Morin, and Martin Turbet. The habitability of Proxima Centauri b. I. Irradiation, rotation and volatile inventory from formation to the present. *Astronomy and Astrophysics*, 596:A111, Dec 2016. doi: 10.1051/0004-6361/201629576.

J. Salmon and R. M. Canup. Accretion of Saturn's Inner Mid-sized Moons from a Massive Primordial Ice Ring. *The Astrophysical Journal*, 836 (1):109, Feb 2017. doi: 10.3847/1538-4357/836/1/109.

M. Segatz, T. Spohn, M. N. Ross, and G. Schubert. Tidal dissipation, surface heat flow, and figure of viscoelastic models of Io. *Icarus*, 75(2): 187–206, Aug 1988. doi: 10.1016/0019-1035(88)90001-2.

D. Shoji, H. Hussmann, K. Kurita, and F. Sohl. Ice rheology and tidal heating of Enceladus. *Icarus*, 226(1):10–19, Sep 2013. doi: 10.1016/j.icarus. 2013.05.004.

Ondřej Souček, Marie Běhounková, Ondřej Čadek, Jaroslav Hron, Gabriel Tobie, and Gaël Choblet. Tidal dissipation in Enceladus' uneven, fractured ice shell. *Icarus*, 328:218–231, Aug 2019. doi: 10.1016/j.icarus.2019.02.012.

F. D. Stacey. The theory of creep in rocks and the problem of convection in the Earth's mantle. *Icarus*, 1(1-6):304–312, Jan 1963. doi: 10.1016/ 0019-1035(62)90028-3.

P. C. Thomas, R. Tajeddine, M. S. Tiscareno, J. A. Burns, J. Joseph, T. J. Loredo, P. Helfenstein, and C. Porco. Enceladus's measured physical libration requires a global subsurface ocean. *Icarus*, 264:37–47, Jan 2016. doi: 10.1016/j.icarus.2015.08.037.

T. Van Hoolst, R. Baland, and A. Trinh. On the librations and tides of large icy satellites. *Icarus*, 226(1):299–315, 2013. doi: https://doi.org/10. 1016/j.icarus.2013.05.036. URL http://www.sciencedirect.com/science/article/pii/S0019103513002364.

M. Walterová and M. Běhounková. Tidal effects in differentiated viscoelastic bodies: a numerical approach. *Celestial Mechanics and Dynamical Astronomy*, 129(1-2):235–256, Sep 2017. doi: 10.1007/s10569-017-9772-x.

James G. Williams and Dale. H. Boggs. Tides on the Moon: Theory and determination of dissipation. *Journal of Geophysical Research (Planets)*, 120(4):689–724, Apr 2015. doi: 10.1002/2014JE004755.

James G. Williams, Alexander S. Konopliv, Dale H. Boggs, Ryan S. Park, Dah-Ning Yuan, Frank G. Lemoine, Sander Goossens, Erwan Mazarico, Francis Nimmo, and Renee C. Weber. Lunar interior properties from the GRAIL mission. *Journal of Geophysical Research (Planets)*, 119(7): 1546–1578, Jul 2014. doi: 10.1002/2013JE004559.

C. F. Yoder. *Astrometric and geodetic properties of Earth and the Solar System*, volume 1. Wiley Online Library, 1995.