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 (2016) and Ragazzo and Ruiz (2017). The singular limit as the inertia tends to zero is analysed 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)) (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 to this Love number, which is proportional to the deformation inertia, is not only relevant but dominant over the 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. In Love (1911) chapter 5, Love 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. 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 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.
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 to 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 time domain 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 to tide deformations in the later. 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 Sec. 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 Sec. 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 Sec. 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 Sec. 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. Sec. 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 Sec. 4. Sec. 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 are computed and compared. Sec. 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.

Equations of motion
In this section we present the equations of motion, as derived in Ragazzo and Ruiz (2015) and Ragazzo and Ruiz (2017). These equations are valid under the following hypotheses: (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 to the motion of K is given by 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 experience 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 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 self-gravitation, there must be an inertia associated to 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 to 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 to 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 to 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 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 Sec 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 steady 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 γ. Eqs. (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 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 selfgravity. At rest the liquid has a spherical shape and moment of inertia I • = 0.4 mR 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) Table 1 Values of γ/γ I from Ragazzo and Ruiz (2015) Table 1. 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 • /mR 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 ). 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. (9) 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 3 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,t) → x = w + εh(w,t) + O(ε 2 ), which is ε-close to the constraint x = w valid for ε = 0. If the expression x = w + εh(w,t) + O(ε 2 ) 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 to 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 to 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 Eqs. (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.
4 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 3 Notice that the ε defined in Eq. (23) can also be written as ε = τ 2 0 /τ 2 2 , since for the Kelvin-Voigt rheology τ 1 = 0 and τ 2 = τ 1 + τ 3 = τ 3 . So, there is no ambiguity in using the same letter in both Eqs. (23) and (26).
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 showed in Fig. 3 the complex compliance is (Bland (1960) chapter 1 eq. (27)) and for the generalized Maxwell model showed in Fig. 4 the complex rigidity is (Bland (1960) chapter 1 eq. (29)) As discussed in Sect. 2, the effect of self-gravitation can be added to the springdashpot system by placing a spring with elastic constant γ in parallel to 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 J −1 (ω)x(ω) =f r (ω), where J −1 (ω) 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 Section 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 spring-dashpot 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 The spring with elastic constant γ represents self-gravity and the mass µ represents deformation inertia.
b andf is the same as that betweenx andf determined by the oscillator in Fig 7, As in Eq. (11), the complex Love number associated to 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, 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 Eqs. (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 equations (10) and (11) for an explanation). Eq. (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 spin-orbit 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. Eq. (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 to a generalized Voigt model, as given in Eq. (29), or to a generalized Maxwell model, as given in Eq. (30), the following inequality holds: |ζ (ω)| ≤ 1, for all ω ∈ R.
Proof: Notice that for the generalized Voigt model J −1 (0) = 0 and for the general- where a and b are real valued, then from Eq. (34) .
and again the positivity of the α ′ s imply a(ω) ≥ 0. The proposition is proved. Proposition 1 implies that the difference in Eq. (40) is small if τ 2 0 ω 2 is small. Using that τ 0 = 1/ω 0 we obtain So, Deformation inertia can be neglected if ω ≪ ω 0 .
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), 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 become: 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 0 0  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 Sec. 6. 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 non linearly 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 Fig. 8, Fig. 9, Fig. 10, and Fig. 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 analysed, 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 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 the 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 were 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. Fig. 12 The graphs of the semi-major axis (top -left), the eccentricity (top -right), the spin-orbit ratio (bottom -left) and the obliquity (bottom -right) as a function of time (years), for the Kelvin-Voigt rheology, τ 1 = 0, and τ 2 = 10 −4 (year). The values ε = 0 (µ = 0), ε = 0.16 (µ = 4), ε = 0.81 (µ = 20), ε = 3.23 (µ = 80), ε = 20.25 (µ = 500), ε = 81 (µ = 2000), and ε = 243 (µ = 6000), correspond to those used in Figs. 8 and 9. 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 Fig. 14 (bottom-right).
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 equation (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 (bottomright), 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 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, (bottomright) 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. secular torque, then we will skip it henceforward. Comparing with the notation of Correia et al (2014), we see that the coordinates of the matrix (50) are related to the 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 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)).
where X −3,2 k (e) are Hansen coefficients, polynomials in the eccentricity of the orbit as described in the appendix C of Correia et al (2014), and a is the semi-major axis of the orbit.
Thus, the non constant stationary solution of Eq.(49) can be written as where ω k = 2θ − kM.
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 equation (57) over M, i.e., by considering only the term j = 0. Hence, the secular torque T is given by where In the 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 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 to a generalized Voigt model Eq. (29) 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 to a given rheology, then ω 2 ω 2 0 < 0.1p =⇒ k(ω, µ) − k(ω, 0) k(ω, 0) < 0.11p for any p ∈ [0, 1].
As argued in Sec. 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 Sec. 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, 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 can not 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, what 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.
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 Sec. 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 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. 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 to the spheroidal mode S 0 2 , is 54 minutes (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 minutes 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.