The effects of deformation inertia (kinetic energy) in the orbital and spin evolution of close-in bodies

The purpose of this work is to evaluate the effect of deformation inertia on tide dynamics, particularly within the context of the tide response equations proposed independently by Boué et al. (Celest Mech Dyn Astron 126:31–60, 2016) and Ragazzo and Ruiz (Celest Mech Dyn Astron 128(1):19–59, 2017). The singular limit as the inertia tends to zero is analyzed, and equations for the small inertia regime are proposed. The analysis of Love numbers shows that, independently of the rheology, deformation inertia can be neglected if the tide-forcing frequency is much smaller than the frequency of small oscillations of an ideal body made of a perfect (inviscid) fluid with the same inertial and gravitational properties of the original body. Finally, numerical integration of the full set of equations, which couples tide, spin and orbit, is used to evaluate the effect of inertia on the overall motion. The results are consistent with those obtained from the Love number analysis. The conclusion is that, from the point of view of orbital evolution of celestial bodies, deformation inertia can be safely neglected. (Exceptions may occur when a higher-order harmonic of the tide forcing has a high amplitude.)


Introduction
The current theories of tides are based on the fundamental works of Newton, Laplace, Kelvin and Darwin; different theories being distinguished by the way extended bodies deform under a given tide-raising force. Following George Darwin, the usual approach to the deformation of a body is: to spatially decompose the tide-raising potential into spherical harmonics (as was first done by Laplace), to time decompose the tide-raising potential into Fourier modes, to apply the same decomposition to the gravitational field of the deformed body, and then to postulate a linear relation between corresponding tide-raising and tide-response terms; the linear factor being the "Love number" (or "real Love number"). Dissipation is introduced by means of a time lag between the instantaneous force and the delayed tide response. In principle, different time delays can be chosen for different space and time harmonic modes. The complex representation of each Fourier mode e iωt implies that the time delay τ can be represented by a complex multiplicative factor e −iωτ . The multiplication of the real Love number (amplitude factor) by the complex phase factor gives rise to the so-called complex Love number. The real Love number can be interpreted as a measure of the rigidity of the body and the time delay as a measure of the viscosity of the body. The question to be investigated in this paper is: Does "tide deformation inertia," or just "deformation inertia," have any effect on the Love numbers? By "tide deformation inertia," we mean the resistance of the extended body to any change in its motion of tide deformation. Clearly, deformation inertia is significant only when tide forcing undergoes a "sufficiently fast" time variation. So, the above question can be rephrased as: How fast must be the variation of the tide forcing for deformation inertia to have an effect on the Love numbers?
There are observational data showing that inertia plays an important role in the Earth ocean tides. The Love numbers k of the Earth can be split into two parts, k = k S + k O , the first due to the deformation of the solid part k S and the second due to the deformation of the oceans k O [see for instance Ray et al. (2001) or Petit and Luzum (2010, section 6)]. For the most important Love numbers of the Earth, the real part of k S are positive, while the real part of k O are negative. In Ray et al. (2001), for instance, for the principal lunar semi-diurnal tide M 2 , k S and k O are written as [see equations (2) to (7) in Ray et al. (2001)] H D e −iψ = |k O |e −iψ (equation (9) in Ray et al. (2001)) where k 2 = 0.302 is the solid Earth elastic Love number [ Table 2 in Ray et al. (2001)], α(1+k ) H = 0.9590 m −1 [equation (9) in Ray et al. (2001)], D = 3.295 cm is an ocean amplitude and ψ = 128.69 • [D and ψ are given in Equation (20) in Ray et al. (2001)]. So, for the principal lunar semi-diurnal tide M 2 , k S = 0.302 [the imaginary part of k S was estimated as 10 −1 in Ray et al. (2001) and has been neglected here] k O = −0.020 − i0.025 and k = 0.282 − i0.025 [this is essentially the same value obtained from Petit and Luzum (2010) and reported in Table 1 in Ragazzo and Ruiz (2017)]. The negativeness of the real part of a Love number shows that the tide kinetic energy associated with this Love number, which is proportional to the deformation inertia, is not only relevant but dominant over the Fig. 1 Amplitude variation of the steady-state solution to the equation μẍ+ ηẋ +γ x = cos(ωt) as a function of ω/ω 0 , where ω 0 = √ γ /μ is the undamped frequency of oscillations and ξ = ηω 0 /(2γ ) is the damping ratio tide potential energy (self-gravity and elastic), see 1 Ragazzo and Ruiz (2017, Appendix 2). So, deformation inertia has an important effect at least on the Earth ocean tides.
The study of the effect of inertia on the solid part of the Earth seems to be pioneered by Love. Love (1911, chapter 5) estimated the correction due to inertia to the "purely statical problem of a homogeneous incompressible sphere held strained by the forces derived from a tide generating potential" (Love 1911, paragraph 67). If the shear modulus of the Earth is neglected (μ = 0 according to Love's notation), then from equations (13), (20) and (49) in chapter 5 of Love (1911) we obtain that deformation inertia accounts for 0.5% of the Earth's response (or 0.5% of the Love number k 2 ) to the principal lunar semi-diurnal tide forcing. So, in this case, deformation inertia is negligible. Love states: If the tide-forcing frequency is much smaller than the frequency of free oscillations of the deformable body, then the effect of inertia on tides can be neglected.
The question of the effect of deformation inertia on tides would be fully settled by Love's result if viscosity would not exist. As shown in the example in Fig. 1, energy dissipation shifts the resonance frequency and sets its amplitude and width. In principle, it could be that some special rheology could decrease the resonance frequency of the deformable body in such a way that deformation inertia would become relevant.
From a theoretical point of view, tide models parameterized by Love numbers have a drawback: Infinitely many complex parameters, the Love numbers, must be chosen. One way to overcome this difficulty is to assume, as Kelvin, Love and others, that the extended body has a very simple structure as, for instance, being a homogeneous fluid or solid. Since the Love numbers encapsulates the structure of the extended body the simplifying assumption may lead to controversial results.
There are several simplified tide models which are not explicitly parameterized by Love numbers [see, for instance, Mignard (1979), Efroimsky andWilliams (2009), Zlenko (2014), Celletti (1990), Antognini et al. (2014), Bambusi and Haus (2015), Ferraz-Mello (2013), Correia et al. (2014), Boué et al. (2016), Wisdom and Meyer (2016), Ragazzo and Ruiz (2015), Ragazzo and Ruiz (2017)], two of them will be important in this paper: the first 1 In order to understand this claim, it is enough to analyze the one degree of freedom harmonic oscillator μẍ + ηẋ +γ x = cos(ωt), with solution x(t) = a cos(ωt)−b sin(ωt), where a+ib plays the role of the Love number. Multiplying both sides of the equation by x and time-averaging gives lim T →∞ dt that is the time-average balance of kinetic and potential energy. This same reasoning applied to the linear equations of tide dynamics, equation (47) in Ragazzo and Ruiz (2017), shows that the time-average tidal balance of energy is proportional to the real part of the Love number. presented in Boué et al. (2016) after previous work in Ferraz-Mello (2013) and Correia et al. (2014) and the second presented in Ragazzo and Ruiz (2017) after previous work in Ragazzo and Ruiz (2015).
Both works Boué et al. (2016) and Ragazzo and Ruiz (2017) present three-dimensional equations for the spin and orbital motion of linear viscoelastic bodies interacting under gravity, the main differences being: (1) in Boué et al. (2016) the full multipole expansion of a gravitational field can be taken into account, while in Ragazzo and Ruiz (2017) the expansion is truncated to quadrupole order; (2) in Boué et al. (2016) a linear viscoelastic rheology is introduced by means of a Correspondence Principle [see section 4.3 and Appendix B of Efroimsky (2012)], while in Ragazzo and Ruiz (2017) it is introduced by means of an "Association Principle," to be explained in the following sections; and (3) in Boué et al. (2016) no inertia is associated with tidal deformations, while the opposite is true in Ragazzo and Ruiz (2017).
As in the examples in Correia et al. (2014) and Boué et al. (2016), in the present paper, the multipole expansion of the gravitational field of a deformed body is truncated to the quadrupole order and only the Maxwell and the Voigt rheologies are considered in the timedomain simulations. Under these conditions, the only difference between the equations in Boué et al. (2016) and those in Ragazzo and Ruiz (2017) is the presence of an inertia parameter, μ ≥ 0, associated with tide deformations in the latter. Therefore, in this paper we follow the presentation of the equations of motion as in Ragazzo and Ruiz (2017) that allows for the recovering of the equations in Boué et al. (2016) simply by making μ = 0.
The paper is organized as follows. In Sect. 2, the mathematical model to be used in the paper is presented and some of the parameters in the model are estimated. This section is a summary of part of Ragazzo and Ruiz (2017). In Sect. 3, nondimensional variables are introduced for the particular cases of the Kelvin-Voigt and Maxwell rheologies. A small inertia parameter induces high-frequency oscillations on the deformable bodies, and this slows down considerably the numerical simulations. For this reason, in Sect. 3 and Appendix A we present an approximation to the equations of motion that smooths out these oscillations by means of an averaging. This is the first new result in the paper. In Sect. 4, we present the expressions for the Love numbers for a quite broad class of rheologies and compare the usual "Correspondence Principle" to the "Association Principle" derived in Ragazzo and Ruiz (2017) and used in this paper; both lead to the same Love numbers. Section 5 contains the main contribution of this paper, a quantitative upper bound for the difference between the Love numbers in the case μ > 0 and μ = 0. This estimate is valid for the broad class of rheologies of Sect. 4. Section 6 contains the time-domain simulations where the systems with and without deformation inertia are compared; only the Maxwell and the Kelvin-Voigt rheologies are considered. The average secular torque with and without inertia is computed and compared. Section 7 is a conclusion where the main results in the paper are summarized. In this section, notwithstanding the negligible effect of deformation inertia in the orbital evolution of celestial bodies, we present a potential application of tide models with deformation inertia to obtain information on the rheology of planets at high frequencies.
(a) When the body is at rest (no rotational motion), its distribution of mass is spherically symmetric. (b) When the body has a rotational motion, its distribution of mass is almost spherically symmetric in the sense that the level sets of the density function are approximately ellipsoidal shells of small eccentricities. (c) The body material has an incompressible behavior under small deformations. (d) The body internal forces are such that the motion preserves the total angular momentum. (e) The mechanical deformation response of the body is linear viscoelastic.
For systems of deformable bodies, we must include the hypothesis: (f) The minimum distance between two bodies is sufficiently large such that the almost sphericity hypothesis (b) still holds.
Let κ be an inertial reference frame and consider a system of N deformable bodies. The kinematics of each body, labeled by an index i, is characterized by: the position x i of the center of mass with respect to κ; an orientation Y i : K i → κ, where Y i is a rotation matrix and K i is a Tisserand frame (defined below); and the traceless part of the moment of inertia operator of the body with respect to the frame K i , represented by a 3 × 3 symmetric traceless matrix B i . The choices of the Tisserand frame and the traceless part of the moment of inertia require explanations.
One of the main issues in the dynamics of deformable bodies is the choice of a "body frame," namely the choice of a moving reference frame that is either unambiguously defined (the principal axis of inertia), or in which the body is as at rest as possible (the Tisserand's mean axis). A Tisserand's frame K is an orthogonal moving frame with the origin at the center of mass of the body 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 (see Munk and MacDonald 1960). For a rigid body, a moving frame K is a Tisserand frame if, and only if, the body remains at rest with respect to K. A Tisserand frame is characterized by its angular velocity in the following way. Let l ∈ κ be the angular momentum vector of the body with respect to the inertial frame and L = Y −1 l ∈ K be its representation in a moving frame K. Let I : K → K denote the instantaneous moment of inertia operator of the body. K is a Tisserand frame if, and only if, L = IΩ, where Ω ∈ K is the angular velocity of the frame K with respect to κ. The angular velocity operator Ω : K → κ associated with the motion of K is given by whereẎ is the time derivative of Y. If Ω(t) is known and Y(0) is chosen, then Y(t) is uniquely determined by the integration ofẎ = YΩ. Let I : K → K denote the moment of inertia of a body with respect to a Tisserand frame (the following arguments are true with respect to any reference frame) and I • be the trace of I over 3. If a deformable body is incompressible, then I • does not change under small deformations (as shown by Darwin, see Rochester and Smylie 1974). So, the time variation of the moment of inertia is determined by the time variation of its traceless nondimensional part B defined by where I is the identity matrix. The following well-known fact is crucial to the model: The gravitational moment of quadrupole of the body Q : K → K is given by Let r denote the distance to the origin of K. Since the frame K is centered at the center of mass of the body, the gravitational field of the body with respect to K is given by the monopole part of order 1/r , plus the quadrupole pole part of order 1/r 3 , plus terms of order 1/r 4 that will be neglected. So, within the approximation considered in this paper, the gravitational field of a body is determined by the mass and I • , which are constants, and by the time-varying matrix B. The small distortion hypothesis (b) and I = I • (I − B) imply that B 1.
The differential equations for the position and the orientation of each body as a function of B andḂ can be derived from the usual laws of mechanics. The derivation of the equations for the deformation variables B is more difficult and requires hypotheses on the rheology of each body. The small distortion hypotheses (b) and (f) suggest that the equations for B must be linear about the equilibrium B = 0 that must exist when the system is at rest, due to hypothesis (a). We further suppose that B has a linear viscoelastic behavior, hypothesis (e). These hypotheses lead to differential equations for B by means of an "Association Principle" derived in Ragazzo and Ruiz (2017) and explained in the following.
The viscoelastic response of a linear material (or its force-extension equations) is analogous to the mechanical response of a one-dimensional spring-dashpot system (Bland 1960). The two simplest mechanical systems that model the viscoelastic behavior of a material are the Maxwell element and the Voigt (or Kelvin-Voigt) element shown in Fig. 2. From these two elements, more complex models can be constructed as those shown in Figs. 3 and 4 (see also Bland 1960). In most of this paper, we only consider rheologies that are modeled by the Maxwell and the Kelvin-Voigt elements. Notice that a Maxwell element experiences an unbounded deformation under a constant force, which seems unreasonable for setting an analogy between B and the spring-dashpot deformation. This apparent difficulty is overcome x by the fact that B is additionally under the effect of self-gravity that tends to restore any offset from the spherical shape. Due to the small deformation hypothesis, the self-gravitational force can also be linearized about the spherical shape represented by B = 0. Besides selfgravitation, there must be an inertia associated with time variations of the deformation of the body and therefore to the time variations of B, this is what we called deformation inertia in the introduction. Both effects, self-gravitation and deformation inertia, can be represented in the spring-dashpot system by an additional spring, with elastic constant γ , in parallel with the rheological elements and by an additional mass element μ that goes along with the deformation. This construction gives rise to the oscillators shown in Figs. 5 and 6. The equation of motion associated with the Kelvin oscillator is where γ is the elastic constant due to self-gravity plus the elastic constant of the spring in the Kelvin-Voigt element, η is the viscosity of the dashpot, and F(t) is the external force. The equations of motion associated with the Maxwell oscillator are where λ is the force acting upon the Maxwell element and α is the elastic constant of the Maxwell element. Using the constraint x = x 1 + x 2 , these equations can be written as The Maxwell oscillator formulated from the Maxwell element in Fig. 2. Remark that the Kelvin oscillator in Fig. 5 can be understood as this Maxwell oscillator in the limit as α → ∞ A given viscoelastic rheology is always associated with a spring-dashpot oscillator for a linear displacement x. The Association Principle in Ragazzo and Ruiz (2017) states that: "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) The justification of this principle relies upon: a Lagrangian formulation of the problem, the isotropy of space, and the incompressibility hypothesis (c) (see Ragazzo and Ruiz 2017, section 3.2). The comparison between this Association Principle and the mostly used Correspondence Principle is presented in Sect. 4.
In this paper, we are only interested in the system formed by a deformable body plus a point mass. For this system, let κ be an inertial reference frame at the center of mass of the system that is supposed at rest. Let x 1 ∈ κ and x 2 ∈ κ be the positions and m 1 and m 2 be the masses of the deformable body and the point mass, respectively. For the deformable body, let K be a Tisserand frame, Y : K → κ be the orientation matrix, Ω = Y −1Ẏ : K → K be the angular velocity operator, and B : K → K be the nondimensional traceless part of the moment of inertia operator, as given in Eq. (1). Let q = Y −1 (x 2 − x 1 ) be the relative position of the point mass with respect to the rotating reference frame K. For this system, the force F that acts upon B and that replaces the forces F(t) in the Kelvin oscillator (2) and in the Maxwell oscillator (3) is given by (see Ragazzo and Ruiz 2017 for details) where q ⊗ q denotes the matrix with entries (q ⊗ q) i j = q i q j . The equations for the orbital motion and the spin of the extended body arė where M = m 1 + m 2 and [B, Ω] = BΩ − ΩB denotes the commutator of matrices. Finally, due to the Association Principle, in the case of the Kelvin-Voigt rheology the equation for B is given by and in the case of the Maxwell rheology the equation for B is given by where F is given in Eq. (5), μ is a nondimensional parameter representing the inertia of B, γ (dimension sec −2 ) and α (dimension sec −2 ) are stiffness coefficients, η (dimension sec −1 ) is a damping coefficient, and Λ is a symmetric matrix (dimension sec −2 ) that represents the stress acting upon the Maxwell element, see Ragazzo and Ruiz (2017 for details). We recall that for μ = 0 and for the Maxwell rheology the above equations coincide with those presented in Boué et al. (2016), after a redefinition of parameters (see Ragazzo and Ruiz 2017, section 8 footnote 8).

The coefficient and an upper bound for the spin rate
The stiffness coefficient γ has an unusual dimension of sec −2 , which can be explained as follows. (The same type of argument explains the unusual dimensions of the other constants.) If the extended body is isolated and has constant spin Ω, then Eqs. (7) and (8) have a steadystate solution given by If I i j (spin) denotes the change of the inertia tensor due to the planet spin Ω, R is the planet volumetric radius, and k • is the planet secular Love number, then where G is the gravitational constant [see, for instance, Williams et al. (2001, Eq. 11)]. Therefore, the moment of inertia strain This explains the unusual dimension of sec −2 of the stiffness coefficient γ . Equations (9) and (10) We remark that some geometric quantities of the bodies, as the geometric radius R, will appear in this paper only because they are used in the definition of several astronomical figures as, for instance, the Love numbers. We stress that all these geometric quantities are meaningless within our theory. For a body of mass m and mean moment of inertia I • = (Tr I)/3, it is possible to construct a fictitious "radius," the mean moment of inertia radius, or simply the inertial radius [sometimes called "radius of gyration" (Chandrasekhar 1987)], defined by The ratio R I /R is the "inertial radius"/ "volumetric mean radius" of the body (0.4(R I /R) 2 is the usual quantity I • /m R 2 ). The inertial density of the Earth is ρ I E = 7332 kg/m 3 (the "geometric" mean density of the Earth is 5514 kg/m 3 ) This is the radius of a homogeneous ball of mass m and moment of inertia tensor I • I. If the body is incompressible and the deformations are small, then R I is constant in time. Using the inertial radius, an inertial density can be defined as As an example, consider a mass m of homogeneous inviscid liquid under self-gravity. At rest the liquid has a spherical shape and moment of inertia I • = 0.4 m R 2 , where R is the radius of equilibrium, which coincides with the inertial radius R I . For a homogeneous fluid body in hydrostatic equilibrium, it has been shown by Kelvin that k • = k f = 3/2 (Munk and MacDonald 1960, p. 26), where k f is the fluid Love number. From Eq. (11) The square of the angular frequency of free oscillations of the spherical mass of fluid, denoted as ω 2 f , is 4 5 Gm/R 3 [Lamb 1932, paragraph 262 Eq. (10)]; so γ f = ω 2 f . Within our theory, the free motion of a self-gravitating mass of perfect fluid is given by Eq. (7) with η = 0 and F = 0, μB = −γ B, which implies the frequency of free oscillations ω = √ γ /μ. So, for a homogeneous body made of a perfect fluid μ = 1. For any given body with inertial constants m and I • , one may consider a homogeneous mass of fluid of radius R = R I to define the values that can be taken as reference values for γ and μ, respectively. In Ragazzo and Ruiz (2015), the values of γ for the Sun and several planets were estimated using the spin rate and the dynamic form factor. The resulting values of γ divided by γ I are shown in Table 1. 2 This table shows that the approximation γ ≈ γ I is good for several celestial bodies, which can be explained by the fact that for large bodies self-gravity is dominant over rheological stress. If the approximation γ ≈ γ I is assumed and the spin rate of the body is constant and equal to Ω, then Eq.
where the direction of rotation was taken as e 3 . The small distortion hypothesis (b) implies that B 1, so a necessary condition for the validity of our theory is

Nondimensional parameters and the small inertia limit
In this section, we study the tide response equations (7) and (8) when μ is small, where F(t) is supposed to be time oscillating with a typical angular frequency ω. The goal is to compare the tide response when μ = 0 with that for μ > 0. Inertia, as small as it can be, is clearly important and dominant for "high-frequency" forcing. So, the first question to be answered in this section is: What does it mean "high-frequency" tide forcing?
In the previous section, we defined ω I as the natural frequency of oscillations of a mass m of homogeneous inviscid liquid under self-gravity. These oscillations arise from the exchange of gravitational energy and kinetic energy of deformation and, so, they are entirely due to the inertial properties of the body. The frequency ω I is also the frequency of resonance of the ideal body: If the forcing frequency ω is below ω I , then the motion is dominated by the self-gravitational force; if ω > ω I , then the motion is dominated by the inertial force; and if ω = ω I the balance of gravitational and inertial effects leads to unbounded oscillations. The scenario for a nonhomogeneous body is the same, in this case the natural frequency of oscillations is ω 0 = √ γ /μ. These considerations lead to the following definition: A tide forcing frequency ω is of "high-frequency" if ω > ω 0 = γ μ .
We expect that ω 0 be of the order of magnitude of ω I given in Eq. (15). If the angular velocity of the tide-raising body in the inertial frame is much smaller than the spin angular velocity Ω of the deformable body, then either ω ≈ jΩ (diurnal modes) or ω ≈ 2 jΩ (semi-diurnal modes), where j > 0 is an integer. So, due to the spin rate limit in Eq. (16), the fundamental frequency ( j = 1) of the tide forcing must always be of low frequency in order to the theory in this paper to be applicable. Notice that forcing frequencies ω = jΩ much larger than ω I , for j 1, are allowed within the theory as far as they do not induce large amplitude oscillations on B.
In order to evaluate the smallness of the deformation inertia, it is convenient to write Eqs. (7) and (8) in nondimensional form. The division of Eq. (7) by γ gives are characteristic times. As mentioned above, if the deformable body is made of an inviscid material, η = 0 in Fig. 5, then the equation for the free oscillations of B becomes μB+γ B = 0 that implies the natural period of free oscillations 2πτ 0 = 2π/ω 0 . The characteristic time τ 3 = η/γ is the relaxation time of the homogeneous equation ηḂ = −γ B obtained when μ = 0.
For the Maxwell rheology, the characteristic times are again obtained dividing Eq. (8) by γ : where τ 0 and τ 3 are defined as in Eq. (19) (in this case τ 3 is the relaxation time when μ = 0 and α = ∞) and is the relaxation time for the Maxwell element. In the small-μ regime, an important characteristic time is that of the homogeneous equation with μ = 0, that is In the next two subsections, we consider the small inertia limit for the Kelvin-Voigt and the Maxwell rheologies separately. Since the equation for each element of the matrix B is uncoupled from the others, it is possible to study the small inertia limit for each element of the matrix separately. So, it is enough to consider the equation for the Kelvin oscillator (2) in the case of the Kelvin-Voigt rheology and the equation for the Maxwell oscillator (3) in the case of the Maxwell rheology.

Kelvin-Voigt rheology
In order to study the singular limit as μ → 0 of the Kelvin oscillator (2), it is convenient to define a small nondimensional inertia parameter as Then, Eq. (2) for the Kelvin oscillator becomes τ 2 The small inertia regime can be studied using standard singular perturbation methods (Kaper 1999). For ε > 0 sufficiently small, it is possible to show the existence of an invariant manifold given by the graph of a function (w, is substituted into the first equation of (24) and the second equation is used to eliminateẇ then we obtain h(w, t) = w − f (t)/γ and x = w +ε(w − f (t)/γ )+O(ε 2 ). Substituting this expression for x into the first equation in (24), we obtain the differential equation τ 3 1+εẇ = −w + f (t) γ , which is similar to the equation obtained for ε = 0. So, for ε small, the solution B to the second-order equations (7) are approximately given by the solution to the following first-order equation:

Maxwell rheology
The analysis of the singular limit μ → 0 of the Maxwell oscillator is different from that of the Kelvin oscillator. While the characteristic equation associated with the Kelvin oscillator in Eq.
(2) has a real eigenvalue proportional to −1/μ, which implies fast attraction to the invariant manifold (w, t) → x = w + εh(w, t) + O(ε 2 ) of the previous paragraph, the characteristic equation associated with the Maxwell oscillator in Eq.
(3) has a pair of complex eigenvalues proportional to ±i/ √ μ. Therefore, as μ → 0 most solutions to Eq. (3) have a fast oscillatory component and the analysis of the singular limit, to be made in Appendix A, is by the method of averaging (Hale 1980;Sanders et al. 2007). The result of the averaging procedure is presented in the following. Let ε be the nondimensional small inertia parameter defined as 3 After three steps of averaging, see Appendix A, we obtain that for ε small, and after a transient where the fast oscillations of B decay as exp[−tτ 3 /(2τ 1 τ 2 )], the solutions B to Eq. (8) are approximately given by: We remark that the first and second derivatives of F must be of the order of magnitude of F in order to guarantee the accuracy of the approximation.

Love numbers, the Correspondence Principle and the Association Principle
In this section, we study the tides on a deformable body induced by a single harmonic force due to a moving point mass. All the results in this section hold for quite general rheologies, which include the Maxwell and the Kelvin-Voigt rheologies previously discussed. So, we start with a brief presentation of the concept of complex compliance in linear viscoelasticity. The one-dimensional force-extension relation for a spring-dashpot system, which models the rheology of a linear viscoelastic material, is a linear differential equation relating force f and extension x. For instance, the force-extension relation of the Maxwell rheology represented in Fig. 2 where α represents the spring rigidity and η the damping coefficient of the dashpot. Imposing a harmonic force f (t) =f (ω)e iωt and a harmonic deformation x(t) =x(ω)e iωt , the force-extension relation implies the equation where the function J (ω) is called the complex compliance (1/J (ω) is the complex rigidity). For instance, for the Maxwell rheology the complex compliance is J (ω) = α −1 + (iωη) −1 .
For the generalized Voigt model shown in Fig. 3 the complex compliance is [Bland 1960, chapter 1 Eq. (27)] and for the generalized Maxwell model shown in Fig. 4 the complex rigidity is [Bland 1960, As discussed in Sect. 2, the effect of self-gravitation can be added to the spring-dashpot system by placing a spring with elastic constant γ in parallel with the part of the system that represents the rheology. Suppose that a harmonic forcing f =f (ω)e iωt acts upon both the spring γ and the subsystem that represents the rheology generating a harmonic displacement x =x(ω)e iωt . Iff γ (ω)e iωt is the part of f that acts upon the spring, then γx(ω) =f γ (ω). Iff r (ω)e iωt is the part of f that acts upon the subsystem that represents the rheology, then is the complex rigidity of the system that represents the rheology. Therefore,x(ω) =f (ω)/[γ + J −1 (ω)] and the complex compliance of the whole system becomes 1 γ + J −1 (ω) .
As argued in Sect. 2, the effect of deformation inertia can be added to the spring-dashpot system of the previous paragraph by placing a mass element μ in series to the system. The construction gives rise to an oscillator as shown in Fig. 7. The dynamics of the deformation matrix B can be obtained from the dynamics of the oscillator by means of the Association Principle in (4): "The differential equation for B is the differential equation for the springdashpot oscillator after replacing x by B." If the deformable body spins steadily, Ω is constant, then the tidal force that drives the motion of B is given by the last term in Eq. (5), namely Each element of the matrix F can be decomposed into Fourier modes of the formf e iωt . Letbe iωt be the component of the matrix B that corresponds tof e iωt by means of the tide response equations. Due to the Association Principle, the relation betweenb andf is the same as that betweenx andf determined by the oscillator in Fig. 7 As in Eq. (11), the complex Love number associated with the frequency ω is given by the following nondimensionalization of the ratiof /b [see Ragazzo and Ruiz 2017, equation (46)] where R is the geometric (volumetric) radius of the body. If the secular Love number is defined as the value of k(ω) at ω = 0, The spring with elastic constant γ represents self-gravity, and the mass μ represents deformation inertia where Notice that for the Maxwell and the Kelvin-Voigt rheologies τ 0 is the characteristic time defined in Eq. (19). According to Efroimsky (2012, section 4.3 and Appendix B), the Correspondence Principle of linear viscoelasticity was pioneered by G. Darwin and, within the context of Love numbers, can be explained as follows (see Efroimsky 2012 for details). The one-dimensional stress-strain constitutive relation of a linear viscoelastic material is a linear differential equation relating the stress σ and the strain ε. For instance, the constitutive relation for the Maxwell rheology isε(t) = α −1σ (t) + η −1 σ (t), where α is a modulus of elasticity and η a coefficient of viscosity. Imposing a harmonic stress σ (t) =σ (ω)e iωt and a harmonic strain ε(t) =ε(ω)e iωt , the constitutive relation implies a complex relation of the form where the functionJ (ω) is a complex compliance (1/J (ω) is a complex rigidity) and ε i j and σ i j are the deviatoric components of the strain and stress tensors, respectively. The static Love number of degree 2 of a homogeneous, incompressible, and elastic spherical body of radius R, density ρ, surface gravity g, and elastic rigidityμ, is given by Within the context of Love numbers, the Correspondence Principle of linear viscoelasticity amounts to change the static rigidityμ in Eq. (36) by the complex rigidityJ −1 (ω), so that the "dynamic" complex Love number k 2 (ω) becomes In order to compare Eqs. (31) and (37), it is necessary to make μ = 0 in Eq. (31), since deformation inertia is neglected in Eq. (37). Moreover, it is necessary to suppose that γ in Eq. (31) is that of a homogeneous mass m of inviscid liquid, see Eq. (14), where k f = 3/2 is the fluid Love number. The substitutions γ = γ f and μ = 0 into Eq. (31) give Equations (37) and (38) give the same Love number after the identification Notice thatJ −1 has the same dimension kg/(m s 2 ) as the shear stress, while J −1 has the same dimension 1/s 2 as the moment of inertia stress [see text between Eqs. (10) and (11) for an explanation]. Equation (39) allows for the comparison of the rigidity and damping coefficients used in the Association Principle formalism and the Correspondence Principle formalism.
The main difference of our 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. In the case of the Association Principle, the scalar Love numbers are easily obtained and are the same as those obtained from the Correspondence Principle. In the case of the Correspondence Principle, except for simple rheologies, it is not obvious how to obtain from the scalar Love numbers the three-dimensional time-domain equations for the coupled spinorbit motion of many bodies. As far as we know, the only papers to fully accomplish this task, exclusively for the Maxwell and the Kelvin-Voigt rheologies, are Correia et al. (2014) and Boué et al. (2016), for the two-and three-dimensional cases, respectively.

Love numbers and the effect of inertia
The main goal of this paper is to compare the tide response when μ = 0 with that when μ > 0. Equation (33) implies that k(ω, μ = 0)/k • = ζ(ω); therefore, this comparison amounts to evaluate the relative difference The following proposition is crucial.
Proposition 1 For any rheology associated with a generalized Voigt model, as given in Eq.
The same estimate in Eq. (42) holds for rheologies with a continuous distribution of characteristic times, as for instance the Andrade rheology, for the proof being sufficient to take the limit as n → ∞ in either Eq. (29) or (30) in the appropriate way (see Bland 1960, chapter 1.2]. Using the definition of high-frequency tidal forcing in Eq. (17), we conclude, as expected, that inertia is negligible for low-frequency tidal forcing. The Love number for the Maxwell rheology is obtained from the generalized Voigt rheology making n = 0 in Eq. (29), where Using that τ 2 = τ 1 + τ 3 , with τ 1 = η/α and τ 3 = η/γ , and that the Kelvin-Voigt rheology is equal to the Maxwell rheology when τ 1 = 0, we obtain the Love number for the Kelvin-Voigt rheology: In order to present graphs of Love numbers, it is convenient to define a nondimensional frequency ωτ 2 . In this case, the expression for the Love numbers becomes: Figures 8 and 9 show the Love numbers for the Kelvin-Voigt rheology, τ 1 = 0. In this case, ε = τ 2 0 (τ 1 +τ 3 ) 2 = τ 2 0 τ 2 3 . Figures 10 and 11 show the Love numbers for the Maxwell rheology and for the ratio τ 1 /τ 2 = 1/3 also used in Correia et al. (2014). Both the absolute value of k/k •  Fig. 8 Absolute value of the Love number |k/k • | for τ 1 = 0 (Kelvin-Voigt rheology), as a function of the nondimensional angular frequency ωτ 2 , and for several values of ε = τ 2 0 /τ 2 2 . The vertical lines indicate the value of twice the initial spin angular speed 2Ω in (initial semi-diurnal forcing frequency) for the two values τ 2 = 10 −2 and τ 2 = 10 −4 (years) used in the numerical simulations in Sect. 6  (47) Fig. 11 Quality factor Q, as defined in Eq. (47), for the Maxwell rheology (τ 1 /τ 2 = 1/3) and the same values of ε s as in Fig. 10. The vertical lines have the same meaning as in Fig. 8

Application to HD 8060b
The results in Sect. 5 give an answer to the question on "how relevant is deformation inertia" in the case the tide forcing is monochromatic and constant in time. The analysis of the same question for the full set of equations, which couples tide and orbital elements, is much harder.
There are two types of difficulties. The first is due to the harmonic decomposition of the tide forcing, unless the orbit of the tide-raising body is circular, there are always harmonics of arbitrarily large frequency. Although inertia is definitely important at high frequencies, this importance may be neutralized by the smallness of the amplitude of the high-frequency harmonics. The second issue is that frequencies and amplitudes vary with time. They are damped by tidal effects and nonlinearly coupled by the equations for the orbit evolution. As a consequence, the dominance of different harmonic modes may change with time. It is difficult to analytically study the dynamics of the full set of equations. For this reason, in this section, a numerical analysis of Eqs. (6), (5), and either (7) or (8) is carried on. The deformable body that will be used as a model in the simulations is the HD 80606 b, adopting the same parameters presented in the integrations of Correia et al. (2014) and Boué et al. (2016), except by the new parameter μ. These values were already used to make the graphs in Figs. 8, 9, 10, and 11. The value of μ will be sometimes exaggerated with respect to the reference value μ = 1. In some situations, this may be compensated by other unrealistic values of τ 2 , but this will not be explicitly analyzed, since the main goal in this section is to test the validity of the criterion in (42), for neglecting inertia, under the dynamics of the fully coupled equations of motion. Following Correia et al. (2014), we use the orbital parameters of the HD 80606 b to give the initial conditions of the system. The initial semi-major axis and eccentricity are a = 0.455 au and e = 0.9330, respectively. The initial rotation period is 0.5 day. We choose a fixed initial obliquity θ = 72 • . Using these initial conditions, we get the graphs in Figs. 12 and 13 (τ 3 = 10 −4 year). The physical parameters of the HD 80606 b (Correia et al. 2014) are: mass m 2 = 1297. m ⊕ = 7.746·10 28 kg, radius R = 10.75R ⊕ = 68488.3 km, where the symbol ⊕ denotes that the parameters are relative to the Earth. The moment of inertia of the planet is given by I • = ξ m 1 R 2 = 8.1527 · 10 40 kg m 2 , where the parameter ξ is computed through the Darwin-Radau equation . For the value μ = 4, we present both the exact and the approximated solutions verifying that they coincide in this scale. We remark that the obliquity of the orbit stabilizes in θ = 0 (deg) in a very short interval of time relatively to the other quantities, even for μ = 6000. Moreover, the dissipation rate does not depend monotonically on μ. Notice that even for μ = 4 the solution is still close to that for μ = 0 Fig. 13 The graph of the eccentricity as function of time for τ 3 = 10 −4 (year). The continuum lines represent the Kelvin-Voigt rheology (τ 1 = 0), and the dashed lines represent the Maxwell rheology (τ 1 /τ 2 = 1/3). We can observe that the power initially dissipated by the Kelvin-Voigt rheology is larger than that for the Maxwell rheology; the same can be observed in  where we choose the value k • = 0.5, which is the fluid Love number of Jupiter [see page 9 of Correia et al. (2014)]. So, using Eq. (11), we compute the value of γ = 2.48041 · 10 9 (1/year 2 ). The mass of the central star for this planet is m 1 = 2008.9 · 10 30 kg The initial conditions adopted in Fig. 14 (τ 3 = 10 −2 year) were obtained integrating the system with μ = 0 and the initial conditions described above until the obliquity goes to zero (planar motion), and the eccentricity reaches e ≈ 0.7. We used this initial condition to make this integration easy and to focus on the spin-orbit resonances. The ratio τ 3 /τ 2 = 2/3 was used in all simulations for the Maxwell rheology; this value was taken from Correia et al. (2014, see Fig. 2(a)) and it is based on the values for the Earth.
Some comments on the choices of μ are necessary at this point. For τ 3 = 10 −2 and τ 1 = 0 (year) (Kelvin-Voigt rheology), it is difficult to distinguish the solutions for μ ≈ 1 from the solutions for μ = 0, so we chose an exaggerated value μ = 62010 that corresponds to the critical value for the underdamped/overdamped transition of the harmonic oscillator Eq. (2). The solutions for μ = 62010 deviate considerably from those with μ = 0. Later we verified that even 10 % of this value, μ = 600, changes visibly the solutions for the Maxwell rheology but not for Kelvin-Voigt rheology. This fact can be seen by the direct integration presented in Fig. 14 (bottom-right), where we compare the results for μ = 600 and μ = 62010 and Kelvin-Voigt rheology with the result for μ = 600 and Maxwell rheology, since we can see that the graph for the Maxwell rheology is closer to the graph μ = 62010 than to the corresponding μ = 600 in Kelvin rheology. This particular case is also discussed in the next section, Fig. 16. However, for τ 3 = 10 −4 year, we can observe in Fig. 13 that the curves of Kelvin-Voigt (continuous lines) and Maxwell (dashed lines) are very close for each fixed μ, but as we increase these values of μ the solutions for both rheologies present considerable variations, until the system begins to present captures into spin-orbit resonances for μ = 2000.
A comparison between the eccentricities obtained integrating the exact equations and the small inertia approximations given in Eq. (25) (Kelvin-Voigt rheology) and Eq. (27) (Maxwell rheology) is shown in Fig. 15.

The secular torque
In this section, using essentially the same approach as in Correia et al. (2014), we study how the (planar) secular torque is modified by the inertia coefficient for the HD 80606 b.
The equation of motion for the deformation variable, Eq. (20), can be written as a single third-order equation Our objective is to consider the planar motion of the system, Fourier expand the term F, solve Eq. (49) for B, substitute the result into the equation of the angular velocity [first of Eq. (6)] and take the average of this last expression. The resulting term of the right side of this equation is called the secular torque of the system. Considering the expression (5) for F, we assume that the variable q has a much faster variation than Ω, in a fixed orbit of the Kepler problem. Denoting q = r (cos ϕ, sin ϕ, 0), where ϕ is the polar angle in the body frame, and assuming that the spin is constant and aligned with the z-axis, we get the formula where C is a constant, diagonal matrix. This term generates a constant and diagonal stationary solution that commutes with q ⊗q and so do not contribute to the secular torque; then we will skip it henceforward. Compared with the notation of Correia et al. (2014), we see that the Fig. 14 The graphs of the semi-major axis (top-left), the eccentricity (top-right) and the spin-orbit ratio, Ω/n (bottom-left) as a function of time (years), for the Kelvin-Voigt rheology (τ 1 = 0) with τ 3 = 10 −2 (year) and for ε = 0 (μ = 0), ε = 4 × 10 −6 (μ = 1), ε = 2.4 × 10 −3 (μ = 600), and ε = 0.25 (μ = 62010). For the values μ = 1 and μ = 600, only the approximated equations (25) were integrated, since the integration of the exact equations is prohibitively slow. Notice that inertia is relevant only for the exaggerated value μ = 62010. The graph of the spin-orbit ratio shows that for μ = 0 the system is captured into spin-orbit resonances, as observed in Boué et al. (2016). Such captures do not disappear with the addition of the inertia term, even for large values of μ, as μ = 62010. In the last graph, (bottom-right), we see a large difference between the Maxwell and Kelvin-Voigt rheologies in the evolution of the eccentricity for μ = 600; this difference is suggested by the corresponding secular torques in Fig. 16 Kelvin Maxwell Fig. 15 The graphs of the exact and approximated solutions for the Kelvin-Voigt (left), τ 1 = 0, and Maxwell (right), τ 1 /τ 2 = 1/3, rheologies, for τ 3 = 10 −4 (year). Here, we integrate all the systems beginning with e = 0.2 and no obliquity, including μ = 0, and so we plot the value e = e(μ) − e(0) as function of time, where e(μ) is the eccentricity of the orbit for the given value of μ. We remark that the approximation is good even for μ = 4. The μ values corresponding to ε = 0.01, 0.02, 0.04 are μ = 0.25, 0.5, 1, respectively (Eq. 23), and the μ values corresponding to ε = 0.02, 0.07 are μ = 1, 4, respectively (Eq. 26) coordinates of the matrix (50) are related to formulas (21) and (22) where M is the mean anomaly,θ = θ − , with θ denoting the rotation angle of the planet in an inertial frame and the longitude of the pericenter, see Fig. 1 of Correia et al. (2014). The coefficients β k are given by where X −3,2 k (e) are Hansen coefficients, polynomials in the eccentricity of the orbit as described in Appendix C of Correia et al. (2014), and a is the semi-major axis of the orbit.
Thus, the nonconstant stationary solution of Eq. (49) can be written as where ω k = 2θ − k M. Using these simplifications and omitting the negligible terms of the first equation of (6), we getΩ and so we can reduce this equation to a single scalar ODË Equation (57) is only a generalization of Eq. (39) of Correia et al. (2014). The secular torque is obtained by averaging Eq. (57) over M, i.e., by considering only the term j = 0. Hence, the secular torque T is given by where year 3 Fig. 16 The graphs of the normalized torque, Eq. (58), as a function of the spin-orbit ratio Ω/n, for τ 3 = 10 −2 (year). We present the graphs of torque for three different values of eccentricity e = 0.2, e = 0.4, e = 0.6 and the same values of μ as in Fig. 14. We compare the differences between the Kelvin-Voigt rheology (left) and the Maxwell rheology (right). For the Kelvin-Voigt rheology, we cannot see any difference between the curves of μ = 0, 1, 600, for all eccentricities, as verified in the three first graphs of Fig. 14. However, for the Maxwell rheology we see a large difference for the values μ ≥ 600, as verified in the last graph (bottom-right) of Fig. 14 for the case μ = 600. We recall that each point where the graph crosses the horizontal axis, with negative derivative, corresponds to a spin-orbit capture, which is compatible with the captures plotted in the (bottom-left) frame of Fig. 14; however, we remind that such captures observed in the integrations take place at different values of e, not represented in only one of the graphs above In Figs. 16 and 17, we present the graphs of the normalized torque, adopting for simplicity K = 1, as a function of the spin-orbit ratio Ω/n, for τ 3 = 10 −2 (year) and τ 3 = 10 −4 (year), respectively. In order to express the torque (58) as a function of this ratio, we need to use explicitly the value of n(e), the orbital angular velocity (derivative of the mean anomaly M). Since this angular velocity is not constant along the motion, we assume that the orbital angular momentum is approximately conserved and hence n becomes a function only on the eccentricity year 3 Fig. 17 This figure is similar to Fig. 16 but for τ 3 = 10 −4 (year). The values of eccentricity are the same as in Fig. 16, but the μ s are different: μ = 20, 80, 500, 2000, 6000. The behavior for μ = 0, 1, 4 is the same as that for μ = 20, i.e., the curve intersects the horizontal axis at only one point (Ω/n > 1), corresponding to the pseudo-synchronous equilibrium. Such curves, for μ < 20, are not presented in this figure because they are not visible in the present scale. Both the Kelvin-Voigt and Maxwell rheologies have roughly the same aspects, except by the intensity of the torque, which is bigger in the Maxwell case. The appearance of roots for μ ≥ 500 is the main point of this set of graphs, in agreement with the graph (bottom-left) of Fig. 12, where we see that for μ = 2000, 6000 the solutions exhibit resonances that are not present for μ ≤ 500 n = n(e) = (Gm 1 m 2 ) 2 where L = |x ×ẋ| here denotes the normalized initial orbital angular momentum used in the integration.

Conclusion
The main result in this paper is that for a very broad class of rheologies; namely, those associated with a generalized Voigt model Eq. (29), or to a generalized Maxwell model Eq. (30), or even to a rheology with a continuous distribution of characteristic times, as for instance the Andrade rheology, obtained from continuous limits of either a generalized Voigt or Maxwell model; deformation inertia is negligible if where ω 0 = √ γ /μ is the natural frequency of oscillation of the system when damping is neglected and ω is the angular frequency of a harmonic tidal force. Moreover, a quantitative statement of this result holds: If k(ω, μ) is the complex Love number associated with a given rheology, then ω 2 ω 2 0 < 0.1 p ⇒ k(ω, μ) − k(ω, 0) k(ω, 0) < 0.11 p for any p ∈ [0, 1].
As argued in Sect. 2.1, given the mass m and the moment of inertia I • of the deformable body, the value of ω 0 shall be of the order of magnitude of the reference value given in Eq. (15), As mentioned in Sect. 2.1, if the angular velocity of the tide-raising body in the inertial frame is much smaller than the spin angular velocity Ω of the deformable body, then the most relevant tide-forcing frequency is the semi-diurnal ω ≈ 2Ω. Since Ω 2 must be much smaller than γ ≈ γ I = ω I , Eq. (16), in order to guarantee a small flatness of the body due to the spin, we conclude that for the main semi-diurnal frequency inertia is always negligible. Moreover, our simulations in Sect. 6 indicate that higher harmonic tide forcing that eventually appear due to the nonlinearity of the full equations of motion do not trigger any significant effect of deformation inertia. So, from the point of view of orbital evolution of celestial bodies, inertia can be safely neglected. (Exceptions may occur when a higher-order harmonic of the tide forcing have a high amplitude.) The above considerations may lead to the deceptive conclusion that deformation inertia is irrelevant to the whole issue of tides in celestial mechanics. In the following, we argue that it may be important in some aspects of the theory.
Inertia introduces high-frequency oscillations near the resonance frequency 4 ω 0 . Energy dissipation shifts the resonance frequency and sets its amplitude and width. In principle, the oscillations at resonance frequency may be forced by small amplitude higher harmonics of the tide forcing. If it were possible to estimate accurately, the amplitude of the tide forcing at this frequency, since γ and μ can be well estimated from the density radial profile, it would be possible to use Eq. (31) to obtain information about the complex compliance J of the deformable body at the resonance frequency.
It has been observed that normal modes of the Earth (within the 2 to 20 mHz frequency band) are continuously excited, but the sources of this background oscillations, the so-called Earth hum, have remained unresolved (Nishida 2013). The measured frequency of free oscillations of the Earth, excited by large earthquakes and associated with the spheroidal mode S 0 2 , is 54 min (0.309 Hz). This frequency is the one that corresponds to the resonance frequency (including the effect of damping) in our model. Unfortunately, the 0.309 Hz frequency is outside the frequency range measured for the Earth's hum. If measurements of the Earth background oscillations are improved to include the 0.309 Hz frequency, then we may be able to improve our knowledge of the rheology of the Earth using Eq. (31) and the estimates of μ and γ . Moreover, the recorded transient motion of the 54-min oscillations after large earthquakes can be used in conjunction with the time-domain equations of motion derived in Ragazzo and Ruiz (2017) to better understand the rheology of the Earth.