© ESO, 2014

1. Introduction

The occurrence, on most open ocean coasts, of high sea tide at about the time of the Moon’s passage across the meridian, gave ancient observers the idea that our satellite exerts an attraction on the water, although the occurrence of a second high tide when the Moon is on the opposite meridian was a great puzzle. The correct explanation for tides was first indicated by Newton (1760) in the Philosophiæ Naturalis Principia Mathematica. They result from a differential effect of the gravitational force acting in accordance with laws of mechanics, since this force is not constant across the diameter of the Earth. Kant (1754) suggested that the Moon not only pulls the Earth, but also exerts a retarding torque upon its surface, this torque slowing down the Earth’s rotation until terrestrial days become as long as lunar months.

The estimates for the tidal deformation of a body are based on a very general formulation of the tidal potential, initiated by Darwin (1879). Assuming a homogeneous Earth consisting of an incompressible fluid with constant viscosity, Darwin derived a tide-generated disturbing potential expanded into a Fourier series. Thus, he was able to obtain the tidal contribution to the spin and orbital evolution, confirming Kant’s prediction (Darwin 1880). Ever since, many authors improved Darwin’s work, leading to a better understanding of this problem (e.g., Kaula 1964; MacDonald 1964; Alexander 1973; Singer 1968; Mignard 1979; Efroimsky & Williams 2009; Ferraz-Mello 2013).

All previous studies have shown the existence of a stationary solution for the rotation rate, which is synchronous when the two bodies move in circular orbits, but becomes super-synchronous for elliptical orbits. Standard tidal models consider an elastic tide and delay the tidal bulge by an ad hoc phase lag in order to take into account the body anelasticity (e.g., Kaula 1964; MacDonald 1964). However, these models are unable to explain the dependency of phase lag with the tidal frequency, which is a crucial point when we aim to study the long-term secular evolution of the planets.

If one adopts a linear model where the phase lag is proportional to the tidal frequency, i.e, a viscous model, the ratio between the orbital and rotational periods have an excess proportional to the square of the eccentricity, an equilibrium often called pseudo-synchronous rotation (e.g., Mignard 1979). However, this prediction is not confirmed by observations of the rotation of most planetary satellites (which are synchronous in non-eccentric orbits) or of the planet Mercury. One way to obtain a synchronous stationary solution is to assume that an additional torque, due to a permanent equatorial asymmetry (a quadrupolar approximation for the body’s figure), is acting on the body to counterbalance the tidal torque (Colombo 1965; Goldreich & Peale 1966). The stationary solutions obtained from the combination of both torques are called spin-orbit resonances, for which the synchronous motion is a particular case.

A more realistic approach to deal with the dependency of the phase lag with the tidal frequency is to assume a viscoelastic rheological model. One of the simplest models of this kind is to consider that the planet behaves like a Maxwell material, i.e., the material is represented by a purely viscous damper and a purely elastic spring connected in series (e.g., Turcotte & Schubert 2002). In this case, the planet can respond as an elastic solid or as a viscous fluid, depending on the frequency of the perturbation. Viscoelastic rheologies have been used recently (e.g., Henning et al. 2009; Remus et al. 2012a), since they are able to reproduce the main features of tidal dissipation, including the pseudo-synchronous rotation or the spin-orbit equilibria (e.g. Storch & Lai 2014). For a review of the main viscoelastic models see Henning et al. (2009).

Efroimsky (2012) describes the rheology of a planet using an empirical power scaling law for the phase lag which is consistent with the accumulated geophysical, seismological, and geodetic observational data. It is shown that the best fitted laws are in conformity with Andrade’s model (Andrade 1910), where the phase lag is inversely proportional to the tidal frequency. The Andrade model can be thought of as the Maxwell model equipped with an extra term describing hereditary reaction of strain to stress (Efroimsky 2012). However, Andrade’s model is extremely complex and the main features of tidal evolution are similar to other viscoelastic models.

Ferraz-Mello (2013) abandons the potential description introduced by Darwin, and proposes to study tides directly from the deformation law of the planet. For that purpose, he adopts a simple Newtonian creep model for the surface deformation, where the distance to the equilibrium is considered as proportional to the stress. As a result, the planet’s surface tends to the equilibrium spheroid, but not instantaneously. The resulting stationary rotations have an average excess velocity which depends on a critical frequency, the relaxation factor, which is inversely proportional to the viscosity of the body. Moreover, the dissipation in the pseudo-synchronous solution presents two regimes of variation, depending on the tidal frequency: in the inviscid limit, it is roughly proportional to the frequency (as in standard theories), but that behavior is inverted when the viscosity is high and the tidal frequency larger than the critical frequency, as in the model proposed by Efroimsky (2012).

An important setback in Ferraz-Mello (2013) model is the absence of an elastic tide. Indeed, when the tidal frequency is larger than the critical frequency, the phase lag tends to 90°, which is in contradiction with the observations done for most rocky planets. In order to conciliate the theory and the observed tidal bulges, we have to assume that the tide is not restricted to the component due to the creeping of the body under the tidal action, but has also a pure elastic component. This supposition was already present in Darwin’s original work (Darwin 1879), who realized that viscoelastic deformations of the planet’s surface could be derived from the Newtonian creep law.

In this paper we pursue in the same direction of Ferraz-Mello (2013), but we adopt a deformation law for the planet that is compatible with the viscoeslatic rheology of a Maxwell body (Sect. 2). We then compute the general solution for the deformation of the planet (Sect. 3), and use it to obtain the tidal evolution of the planet (Sect. 4). We also identify the equilibrium positions for the tidal torque (Sect. 5), and inspect the effect from some libration of the rotation rate around these equilibria (Sect. 6). We finally perform numerical simulations for some known close-in exoplanets (Sect. 7) and discuss the results in the last section.

2. Model

Fig. 1 Reference angles and notations. (i,j,k) is a fixed inertial reference frame, while (I,J,K) is a reference frame fixed in the planet. k = K are orthogonal to the orbit, and thus not shown. r is the radius vector from m to m0, P is the direction of pericenter, ϖ is the longitude of pericenter, and v is the true anomaly. The rotation angle θ is referred to the fixed reference direction i. We also use the angles $\begin{array}{c}{}_{\mathrm{˜}}\\ \mathit{\theta }\end{array}\mathrm{=}\mathit{\theta }\mathrm{-}\mathit{\varpi }$, f = v + ϖ, and γ = θ − f.

We consider a system consisting of a central star of mass m₀, and a companion planet of mass m. The planet is considered an oblate ellipsoid with gravity field coefficients given by J₂, C₂₂ and S₂₂, sustained by the reference frame (I,J,K), where K is the axis of maximal inertia (Fig. 1). We furthermore assume that the spin axis of the planet, with rotation rate Ω, is also along K (gyroscopic approximation), and that K is orthogonal to the orbital plane (which corresponds to zero obliquity). The gravitational potential of the planet at the star is then given by (e.g., Correia & Rodríguez 2013) $\begin{array}{ccc}\mathit{V}\mathrm{\left(}{r}\mathrm{\right)}& \mathrm{=}& \end{array}$(1)where $\begin{array}{ccc}\cos \mathrm{2}\mathit{\gamma }\mathrm{=}\mathrm{(}{I}\mathrm{·}{r̂}{\mathrm{)}}^{\mathrm{2}}\mathrm{-}\mathrm{(}{J}\mathrm{·}{r̂}{\mathrm{)}}^{\mathrm{2}} \mathrm{and} \sin \mathrm{2}\mathit{\gamma }\mathrm{=}\mathrm{-}\mathrm{2}\mathrm{(}{I}\mathrm{·}{r̂}\mathrm{)}\mathrm{(}{J}\mathrm{·}{r̂}\mathrm{)}\mathit{,}& & \end{array}$(2)G is the gravitational constant, R is the mean planet radius, r is the position of the star with respect to center of mass of the planet, and is the unit vector. We neglect terms in (R/r)³ (quadrupolar approximation).

2.1. Deformation of the planet

The mass distribution inside the planet, characterized by the gravity field coefficients, is a result of the forces acting on it, that is, self gravity and the planet’s response to any perturbing potential, V_p. Here, we consider that the planet deforms under the action of the centrifugal and tidal potentials. Thus, on the planet’s surface, R, the non-spherical contribution of the perturbing potential is given by (e.g., Correia & Rodríguez 2013) $\begin{array}{ccc}{\mathit{V}}_{\mathrm{p}}\mathrm{\left(}{R}\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}{\mathrm{\Omega }}^{\mathrm{2}}\mathrm{(}{K}\mathrm{·}{R}{\mathrm{)}}^{\mathrm{2}}\mathrm{-}\frac{\mathrm{3}\mathit{G}{\mathit{m}}_{\mathrm{0}}}{\mathrm{2}{\mathit{r}}^{\mathrm{3}}}\mathrm{(}{r̂}\mathrm{·}{R}{\mathrm{)}}^{\mathrm{2}}\mathit{.}& & \end{array}$(3)We let ζ(R) be the radial deformation of the free surface. In the case of a homogeneous incompressible elastic solid sphere, the equilibrium surface deformation ζ_e(R) is given by (Kelvin 1863a,b; Love 1911) $\begin{array}{ccc}{\mathit{\zeta }}_{\mathrm{e}}\mathrm{\left(}{R}\mathrm{\right)}\mathrm{=}\mathrm{-}{\mathit{h}}_{\mathrm{2}}{\mathit{V}}_{\mathrm{p}}\mathrm{\left(}{R}\mathrm{\right)}\mathit{/}\mathit{g,}& & \end{array}$(4)where g = Gm/R² is the surface gravity, and $\begin{array}{ccc}{\mathit{h}}_{\mathrm{2}}\mathrm{=}\frac{\mathrm{5}}{\mathrm{2}}{\left(\mathrm{1}\mathrm{+}\frac{\mathrm{19}\mathit{\mu }}{\mathrm{2}\mathit{g\rho R}}\right)}^{-1}& & \end{array}$(5)is the second Love number associated with the radial displacement, ρ the density, and μ the rigidity (or shear modulus). Darwin (1879) extended the results to a homogeneous incompressible viscous body and obtained the differential equation $\begin{array}{ccc}\mathit{\zeta }\mathrm{+}{\mathit{\tau }}_{\mathit{v}}\mathit{\zeta ̇}\mathrm{=}\mathrm{-}{\mathit{h}}_{\mathrm{f}}{\mathit{V}}_{\mathrm{p}}\mathit{/}\mathit{g,}& & \end{array}$(6)where h_f = 5/2 is the fluid Love number¹, $\begin{array}{ccc}{\mathit{\tau }}_{\mathit{v}}\mathrm{=}\frac{\mathrm{19}\mathit{\eta }}{\mathrm{2}\mathit{g\rho R}}\mathrm{=}\frac{\mathrm{38}\mathit{\pi }{\mathit{R}}^{\mathrm{2}}\mathit{\eta }}{\mathrm{3}\mathit{gm}}& & \end{array}$(7)is the fluid relaxation time, and η is the viscosity. Darwin (1879) also realized that viscoelastic deformations can easily be derived from expression (6) using the substitution $\begin{array}{ccc}\frac{\mathrm{1}}{\mathit{\eta }}\mathrm{\to }\frac{\mathrm{1}}{\mathit{\mu }}\left(\frac{\mathrm{1}}{{\mathit{\tau }}_{\mathrm{e}}}\mathrm{+}\frac{\mathrm{d}}{\mathrm{d}\mathit{t}}\right)\mathit{,}& & \end{array}$(8)where τ_e = η/μ is the elastic or Maxwell relaxation time. As a result, we obtain (Darwin 1879) $\begin{array}{ccc}\mathit{\zeta }\mathrm{+}\mathit{\tau }\mathit{\zeta ̇}\mathrm{=}\mathrm{-}{\mathit{h}}_{\mathrm{f}}\mathrm{(}{\mathit{V}}_{\mathrm{p}}\mathrm{+}{\mathit{\tau }}_{\mathrm{e}}\mathit{V̇}\mathrm{p}\mathrm{)}\mathit{/}\mathit{g,} \mathrm{with} \mathit{\tau }\mathrm{=}{\mathit{\tau }}_{\mathit{v}}\mathrm{+}{\mathit{\tau }}_{\mathrm{e}}\mathit{.}& & \end{array}$(9)Finally, the deformation ζ of the body generates an additional potential V′(r) given by $\begin{array}{ccc}{\mathit{V}}^{\mathrm{\prime }}\mathrm{\left(}{r}\mathrm{\right)}\mathrm{=}\mathrm{-}\mathrm{\int }\frac{\mathit{G\rho \zeta }\mathrm{\left(}{{R}}^{\mathrm{\prime }}\mathrm{\right)}}{\mathrm{|}{r}\mathrm{-}{{R}}^{\mathrm{\prime }}\mathrm{|}}\hspace{0.17em}\mathrm{d}{{R}}^{\mathrm{\prime }}\mathit{.}& & \end{array}$(10)In the case of quadrupolar deformations, as considered here (Eq. (1)), the additional potential exterior to the body (r ≥ R) is $\begin{array}{ccc}{\mathit{V}}^{\mathrm{\prime }}\mathrm{\left(}{r}\mathrm{\right)}\mathrm{=}\mathrm{-}\frac{\mathrm{3}}{\mathrm{5}}\mathit{g}{\left(\frac{\mathit{R}}{\mathit{r}}\right)}^{\mathrm{3}}\mathit{\zeta }\mathrm{\left(}{R}\mathrm{\right)}\mathit{.}& & \end{array}$(11)Thus, at the mean radius (r = R), the additional potential must satisfy $\begin{array}{ccc}{\mathit{V}}^{\mathrm{\prime }}\mathrm{+}\mathit{\tau }\mathit{V̇}\mathrm{\prime }\mathrm{=}{\mathit{k}}_{\mathrm{f}}\mathrm{(}{\mathit{V}}_{\mathrm{p}}\mathrm{+}{\mathit{\tau }}_{\mathrm{e}}\mathit{V̇}\mathrm{p}\mathrm{)}\mathit{,}& & \end{array}$(12)where k_f = 3/2 is the fluid second Love number for potential¹. The equilibrium potential V_e is obtained for static perturbations or instantaneous response (τ = τ_e = 0) $\begin{array}{ccc}{\mathit{V}}_{\mathrm{e}}\mathrm{=}{\mathit{k}}_{\mathrm{f}}{\mathit{V}}_{\mathrm{p}}\mathrm{=}\frac{{\mathit{k}}_{\mathrm{f}}}{\mathrm{2}}{\mathrm{\Omega }}^{\mathrm{2}}\mathrm{(}{K}\mathrm{·}{R}{\mathrm{)}}^{\mathrm{2}}\mathrm{-}{\mathit{k}}_{\mathrm{f}}\frac{\mathrm{3}\mathit{G}{\mathit{m}}_{\mathrm{0}}}{\mathrm{2}{\mathit{r}}^{\mathrm{3}}}\mathrm{(}{r̂}\mathrm{·}{R}{\mathrm{)}}^{\mathrm{2}}\mathit{.}& & \end{array}$(13)The instantaneous variation of the additional potential V′ that results from the perturbation V_p, can then be obtained from the equilibrium potential V_e through expression (12) as $\begin{array}{ccc}{\mathit{V}}^{\mathrm{\prime }}\mathrm{+}\mathit{\tau }\mathit{V̇}\mathrm{\prime }\mathrm{=}{\mathit{V}}_{\mathrm{e}}\mathrm{+}{\mathit{\tau }}_{\mathrm{e}}\mathit{V̇}\mathrm{e}\mathit{.}& & \end{array}$(14)We note that this viscoelastic deformation law is also fully compatible with a viscous rheology (Eq. (6)), provided that we set τ_e = 0. Previous expression could also have been obtained in the Fourier domain (see Appendix A).

2.2. Equations of motion

The force between the star and the planet is easily obtained from the gravitational potential of the system (Eq. (1)) as F = −m₀ × ∇V(r). From this force we can directly obtain the orbital evolution of the system $\begin{array}{c}{}_{{¨}}\\ {r}\end{array}\mathrm{=}{F}\mathit{/}\mathit{\beta }$, and the spin evolution $\begin{array}{c}{}_{\mathrm{¨}}\\ \mathit{\theta }\end{array}\mathrm{=}\mathrm{-}\mathrm{(}{r}\mathrm{\times }{F}\mathrm{)}\mathrm{·}{K}\mathit{/}\mathit{C}$ (using the conservation of the total angular momentum of the system), where β = m₀m/ (m₀ + m) is the reduced mass of the system, and C is the principal inertia moment of the planet. Expressing and I = (cosθ, sin θ), where f is the true longitude and θ is the rotation angle, we thus have $\begin{array}{ccc}\begin{array}{c}{¨}\\ {r}\end{array}\mathrm{=}& \mathrm{-}& \frac{{\mathit{\mu }}_{\mathrm{0}}}{{\mathit{r}}^{\mathrm{2}}}{r̂}\mathrm{-}\frac{\mathrm{3}{\mathit{\mu }}_{\mathrm{0}}{\mathit{R}}^{\mathrm{2}}}{\mathrm{2}{\mathit{r}}^{\mathrm{4}}}{\mathit{J}}_{\mathrm{2}}\hspace{0.17em}{r̂}\mathrm{-}\frac{\mathrm{9}{\mathit{\mu }}_{\mathrm{0}}{\mathit{R}}^{\mathrm{2}}}{{\mathit{r}}^{\mathrm{4}}}\left[{\mathit{C}}_{\mathrm{22}}\cos \mathrm{2}\mathit{\gamma }\mathrm{-}{\mathit{S}}_{\mathrm{22}}\sin \mathrm{2}\mathit{\gamma }\right]{r̂}\\ & \mathrm{+}& \end{array}$(15)and $\begin{array}{ccc}\begin{array}{c}\mathrm{¨}\\ \mathit{\theta }\end{array}& \mathrm{=}& \end{array}$(16)where μ₀ = G(m₀ + m), and γ = θ − f (Fig. 1).

Because the planet is not rigid and can be deformed under the action of a perturbing potential, the gravity field coefficients J₂, C₂₂ and S₂₂ are not constant. Indeed, according to the deformation law given by expression (14) we have

$\begin{array}{ccc}& & {\mathit{J}}_{\mathrm{2}}\mathrm{+}\mathit{\tau }\mathit{J̇}\mathrm{2}\mathrm{=}{\mathit{J}}_{\mathrm{2}}^{\mathrm{e}}\mathrm{+}{\mathit{\tau }}_{\mathrm{e}}\mathit{J̇}\begin{array}{c}\mathrm{e}\\ \mathrm{2}\end{array}\mathit{,}\\ & & {\mathit{C}}_{\mathrm{22}}\mathrm{+}\mathit{\tau }\mathit{Ċ}\mathrm{22}\mathrm{=}{\mathit{C}}_{\mathrm{22}}^{\mathrm{e}}\mathrm{+}{\mathit{\tau }}_{\mathrm{e}}\mathit{Ċ}\begin{array}{c}\mathrm{e}\\ \mathrm{22}\end{array}\mathit{,}\\ & & {\mathit{S}}_{\mathrm{22}}\mathrm{+}\mathit{\tau }\mathit{Ṡ}\mathrm{22}\mathrm{=}{\mathit{S}}_{\mathrm{22}}^{\mathrm{e}}\mathrm{+}{\mathit{\tau }}_{\mathrm{e}}\mathit{Ṡ}\begin{array}{c}\mathrm{e}\\ \mathrm{22}\end{array}\mathit{.}\end{array}$The equilibrium values for each coefficient are easily obtained comparing expressions (1) and (13), which gives (Correia & Rodríguez 2013) $\begin{array}{ccc}& & {\mathit{J}}_{\mathrm{2}}^{\mathrm{e}}\mathrm{=}{\mathit{k}}_{\mathrm{f}}\left[\frac{{\mathrm{\Omega }}^{\mathrm{2}}{\mathit{R}}^{\mathrm{3}}}{\mathrm{3}\mathit{Gm}}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}\frac{{\mathit{m}}_{\mathrm{0}}}{\mathit{m}}{\left(\frac{\mathit{R}}{\mathit{r}}\right)}^{\mathrm{3}}\right]\mathit{,}\\ & & {\mathit{C}}_{\mathrm{22}}^{\mathrm{e}}\mathrm{=}\frac{{\mathit{k}}_{\mathrm{f}}}{\mathrm{4}}\frac{{\mathit{m}}_{\mathrm{0}}}{\mathit{m}}{\left(\frac{\mathit{R}}{\mathit{r}}\right)}^{\mathrm{3}}\cos \mathrm{2}\mathit{\gamma ,}\\ & & {\mathit{S}}_{\mathrm{22}}^{\mathrm{e}}\mathrm{=}\mathrm{-}\frac{{\mathit{k}}_{\mathrm{f}}}{\mathrm{4}}\frac{{\mathit{m}}_{\mathrm{0}}}{\mathit{m}}{\left(\frac{\mathit{R}}{\mathit{r}}\right)}^{\mathrm{3}}\mathit{,}\end{array}$and for their time derivatives, respectively $\begin{array}{ccc}& & \mathit{J̇}\begin{array}{c}\mathrm{e}\\ \mathrm{2}\end{array}\mathrm{=}{\mathit{k}}_{\mathrm{f}}\left[\frac{\mathrm{2}\mathrm{\Omega }{\mathit{R}}^{\mathrm{3}}}{\mathrm{3}\mathit{Gm}}\mathrm{\Omega ̇}\mathrm{-}\frac{\mathrm{3}}{\mathrm{2}}\frac{{\mathit{m}}_{\mathrm{0}}}{\mathit{m}}{\left(\frac{\mathit{R}}{\mathit{r}}\right)}^{\mathrm{3}}\frac{{ṙ}\mathrm{·}{r}}{{\mathit{r}}^{\mathrm{2}}}\right]\mathit{,}\\ & & \mathit{Ċ}\begin{array}{c}\mathrm{e}\\ \mathrm{22}\end{array}\mathrm{=}\mathrm{-}\frac{{\mathit{k}}_{\mathrm{f}}}{\mathrm{4}}\frac{{\mathit{m}}_{\mathrm{0}}}{\mathit{m}}{\left(\frac{\mathit{R}}{\mathit{r}}\right)}^{\mathrm{3}}\left[\mathrm{3}\frac{\mathit{ṙ}}{\mathit{r}}\cos \mathrm{2}\mathit{\gamma }\mathrm{+}\mathrm{2}\mathit{\gamma ̇}\sin \mathrm{2}\mathit{\gamma }\right]\mathit{,}\\ & & \mathit{Ṡ}\begin{array}{c}\mathrm{e}\\ \mathrm{22}\end{array}\mathrm{=}\frac{{\mathit{k}}_{\mathrm{f}}}{\mathrm{4}}\frac{{\mathit{m}}_{\mathrm{0}}}{\mathit{m}}{\left(\frac{\mathit{R}}{\mathit{r}}\right)}^{\mathrm{3}}\left[\mathrm{3}\frac{\mathit{ṙ}}{\mathit{r}}\sin \mathrm{2}\mathit{\gamma }\mathrm{-}\mathrm{2}\mathit{\gamma ̇}\cos \mathrm{2}\mathit{\gamma }\right]\mathit{.}\end{array}$It should be noted that for the equilibrium solution (Eqs. (21), (22)), we have $\begin{array}{c}{}_{\mathrm{¨}}\\ \mathit{\theta }\end{array}\mathrm{=}\mathrm{0}$ (Eq. (16)).

3. General solution

Introducing the complex notations $\begin{array}{ccc}\mathit{Z}\mathrm{=}{\mathit{C}}_{\mathrm{22}}\mathrm{-}\mathrm{i}\hspace{0.17em}{\mathit{S}}_{\mathrm{22}} \mathrm{and} {\mathit{Z}}^{\mathrm{e}}\mathrm{=}{\mathit{C}}_{\mathrm{22}}^{\mathrm{e}}\mathrm{-}\mathrm{i}\hspace{0.17em}{\mathit{S}}_{\mathrm{22}}^{\mathrm{e}}\mathit{,}& & \end{array}$(26)we can rewrite Eqs. (18) and (19) together as $\begin{array}{ccc}\mathit{Z}\mathrm{+}\mathit{\tau }\mathit{Ż}\mathrm{=}{\mathit{Z}}^{\mathrm{e}}\mathrm{+}{\mathit{\tau }}_{\mathrm{e}}\dot{{\mathit{Z}}^{\mathrm{e}}}\mathit{.}& & \end{array}$(27)

The solution of Eq. (27) without the second hand member is Z = Z₀e^−t/τ, and the variation of the constant Z₀ gives the general solution of Eq. (27) $\begin{array}{ccc}\mathit{Z}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{=}\frac{{\mathit{\tau }}_{\mathrm{e}}}{\mathit{\tau }}\hspace{0.17em}{\mathit{Z}}^{\mathrm{e}}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{+}\frac{\mathrm{1}}{\mathit{\tau }}(\mathrm{1}\mathrm{-}\frac{{\mathit{\tau }}_{\mathrm{e}}}{\mathit{\tau }}){\mathrm{\int }}_{\mathrm{0}}^{\mathit{t}}{\mathit{Z}}^{\mathrm{e}}\mathrm{\left(}{\mathit{t}}^{\mathrm{\prime }}\mathrm{\right)}\hspace{0.17em}{\mathrm{e}}^{\mathrm{(}{\mathit{t}}^{\mathrm{\prime }}\mathrm{-}\mathit{t}\mathrm{)}\mathit{/}\mathit{\tau }}\mathrm{d}{\mathit{t}}^{\mathrm{\prime }}\mathrm{+}\mathrm{𝒞}\hspace{0.17em}{\mathrm{e}}^{\mathrm{-}\mathit{t}\mathit{/}\mathit{\tau }}\mathit{,}& & \end{array}$(28)where is an integration constant. Although one can search for a solution of this equation in the precise setting of this paper, it is much more efficient to solve it first in a very general setting where Z^e is a general almost periodic function, that is an expression of the form $\begin{array}{ccc}{\mathit{Z}}^{\mathrm{e}}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{=}\sum _{\mathit{k}}{\mathit{\beta }}_{\mathit{k}}\hspace{0.17em}{\mathrm{e}}^{\mathrm{i}{\mathit{\omega }}_{\mathit{k}}\mathit{t}}\mathit{,}& & \end{array}$(29)where β_k are constants, and ω_k are a numerable set of fixed frequencies. The result is $\begin{array}{ccc}\mathit{Z}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{=}\sum _{\mathit{k}}{\mathit{\beta }}_{\mathit{k}}\frac{\mathrm{1}\mathrm{+}\mathrm{i}{\mathit{\tau }}_{\mathrm{e}}{\mathit{\omega }}_{\mathit{k}}}{\mathrm{1}\mathrm{+}\mathrm{i}\mathit{\tau }{\mathit{\omega }}_{\mathit{k}}}\hspace{0.17em}{\mathrm{e}}^{\mathrm{i}{\mathit{\omega }}_{\mathit{k}}\mathit{t}}\mathrm{+}\mathrm{𝒞}\hspace{0.17em}{\mathrm{e}}^{\mathrm{-}\mathit{t}\mathit{/}\mathit{\tau }}\mathit{.}& & \end{array}$(30)The first term is the stationnary mode and the second one the transient mode that will decay rapidly to zero with time. In general, the constant in Eq. (30) is different from that of Eq. (28). It should be noted that here we do not have small divisors problems as the denominators 1 + iτω_k never vanish. A constant term β_k0 in Z^e (corresponding to ω_k0 = 0) will add an equivalent constant contribution β_k0 to Z.

4. Tidal evolution

From the general solution of Eq. (27), we can study more specific solutions. For this, we will look to the explicit expression of Z^e. Using expressions (21) and (22), we can write $\begin{array}{ccc}{\mathit{Z}}^{\mathrm{e}}\mathrm{=}\mathrm{𝒜}{\left(\frac{\mathit{a}}{\mathit{r}}\right)}^{\mathrm{3}}{\mathrm{e}}^{\mathrm{i}\mathrm{2}\mathrm{(}\begin{array}{c}\mathrm{˜}\\ \mathit{\theta }\end{array}\mathrm{-}\mathit{v}\mathrm{)}}\mathit{,} \mathrm{with} \mathrm{𝒜}\mathrm{=}\frac{{\mathit{k}}_{\mathrm{f}}}{\mathrm{4}}\frac{{\mathit{m}}_{\mathrm{0}}}{\mathit{m}}{\left(\frac{\mathit{R}}{\mathit{a}}\right)}^{\mathrm{3}}\mathit{,}& & \end{array}$(31)where a is the semi-major axis, $\begin{array}{c}{}_{\mathrm{˜}}\\ \mathit{\theta }\end{array}\mathrm{=}\mathit{\theta }\mathrm{-}\mathit{\varpi }$, ϖ is the longitude of the pericenter, and v is the true anomaly (Fig. 1). Using the expansion of (a/r)³exp(−2iv) in terms of Hansen coefficients ${\mathit{X}}_{\mathit{k}}^{\mathit{l,m}}\mathrm{\left(}\mathit{e}\mathrm{\right)}$ (e.g., Hughes 1981) (see also Appendix C), $\begin{array}{ccc}{\left(\frac{\mathit{r}}{\mathit{a}}\right)}^{\mathit{l}}{\mathrm{e}}^{\mathrm{i}\mathit{mv}}\mathrm{=}\sum _{\mathit{k}\mathrm{=}\mathrm{-}\mathrm{\infty }}^{\mathrm{+}\mathrm{\infty }}{\mathit{X}}_{\mathit{k}}^{\mathit{l,m}}\mathrm{\left(}\mathit{e}\mathrm{\right)}\hspace{0.17em}{\mathrm{e}}^{\mathrm{i}\mathit{kM}}\mathit{,}& & \end{array}$(32)we obtain $\begin{array}{ccc}{\mathit{Z}}^{\mathrm{e}}\mathrm{=}\sum _{\mathit{k}\mathrm{=}\mathrm{-}\mathrm{\infty }}^{\mathrm{+}\mathrm{\infty }}{\mathit{\beta }}_{\mathit{k}}\hspace{0.17em}{\mathrm{e}}^{\mathrm{i}\mathrm{(}\mathrm{2}\begin{array}{c}\mathrm{˜}\\ \mathit{\theta }\end{array}\mathrm{-}\mathit{kM}\mathrm{)}}\mathit{,} \mathrm{with} {\mathit{\beta }}_{\mathit{k}}\mathrm{=}\mathrm{𝒜}\hspace{0.17em}{\mathit{X}}_{\mathit{k}}^{\mathrm{-}\mathrm{3}\mathit{,}\mathrm{2}}\mathrm{\left(}\mathit{e}\mathrm{\right)}\mathit{,}& & \end{array}$(33)where M is the mean anomaly, and e is the eccentricity. β_k depends only on the orbital parameters that are assumed to be constant over one revolution. The evolution of θ is given by Eq. (16), but since tidal effects are usually very weak, in a zero order solution, we can also neglect $\begin{array}{c}{\mathrm{¨}}_{}\\ \mathit{\theta }\end{array}$, and assume that the rotation rate $\mathrm{\Omega }\mathrm{=}\mathrm{d}\begin{array}{c}{}_{\mathrm{˜}}\\ \mathit{\theta }\end{array}\mathit{/}\mathrm{d}\mathit{t}$ is constant over one orbital revolution. Thus, we can write $\begin{array}{ccc}{\mathit{\omega }}_{\mathit{k}}\mathrm{=}\mathrm{2}\mathrm{\Omega }\mathrm{-}\mathit{kn}\mathrm{=}\frac{\mathrm{d}}{\mathrm{d}\mathit{t}}\mathrm{(}\mathrm{2}\begin{array}{c}\mathrm{˜}\\ \mathit{\theta }\end{array}\mathrm{-}\mathit{kM}\mathrm{)}\mathit{,}& & \end{array}$(34)where n = Ṁ, and ω_k is a constant frequency. As the β_k are constant in Eq. (29), we are in a special case of the application of the results of Sect. 3, and we obtain, after neglecting the transtient part of expression (30) that decays to zero for t ≫ τ, $\begin{array}{ccc}\mathit{Z}\mathrm{=}\sum _{\mathit{k}\mathrm{=}\mathrm{-}\mathrm{\infty }}^{\mathrm{+}\mathrm{\infty }}{\mathit{\beta }}_{\mathit{k}}\frac{\mathrm{1}\mathrm{+}\mathrm{i}{\mathit{\tau }}_{\mathrm{e}}{\mathit{\omega }}_{\mathit{k}}}{\mathrm{1}\mathrm{+}\mathrm{i}\mathit{\tau }{\mathit{\omega }}_{\mathit{k}}}\hspace{0.17em}{\mathrm{e}}^{\mathrm{i}\mathrm{(}\mathrm{2}\begin{array}{c}\mathrm{˜}\\ \mathit{\theta }\end{array}\mathrm{-}\mathit{kM}\mathrm{)}}\mathit{.}& & \end{array}$(35)

Fig. 2 Evolution of the second Love number k2, the dissipation Q-factor, and phase lag k2sinδ as a function of the tidal frequency ω (normalized by the relaxation time τ). For all parameters there is a transition of regime around τω ~ 1. We use here the Earth values for the Love numbers, ke = 0.3 and kf = 0.9 (Yoder 1995).

4.1. Tidal torque

We have already seen that for Z = Z^e, $\begin{array}{c}{}_{\mathrm{¨}}\\ \mathit{\theta }\end{array}\mathrm{=}\mathrm{0}$, as there are no coupling torques in the problem due to alignment of the planet deformation with the star. It is then useful to use the difference $\begin{array}{ccc}\mathrm{\Delta }\mathit{Z}\mathrm{=}\mathit{Z}\mathrm{-}{\mathit{Z}}^{\mathrm{e}}\mathrm{=}\mathrm{-}\left(\mathrm{1}\mathrm{-}\frac{{\mathit{\tau }}_{\mathrm{e}}}{\mathit{\tau }}\right){\mathrm{e}}^{\mathrm{i}\mathrm{2}\begin{array}{c}\mathrm{˜}\\ \mathit{\theta }\end{array}}\sum _{\mathit{k}\mathrm{=}\mathrm{-}\mathrm{\infty }}^{\mathrm{+}\mathrm{\infty }}{\mathit{\beta }}_{\mathit{k}}\frac{\mathrm{i}\mathit{\tau }{\mathit{\omega }}_{\mathit{k}}}{\mathrm{1}\mathrm{+}\mathrm{i}\mathit{\tau }{\mathit{\omega }}_{\mathit{k}}}\hspace{0.17em}{\mathrm{e}}^{\mathrm{-}\mathrm{i}\mathit{kM}}\mathit{.}& & \end{array}$(36)Denoting $\begin{array}{ccc}\mathrm{ℬ}\mathrm{=}\frac{\mathrm{6}\mathit{G}{\mathit{m}}_{\mathrm{0}}\mathit{m}{\mathit{R}}^{\mathrm{2}}}{\mathit{C}{\mathit{a}}^{\mathrm{3}}}\mathit{,}& & \end{array}$(37)we obtain the variations of $\begin{array}{c}{\mathrm{¨}}_{}\\ \mathit{\theta }\end{array}$ in the first order solution $\begin{array}{ccc}\begin{array}{c}\mathrm{¨}\\ \mathit{\theta }\end{array}\mathrm{=}\mathrm{ℬ}\left(\frac{{\mathit{a}}^{\mathrm{3}}}{{\mathit{r}}^{\mathrm{3}}}\right)\mathrm{\Im }\left(\mathrm{\Delta }\mathit{Z}\hspace{0.17em}{\mathrm{e}}^{\mathrm{-}\mathrm{i}\mathrm{2}\mathrm{(}\begin{array}{c}\mathrm{˜}\\ \mathit{\theta }\end{array}\mathrm{-}\mathit{v}\mathrm{)}}\right)\mathrm{=}\frac{\mathrm{ℬ}}{\mathrm{𝒜}}\mathrm{\Im }\mathrm{(}\mathrm{\Delta }\mathit{Z}\hspace{0.17em}{{\mathit{Z}}^{\mathrm{e}}}^{\mathit{\star }}\mathrm{)}\mathit{.}& & \end{array}$(38)where ℑ(z) denotes the imaginary part of z and z^⋆ its complex conjugate. We thus see from expression (36) that θ disappears from the expression of $\begin{array}{c}{\mathrm{¨}}_{}\\ \mathit{\theta }\end{array}$. After an expansion in Hansen coefficients of expression (38), we obtain an explicit expression of $\begin{array}{c}{\mathrm{¨}}_{}\\ \mathit{\theta }\end{array}$ with respect to the mean anomaly M$\begin{array}{ccc}\begin{array}{c}\mathrm{¨}\\ \mathit{\theta }\end{array}& \mathrm{=}& \mathrm{-}\mathrm{𝒦}\hspace{0.17em}\mathrm{\Im }\left(\left(\sum _{\mathit{k}\mathrm{=}\mathrm{-}\mathrm{\infty }}^{\mathrm{+}\mathrm{\infty }}{\mathit{X}}_{\mathit{k}}^{\mathrm{-}\mathrm{3}\mathit{,}\mathrm{2}}\frac{\mathrm{i}\mathit{\tau }{\mathit{\omega }}_{\mathit{k}}}{\mathrm{1}\mathrm{+}\mathrm{i}\mathit{\tau }{\mathit{\omega }}_{\mathit{k}}}{\mathrm{e}}^{\mathrm{-}\mathrm{i}\mathit{kM}}\right)\left(\sum _{\mathit{l}\mathrm{=}\mathrm{-}\mathrm{\infty }}^{\mathrm{+}\mathrm{\infty }}{\mathit{X}}_{\mathit{l}}^{\mathrm{-}\mathrm{3}\mathit{,}\mathrm{2}}{\mathrm{e}}^{\mathrm{i}\mathit{lM}}\right)\right)\\ & \mathrm{=}& \end{array}$(39)where $\begin{array}{ccc}\mathrm{𝒦}\mathrm{=}\mathrm{𝒜ℬ}\left(\mathrm{1}\mathrm{-}\frac{{\mathit{\tau }}_{\mathrm{e}}}{\mathit{\tau }}\right)\mathrm{=}{\mathit{k}}_{\mathrm{f}}\frac{\mathrm{3}\mathit{G}{\mathit{m}}_{\mathrm{0}}^{\mathrm{2}}{\mathit{R}}^{\mathrm{5}}}{\mathrm{2}\mathit{C}{\mathit{a}}^{\mathrm{6}}}\left(\mathrm{1}\mathrm{-}\frac{{\mathit{\tau }}_{\mathrm{e}}}{\mathit{\tau }}\right)\mathrm{·}& & \end{array}$(40)The first order secular evolution of $\begin{array}{c}{\mathrm{¨}}_{}\\ \mathit{\theta }\end{array}$ over one orbital period is easily obtained by averaging expression (39) over M, which is equivalent to keep only the terms with j = 0. We thus get $\mathrm{⟨}\begin{array}{c}{}_{\mathrm{¨}}\\ \mathit{\theta }\end{array}{\mathrm{⟩}}_{\mathit{M}}\mathrm{=}\mathit{T}$, with the secular torque T given by $\begin{array}{ccc}\mathit{T}\mathrm{=}\mathrm{-}\mathrm{𝒦}\sum _{\mathit{k}\mathrm{=}\mathrm{-}\mathrm{\infty }}^{\mathrm{+}\mathrm{\infty }}{\mathrm{(}{{\mathit{X}}_{\mathit{k}}^{\mathrm{-}\mathrm{3}\mathit{,}\mathrm{2}}}^{\mathrm{)}}}^{\mathrm{2}}\frac{\mathit{\tau }{\mathit{\omega }}_{\mathit{k}}}{\mathrm{1}\mathrm{+}{\mathit{\tau }}^{\mathrm{2}}{\mathit{\omega }}_{\mathit{k}}^{\mathrm{2}}}\mathrm{·}& & \end{array}$(41)

4.2. Relation with classical notations

Expressions (35) and (41) can be related to classical tidal parameters, by transforming the complex notation in a real amplitude and a phase lag $\begin{array}{ccc}\mathit{Z}\mathrm{=}\frac{\mathrm{1}}{\mathrm{4}}\frac{{\mathit{m}}_{\mathrm{0}}}{\mathit{m}}{\left(\frac{\mathit{R}}{\mathit{a}}\right)}^{\mathrm{3}}\sum _{\mathit{k}\mathrm{=}\mathrm{-}\mathrm{\infty }}^{\mathrm{+}\mathrm{\infty }}{\mathit{X}}_{\mathit{k}}^{\mathrm{-}\mathrm{3}\mathit{,}\mathrm{2}}\mathrm{\left(}\mathit{e}\mathrm{\right)}\hspace{0.17em}{\mathit{k}}_{\mathrm{2}}\mathrm{\left(}{\mathit{\omega }}_{\mathit{k}}\mathrm{\right)}\hspace{0.17em}{\mathrm{e}}^{\mathrm{i}\mathrm{(}\mathrm{2}\begin{array}{c}\mathrm{˜}\\ \mathit{\theta }\end{array}\mathrm{-}\mathit{kM}\mathrm{-}{\mathit{\delta }}_{\mathit{k}}\mathrm{)}}\mathit{,}& & \end{array}$(42)and $\begin{array}{ccc}\mathit{T}\mathrm{=}\mathrm{-}\frac{\mathrm{3}\mathit{G}{\mathit{m}}_{\mathrm{0}}^{\mathrm{2}}{\mathit{R}}^{\mathrm{5}}}{\mathrm{2}\mathit{C}{\mathit{a}}^{\mathrm{6}}}\sum _{\mathit{k}\mathrm{=}\mathrm{-}\mathrm{\infty }}^{\mathrm{+}\mathrm{\infty }}{\mathrm{(}{{\mathit{X}}_{\mathit{k}}^{\mathrm{-}\mathrm{3}\mathit{,}\mathrm{2}}}^{\mathrm{)}}}^{\mathrm{2}}{\mathit{k}}_{\mathrm{2}}\mathrm{\left(}{\mathit{\omega }}_{\mathit{k}}\mathrm{\right)}\sin {\mathit{\delta }}_{\mathit{k}}\mathit{,}& & \end{array}$(43)with $\begin{array}{ccc}\tan {\mathit{\delta }}_{\mathit{k}}\mathrm{=}\frac{\mathrm{(}\mathit{\tau }\mathrm{-}{\mathit{\tau }}_{\mathrm{e}}\mathrm{)}{\mathit{\omega }}_{\mathit{k}}}{\mathrm{1}\mathrm{+}\mathit{\tau }{\mathit{\tau }}_{\mathrm{e}}{\mathit{\omega }}_{\mathit{k}}^{\mathrm{2}}}\mathit{,}& & \end{array}$(44)and $\begin{array}{ccc}{\mathit{k}}_{\mathrm{2}}\mathrm{\left(}{\mathit{\omega }}_{\mathit{k}}\mathrm{\right)}\mathrm{=}{\mathit{k}}_{\mathrm{f}}\sqrt{\frac{\mathrm{1}\mathrm{+}\mathrm{(}{\mathit{\tau }}_{\mathrm{e}}{\mathit{\omega }}_{\mathit{k}}{\mathrm{)}}^{\mathrm{2}}}{\mathrm{1}\mathrm{+}\mathrm{(}\mathit{\tau }{\mathit{\omega }}_{\mathit{k}}{\mathrm{)}}^{\mathrm{2}}}}\mathrm{·}& & \end{array}$(45)We then conclude that the C₂₂ and S₂₂ are shifted by an angle δ_k with respect to the perturbation, and weakened by a factor of k₂(ω_k) /k_f (Fig. 2). The phase lag δ_k (Eq. (44)) is zero for ω_k = 0 and ω_k → ∞, and its maximal value is obtained at the frequency ω_k = (ττ_e)^−1/2. The amplitude k₂(ω_k) is monotonically decreasing with ω_k with $\begin{array}{ccc}{\mathit{k}}_{\mathrm{f}}\mathrm{=}{\mathit{k}}_{\mathrm{2}}\mathrm{\left(}\mathrm{0}\mathrm{\right)}\mathrm{\ge }{\mathit{k}}_{\mathrm{2}}\mathrm{\left(}{\mathit{\omega }}_{\mathit{k}}\mathrm{\right)}\mathrm{\ge }{\mathit{k}}_{\mathrm{2}}\mathrm{\left(}\mathrm{\infty }\mathrm{\right)}\mathrm{=}{\mathit{k}}_{\mathrm{f}}\frac{{\mathit{\tau }}_{\mathrm{e}}}{\mathit{\tau }}\mathrm{\equiv }{\mathit{k}}_{\mathrm{e}}\mathit{.}& & \end{array}$(46)We call k_e the elastic Love number. One can consider that the transition between the elastic and the viscous regimes occurs for k₂(ω₀) = (k_f + k_e)/2. Because of the expression of k₂ (Eq. (45)), we use ${\mathit{k}}_{\mathrm{2}}\mathrm{\left(}{\mathit{\omega }}_{\mathrm{0}}\mathrm{\right)}\mathrm{=}\sqrt{\mathrm{(}{\mathit{k}}_{\mathrm{f}}^{\mathrm{2}}\mathrm{+}{\mathit{k}}_{\mathrm{e}}^{\mathrm{2}}\mathrm{)}\mathit{/}\mathrm{2}}$, which gives the very simple transition condition $\begin{array}{ccc}\mathit{\tau }{\mathit{\omega }}_{\mathrm{0}}\mathrm{=}\mathrm{1.}& & \end{array}$(47)If one is able to measure k_e, k_f and δ_k for a planet, it becomes possible to determine the relaxation times τ and τ_e.

4.3. Energy dissipation and orbital evolution

Consider that , where Ω is constant and δΩ a small perturbation of Ω. Neglecting terms of second order, the rotational energy dissipated in the planet is then easily obtained as $\begin{array}{ccc}\mathrm{⟨}\mathit{Ė}\mathrm{rot}{\mathrm{⟩}}_{\mathit{M}}\mathrm{=}{⟨\frac{\mathrm{d}}{\mathrm{d}\mathit{t}}\left(\frac{\mathrm{1}}{\mathrm{2}}\mathit{C}\hspace{0.17em}\mathit{\theta ̇}\mathrm{2}\right)⟩}_{\mathit{M}}\mathrm{\approx }\mathit{C}\hspace{0.17em}\mathrm{\Omega }\hspace{0.17em}\mathrm{⟨}\begin{array}{c}\mathrm{¨}\\ \mathit{\theta }\end{array}{\mathrm{⟩}}_{\mathit{M}}\mathrm{=}\mathit{C}\hspace{0.17em}\mathrm{\Omega }\hspace{0.17em}\mathit{T}\mathit{.}& & \end{array}$(48)The variation of orbital energy is obtained as the variation of the potential energy (Eq. (1)) under the variation of the potential coefficients, that is $\begin{array}{ccc}\mathit{Ė}\mathrm{orb}\mathrm{=}{\mathit{m}}_{\mathrm{0}}\frac{\mathit{\partial V}}{\mathit{\partial t}}\mathrm{=}\mathrm{-}\frac{\mathit{G}{\mathit{m}}_{\mathrm{0}}\mathit{m}{\mathit{R}}^{\mathrm{2}}}{\mathrm{2}{\mathit{a}}^{\mathrm{3}}}\left[\frac{{\mathit{a}}^{\mathrm{3}}}{{\mathit{r}}^{\mathrm{3}}}\mathit{J̇}\mathrm{2}\mathrm{+}\mathrm{6}\frac{{\mathit{a}}^{\mathrm{3}}}{{\mathit{r}}^{\mathrm{3}}}\hspace{0.17em}\mathrm{\Re }\left(\mathit{Ż}{\mathrm{e}}^{\mathrm{i}\mathrm{2}\mathrm{(}\mathit{f}\mathrm{-}\begin{array}{c}\mathrm{˜}\\ \mathit{\theta }\end{array}\mathrm{)}}\right)\right]\mathit{,}& & \end{array}$(49)where ℜ(z) is the real part of the complex number z. Using the expression of ${\mathit{J}}_{\mathrm{2}}^{\mathrm{e}}$ (Eq. (20)) and Sect. 3, we have $\begin{array}{ccc}\mathit{J̇}\mathrm{2}\mathrm{=}\mathrm{2}\mathrm{i𝒜}\sum _{\mathit{k}\mathrm{=}\mathrm{-}\mathrm{\infty }}^{\mathrm{+}\mathrm{\infty }}{\mathit{X}}_{\mathit{k}}^{\mathrm{-}\mathrm{3}\mathit{,}\mathrm{0}}\mathit{kn}\frac{\mathrm{1}\mathrm{+}\mathrm{i}{\mathit{\tau }}_{\mathrm{e}}\mathit{kn}}{\mathrm{1}\mathrm{+}\mathrm{i}\mathit{\tau kn}}{\mathrm{e}}^{\mathit{ikM}}\mathit{,}& & \end{array}$(50)and from expression (35) $\begin{array}{ccc}\frac{{\mathit{a}}^{\mathrm{3}}}{{\mathit{r}}^{\mathrm{3}}}\mathit{Ż}{\mathrm{e}}^{\mathrm{i}\mathrm{2}\mathrm{(}\mathit{f}\mathrm{-}\begin{array}{c}\mathrm{˜}\\ \mathit{\theta }\end{array}\mathrm{)}}\mathrm{=}\mathrm{i𝒜}\left(\sum _{\mathit{k}\mathrm{=}\mathrm{-}\mathrm{\infty }}^{\mathrm{+}\mathrm{\infty }}{\mathit{X}}_{\mathit{k}}^{\mathrm{-}\mathrm{3}\mathit{,}\mathrm{2}}{\mathit{\omega }}_{\mathit{k}}\frac{\mathrm{1}\mathrm{+}\mathrm{i}{\mathit{\tau }}_{\mathrm{e}}{\mathit{\omega }}_{\mathit{k}}}{\mathrm{1}\mathrm{+}\mathrm{i}\mathit{\tau }{\mathit{\omega }}_{\mathit{k}}}{\mathrm{e}}^{\mathrm{-}\mathrm{i}\mathit{kM}}\right)\left(\frac{{\mathit{a}}^{\mathrm{3}}}{{\mathit{r}}^{\mathrm{3}}}{\mathrm{e}}^{\mathrm{i}\mathrm{2}\mathit{v}}\right)\mathit{.}& & \end{array}$(51)Finally, after averaging over the mean anomaly M, $\begin{array}{ccc}\mathrm{⟨}\mathit{Ė}\mathrm{orb}{\mathrm{⟩}}_{\mathit{M}}\mathrm{=}\frac{\mathit{C}\mathrm{𝒦}}{\mathrm{6}}\sum _{\mathit{k}\mathrm{=}\mathrm{-}\mathrm{\infty }}^{\mathrm{+}\mathrm{\infty }}\left[{\mathrm{(}{{\mathit{X}}_{\mathit{k}}^{\mathrm{-}\mathrm{3}\mathit{,}\mathrm{0}}}^{\mathrm{)}}}^{\mathrm{2}}\frac{\mathit{\tau }{\mathit{k}}^{\mathrm{2}}{\mathit{n}}^{\mathrm{2}}}{\mathrm{1}\mathrm{+}{\mathit{\tau }}^{\mathrm{2}}{\mathit{k}}^{\mathrm{2}}{\mathit{n}}^{\mathrm{2}}}\mathrm{+}{\mathrm{(}{{\mathit{X}}_{\mathit{k}}^{\mathrm{-}\mathrm{3}\mathit{,}\mathrm{2}}}^{\mathrm{)}}}^{\mathrm{2}}\frac{\mathrm{3}\mathit{\tau }{\mathit{\omega }}_{\mathit{k}}^{\mathrm{2}}}{\mathrm{1}\mathrm{+}{\mathit{\tau }}^{\mathrm{2}}{\mathit{\omega }}_{\mathit{k}}^{\mathrm{2}}}\right]\mathrm{·}& & \end{array}$(52)The last term results form the contribution of C₂₂ and S₂₂, while the first term results from the contribution of the J₂, so it does not depend on the rotation angle (Eq. (20)). We can also easily compute the semi-major axis and eccentricity evolution from the tidal energy (e.g., Correia et al. 2011) $\begin{array}{ccc}\mathit{ȧ}\mathrm{=}\frac{\mathrm{2}\mathit{Ė}\mathrm{orb}}{\mathit{\beta a}{\mathit{n}}^{\mathrm{2}}}\mathit{,} \mathrm{and} \mathit{ė}\mathrm{=}\frac{\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{e}}^{\mathrm{2}}\mathrm{)}}{\mathit{\beta n}{\mathit{a}}^{\mathrm{2}}\mathit{e}}\left(\frac{\mathit{Ė}\mathrm{orb}}{\mathit{n}}\mathrm{+}\frac{\mathit{Ė}\mathrm{rot}}{\mathrm{\Omega }\sqrt{\mathrm{1}\mathrm{-}{\mathit{e}}^{\mathrm{2}}}}\right)\mathrm{·}& & \end{array}$(53)The total energy released inside the body due to tides is then $\begin{array}{ccc}\mathit{Ė}\mathrm{=}\mathrm{-}\mathrm{(}\mathit{Ė}\mathrm{rot}\mathrm{+}\mathit{Ė}\mathrm{orb}\mathrm{)}\mathit{.}& & \end{array}$(54)A commonly used dimensionless measure of the tidal dissipation is the quality factor, $\begin{array}{ccc}{\mathit{Q}}^{-1}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}\mathit{\pi }{\mathit{E}}_{\mathrm{0}}}\mathrm{􏽉}\mathit{Ė}\hspace{0.17em}\mathrm{d}\mathit{t,}& & \end{array}$(55)where the line integral over Ė is the energy dissipated during one period of tidal stress, and E₀ is the peak energy stored in the system during the same period. In general, Q is a function of the frequency ω_k, and it can be related to the phase lag through ${\mathit{Q}}_{\mathit{k}}^{-1}\mathrm{=}\sin {\mathit{\delta }}_{\mathit{k}}$ (Efroimsky 2012, Eq. (141)). In that case, the ratio k₂/Q_k is directly connected with the tidal torque (Eq. (43)) and therefore it is easily measurable $\begin{array}{ccc}{\mathit{k}}_{\mathrm{2}}\mathrm{\left(}{\mathit{\omega }}_{\mathit{k}}\mathrm{\right)}\mathit{/}{\mathit{Q}}_{\mathit{k}}\mathrm{=}{\mathit{k}}_{\mathrm{2}}\mathrm{\left(}{\mathit{\omega }}_{\mathit{k}}\mathrm{\right)}\sin {\mathit{\delta }}_{\mathit{k}}\mathrm{=}{\mathit{k}}_{\mathrm{f}}\frac{\mathrm{(}\mathit{\tau }\mathrm{-}{\mathit{\tau }}_{\mathrm{e}}\mathrm{)}{\mathit{\omega }}_{\mathit{k}}}{\mathrm{1}\mathrm{+}\mathrm{(}\mathit{\tau }{\mathit{\omega }}_{\mathit{k}}{\mathrm{)}}^{\mathrm{2}}}\mathrm{·}& & \end{array}$(56)For static perturbations (ω_k = 0) the deformation of the planet is maximal (Eq. (45)), but the phase lag δ_k is zero, so there is no tidal dissipation. All gravity coefficients are given by their equilibrium values (Eqs. (20)−(22)), and the last term in expression (15) and expression (16) become zero. Then, the only effect of the deformation is to add some precession to the elliptical orbit. The same is true for short period perturbations (ω_k ≫ 1 /τ), since we also have k₂/Q_k → 0 (Fig. 2c). However, in this case the equilibrium figure must be corrected using k_e instead of k_f (Eq. (46)).

5. Equilibrium configurations

5.1. Low frequency regime (| τω_k | ≪ 1 for all k)

When τω_k ≪ 1 we have $\mathrm{1}\mathrm{+}{\mathit{\tau }}^{\mathrm{2}}{\mathit{\omega }}_{\mathit{k}}^{\mathrm{2}}\mathrm{\approx }\mathrm{1}$, and the expression of the secular evolution of the rotation rate (41) simplifies. With ω_k = 2Ω − kn, and using Appendix B we obtain the following expression for the secular tidal torque $\begin{array}{ccc}\mathit{T}\mathrm{=}\mathrm{-}\frac{\mathrm{3}\mathit{G}{\mathit{m}}_{\mathrm{0}}^{\mathrm{2}}{\mathit{R}}^{\mathrm{5}}}{\mathit{C}{\mathit{a}}^{\mathrm{6}}}{\mathit{k}}_{\mathrm{f}}\hspace{0.17em}\mathit{n}\mathrm{\Delta }\mathit{\tau }\left({\mathit{f}}_{\mathrm{1}}\mathrm{\left(}\mathit{e}\mathrm{\right)}\frac{\mathrm{\Omega }}{\mathit{n}}\mathrm{-}{\mathit{f}}_{\mathrm{2}}\mathrm{\left(}\mathit{e}\mathrm{\right)}\right)\mathit{,}& & \end{array}$(57)where Δτ = τ − τ_e = τ_v is a constant time-lag between the perturbation and the maximal deformation, that is equivalent to the fluid or “viscous” relaxation time (Eq. (9)), $\begin{array}{ccc}{\mathit{f}}_{\mathrm{1}}\mathrm{\left(}\mathit{e}\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}\mathrm{+}\mathrm{3}{\mathit{e}}^{\mathrm{2}}\mathrm{+}\frac{\mathrm{3}}{\mathrm{8}}{\mathit{e}}^{\mathrm{4}}}{\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{e}}^{\mathrm{2}}{\mathrm{)}}^{\mathrm{9}\mathit{/}\mathrm{2}}}\mathit{,}& & \end{array}$(58)and $\begin{array}{ccc}{\mathit{f}}_{\mathrm{2}}\mathrm{\left(}\mathit{e}\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}\mathrm{+}\frac{\mathrm{15}}{\mathrm{2}}{\mathit{e}}^{\mathrm{2}}\mathrm{+}\frac{\mathrm{45}}{\mathrm{8}}{\mathit{e}}^{\mathrm{4}}\mathrm{+}\frac{\mathrm{5}}{\mathrm{16}}{\mathit{e}}^{\mathrm{6}}}{\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{e}}^{\mathrm{2}}{\mathrm{)}}^{\mathrm{6}}}\mathrm{·}& & \end{array}$(59)Thus, for a given eccentricity, there is only one possible equilibrium configuration for the rotation rate, obtained when T = 0, $\begin{array}{ccc}{\mathrm{\Omega }}_{\mathrm{e}}\mathit{/}\mathit{n}\mathrm{=}{\mathit{f}}_{\mathrm{2}}\mathrm{\left(}\mathit{e}\mathrm{\right)}\mathit{/}{\mathit{f}}_{\mathrm{1}}\mathrm{\left(}\mathit{e}\mathrm{\right)}\mathit{,}& & \end{array}$(60)which is often known as the pseudo-synchronization. The low frequency regime is widely described in the literature, and it is usually known as the “viscous” or constant time-lag approximation (e.g., Singer 1968; Alexander 1973; Mignard 1979).

5.2. High frequency regime (τn ≳ 1)

In this case, the situation for equilibrium points (T = 0) is more complex than in the low frequency regime as the expression (41) may have many zeros. The situation is very similar to that discussed for the Andrade rheology model (e.g., Makarov & Efroimsky 2013). Indeed, the torque (Eq. (41)) is a sum $\begin{array}{ccc}\mathit{T}\mathrm{=}\mathrm{𝒦}\sum _{\mathit{k}\mathrm{=}\mathrm{-}\mathrm{\infty }}^{\mathrm{\infty }}{\mathrm{(}{{\mathit{X}}_{\mathit{k}}^{\mathrm{-}\mathrm{3}\mathit{,}\mathrm{2}}}^{\mathrm{)}}}^{\mathrm{2}}{\mathit{T}}_{\mathit{k}}^{\mathrm{\left(}\mathit{M}\mathrm{\right)}}& & \end{array}$(61)with components of the form $\begin{array}{ccc}{\mathit{T}}_{\mathit{k}}^{\mathrm{\left(}\mathit{M}\mathrm{\right)}}\mathrm{=}\frac{\mathrm{-}\mathit{\tau }{\mathit{\omega }}_{\mathit{k}}}{\mathrm{1}\mathrm{+}{\mathit{\tau }}^{\mathrm{2}}{\mathit{\omega }}_{\mathit{k}}^{\mathrm{2}}}\mathrm{·}& & \end{array}$(62)With Andrade rheology, the torque has an equivalent expression, where all ${\mathit{T}}_{\mathit{k}}^{\mathrm{\left(}\mathit{M}\mathrm{\right)}}$ are replaced by ${\mathit{T}}_{\mathit{k}}^{\mathrm{\left(}\mathit{A}\mathrm{\right)}}$ defined as (Efroimsky 2012) $\begin{array}{ccc}{\mathit{T}}_{\mathit{k}}^{\mathrm{\left(}\mathit{A}\mathrm{\right)}}\mathrm{=}\frac{\mathrm{ℐ}}{{\mathrm{ℛ}}^{\mathrm{2}}\mathrm{+}{\mathrm{ℐ}}^{\mathrm{2}}}\mathrm{sign}\mathrm{\left(}\mathit{\tau }{\mathit{\omega }}_{\mathit{k}}\mathrm{\right)}& & \end{array}$(63)with $\begin{array}{ccc}\mathrm{ℛ}& \mathrm{=}& \mathrm{1}\mathrm{+}\frac{\mathrm{1}}{\mathrm{|}\mathit{\tau }{\mathit{\omega }}_{\mathit{k}}{\mathrm{|}}^{\mathit{\alpha }}}{\left(\frac{{\mathit{\tau }}_{\mathrm{e}}}{\mathit{\tau }}\right)}^{\mathrm{1}\mathrm{-}\mathit{\alpha }}\cos \left(\frac{\mathit{\alpha \pi }}{\mathrm{2}}\right)\hspace{0.17em}\mathrm{\Gamma }\mathrm{(}\mathit{\alpha }\mathrm{+}\mathrm{1}\mathrm{)}\mathit{,}\\ \mathrm{ℐ}& \mathrm{=}& \mathrm{-}\frac{\mathrm{1}}{\mathrm{\left|}\mathit{\tau }{\mathit{\omega }}_{\mathit{k}}\mathrm{\right|}}\mathrm{-}\frac{\mathrm{1}}{\mathrm{|}\mathit{\tau }{\mathit{\omega }}_{\mathit{k}}{\mathrm{|}}^{\mathit{\alpha }}}{\left(\frac{{\mathit{\tau }}_{\mathrm{e}}}{\mathit{\tau }}\right)}^{\mathrm{1}\mathrm{-}\mathit{\alpha }}\sin \left(\frac{\mathit{\alpha \pi }}{\mathrm{2}}\right)\hspace{0.17em}\mathrm{\Gamma }\mathrm{(}\mathit{\alpha }\mathrm{+}\mathrm{1}\mathrm{)}\mathit{,}\end{array}$(64)where α is an empirical adjustable parameter whose value depends on the material.

Fig. 3 Comparison between Maxwell and Andrade secular tidal torque components Tk. Parameters of Andrade rheology are α = 0.2 and τe/τ = 0.73.

As shown in Fig. 3, both models have qualitatively similar asymptotic behaviors. Thus, as explained in Makarov & Efroimsky (2013), several spin-orbit resonant configurations are at equilibrium T = 0.

Fig. 4 Normalized torque ${\sum }_{\mathit{k}}\hspace{0.17em}\mathrm{(}{\mathit{X}}_{\mathit{k}}^{\mathrm{-}\mathrm{3}\mathit{,}\mathrm{2}}{\mathrm{)}}^{\mathrm{2}}{\mathit{T}}_{\mathit{k}}^{\mathrm{\left(}\mathit{M}\mathrm{\right)}}$ for different values of τn and computed at a fixed eccentricity e = 0.4. The blue dashed line gives the tidal torque corresponding to the linear model (Eq. (57)), whose equilibrium corresponds to the pseudo-synchronization (Eq. (60)). Hansen coefficients have been obtained by FFT as described in Appendix C.

In Fig. 4, we show the evolution of the secular tidal torque as τn increases from 1 to 1000 for a given eccentricity e = 0.4. As long as τn ≲ 1, there is only one equilibrium rotation state T = 0 which corresponds to the pseudo-synchronous state. However, as τn increases beyond ~10, around each spin-orbit resonance the torque is dominated by a single component ${\mathit{T}}_{\mathit{k}}^{\mathrm{\left(}\mathit{M}\mathrm{\right)}}$ with amplitude $\mathrm{(}{\mathit{X}}_{\mathit{k}}^{\mathrm{-}\mathrm{3}\mathit{,}\mathrm{2}}{\mathrm{)}}^{\mathrm{2}}\mathit{/}\mathrm{2}$. This leads to many equilibria, but only rotations simultaneously satisfying T = 0 and dT/ dΩ < 0 are stable. As a result, only equilibria close to spin-orbit resonances are stable (Makarov & Efroimsky 2013; Storch & Lai 2014). Moreover, not all spin-orbit configurations are equilibria, since the torque also has to cross the x-axis. This condition is satisfied around a p = k₀/ 2 spin-orbit resonance if $\begin{array}{ccc}{\left[{\mathit{X}}_{\mathrm{2}\mathit{p}}^{\mathrm{-}\mathrm{3}\mathit{,}\mathrm{2}}\mathrm{\left(}\mathit{e}\mathrm{\right)}\right]}^{\mathrm{2}}\mathrm{\ge }\mathrm{2}\left|\sum _{\mathit{k}\mathrm{=}\mathrm{-}\mathrm{\infty }}^{\mathrm{\infty }}{\mathrm{[}{\mathit{X}}_{\mathit{k}}^{\mathrm{-}\mathrm{3}\mathit{,}\mathrm{2}}\mathrm{(}\mathit{e}{\mathrm{)}}^{\mathrm{]}}}^{\mathrm{2}}\frac{\mathit{\tau }{\mathit{\omega }}_{\mathit{k}}}{\mathrm{1}\mathrm{+}{\mathit{\tau }}^{\mathrm{2}}{\mathit{\omega }}_{\mathit{k}}^{\mathrm{2}}}\right|\mathrm{·}& & \end{array}$(65)In particular, for the spin-orbit resonance p = 3/2 (the last before the synchronization), by truncating the series for terms equal or higher than e², we approximately get stability as long as $\begin{array}{ccc}\mathit{e}\gtrsim \frac{\mathrm{2}}{\mathrm{7}}\sqrt{\frac{\mathrm{2}}{\mathit{\tau n}}}\mathrm{·}& & \end{array}$(66)This last expression can also be used to quickly determine which τ values allow non-synchronous rotation for a given eccentricity. In Fig. 5 we plot the critical eccentricities for some spin-orbit resonances as a function of the product τn given by expressions (65) and (66). We conclude that as the relaxation time τ increases, the spin-orbit equilibria are only broken for lower values of the eccentricity, i.e., spin-orbit resonances become more stable.

Fig. 5 Critical eccentricities for keeping the rotation in a spin-orbit resonance as a function of τn (Eq. (65)). Each color corresponds to a different spin-orbit resonance. The dashed line shows the approximation given by expression (66). Spin-orbit resonances become more stable for higher values of τ.

6. Libration

Most works on tidal dissipation usually assume the zero order solution in θ, the one in the form (Sect. 4). However, although the tidal drift on the rotation rate can be neglected over one orbital period, near spin-orbit resonances the rotation may present some libration around the equilibrium configuration. Therefore, we now compute the effect of a low amplitude libration θ₁ such that $\begin{array}{ccc}\mathit{\theta }\mathrm{=}{\mathit{\theta }}_{\mathrm{0}}\mathrm{+}{\mathit{\theta }}_{\mathrm{1}}& & \end{array}$(67)with , and $\begin{array}{ccc}{\mathit{\theta }}_{\mathrm{1}}\mathrm{=}\mathrm{\Re }\mathrm{(}\mathit{\rho }{\mathrm{e}}^{\mathrm{i}\mathit{\nu t}}\mathrm{)}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}\mathrm{(}\mathit{\rho }{\mathrm{e}}^{\mathrm{i}\mathit{\nu t}}\mathrm{+}{\mathit{\rho }}^{\mathit{\star }}{\mathrm{e}}^{\mathrm{-}\mathrm{i}\mathit{\nu t}}\mathrm{)}\mathit{.}& & \end{array}$(68)By assumption, the libration amplitude is small: | ρ | ≪ 2π; and the frequency ν is secular. We suppose also that there is an integer k₀, for which there is a resonant relation $\begin{array}{ccc}\mathrm{\left|}{\mathit{\omega }}_{{\mathit{k}}_{\mathrm{0}}}\mathrm{\right|}\mathrm{=}\mathrm{|}\mathrm{2}\mathrm{\Omega }\mathrm{-}{\mathit{k}}_{\mathrm{0}}\mathit{n}\mathrm{|}\mathrm{\ll }\mathit{\nu ,}& & \end{array}$(69)while for all k ≠ k₀, | ω_k | ≫ ν. In this case, the expression of Z^e (Eq. (33)) developped up to the first order in | ρ | reads $\begin{array}{ccc}{\mathit{Z}}^{\mathrm{e}}\mathrm{=}\sum _{\mathit{k}\mathrm{=}\mathrm{-}\mathrm{\infty }}^{\mathrm{+}\mathrm{\infty }}{\mathit{\beta }}_{\mathit{k}}\hspace{0.17em}\mathrm{(}\mathrm{1}\mathrm{+}\mathrm{2}\mathrm{i}{\mathit{\theta }}_{\mathrm{1}}\mathrm{)}{\mathrm{e}}^{\mathrm{i}{\mathit{\omega }}_{\mathit{k}}\mathit{t}}\mathrm{·}& & \end{array}$(70)Thus, from Sect. 3, we obtain $\begin{array}{ccc}\mathit{Z}\mathrm{=}\sum _{\mathit{k}\mathrm{=}\mathrm{-}\mathrm{\infty }}^{\mathrm{+}\mathrm{\infty }}{\mathit{\beta }}_{\mathit{k}}\hspace{0.17em}(\frac{\mathrm{1}\mathrm{+}\mathrm{i}{\mathit{\tau }}_{\mathrm{e}}{\mathit{\omega }}_{\mathit{k}}}{\mathrm{1}\mathrm{+}\mathrm{i}\mathit{\tau }{\mathit{\omega }}_{\mathit{k}}}\mathrm{+}\mathrm{i}\frac{\mathrm{1}\mathrm{+}\mathrm{i}{\mathit{\tau }}_{\mathrm{e}}\mathrm{(}{\mathit{\omega }}_{\mathit{k}}\mathrm{+}\mathit{\nu }\mathrm{)}}{\mathrm{1}\mathrm{+}\mathrm{i}\mathit{\tau }\mathrm{(}{\mathit{\omega }}_{\mathit{k}}\mathrm{+}\mathit{\nu }\mathrm{)}}\mathit{\rho }{\mathrm{e}}^{\mathrm{i}\mathit{\nu t}}& & \mathrm{+}\mathrm{i}\frac{\mathrm{1}\mathrm{+}\mathrm{i}{\mathit{\tau }}_{\mathrm{e}}\mathrm{(}{\mathit{\omega }}_{\mathit{k}}\mathrm{-}\mathit{\nu }\mathrm{)}}{\mathrm{1}\mathrm{+}\mathrm{i}\mathit{\tau }\mathrm{(}{\mathit{\omega }}_{\mathit{k}}\mathrm{-}\mathit{\nu }\mathrm{)}}{\mathit{\rho }}^{\mathit{\star }}{\mathrm{e}}^{\mathrm{-}\mathrm{i}\mathit{\nu t}}){\mathrm{e}}^{\mathrm{i}{\mathit{\omega }}_{\mathit{k}}\mathit{t}}\mathit{,}\\ & & \end{array}$(71)and $\begin{array}{ccc}\mathrm{\Delta }\mathit{Z}\mathrm{=}\mathrm{-}\left(\mathrm{1}\mathrm{-}\frac{{\mathit{\tau }}_{\mathrm{e}}}{\mathit{\tau }}\right)\sum _{\mathit{k}\mathrm{=}\mathrm{-}\mathrm{\infty }}^{\mathrm{+}\mathrm{\infty }}{\mathit{\beta }}_{\mathit{k}}\hspace{0.17em}(\frac{\mathrm{i}\mathit{\tau }{\mathit{\omega }}_{\mathit{k}}}{\mathrm{1}\mathrm{+}\mathrm{i}\mathit{\tau }{\mathit{\omega }}_{\mathit{k}}}\mathrm{-}\frac{\mathit{\tau }\mathrm{(}{\mathit{\omega }}_{\mathit{k}}\mathrm{+}\mathit{\nu }\mathrm{)}}{\mathrm{1}\mathrm{+}\mathrm{i}\mathit{\tau }\mathrm{(}{\mathit{\omega }}_{\mathit{k}}\mathrm{+}\mathit{\nu }\mathrm{)}}\mathit{\rho }{\mathrm{e}}^{\mathrm{i}\mathit{\nu t}}& \mathrm{-}\frac{\mathit{\tau }\mathrm{(}{\mathit{\omega }}_{\mathit{k}}\mathrm{-}\mathit{\nu }\mathrm{)}}{\mathrm{1}\mathrm{+}\mathrm{i}\mathit{\tau }\mathrm{(}{\mathit{\omega }}_{\mathit{k}}\mathrm{-}\mathit{\nu }\mathrm{)}}{\mathit{\rho }}^{\mathit{\star }}{\mathrm{e}}^{\mathrm{-}\mathrm{i}\mathit{\nu t}}){\mathrm{e}}^{\mathrm{i}{\mathit{\omega }}_{\mathit{k}}\mathit{t}}\mathit{.}& \\ & & \end{array}$(72)Since (ω_k − ω_l)t = (k − l)M, from expression (38) for $\begin{array}{c}{\mathrm{¨}}_{}\\ \mathit{\theta }\end{array}$, we obtain at first order in ρ$\begin{array}{ccc}\begin{array}{c}\mathrm{¨}\\ \mathit{\theta }\end{array}\mathrm{=}\mathrm{-}\mathrm{𝒦}\sum _{\mathit{k}\mathrm{=}\mathrm{-}\mathrm{\infty }}^{\mathrm{+}\mathrm{\infty }}\sum _{\mathit{l}\mathrm{=}\mathrm{-}\mathrm{\infty }}^{\mathrm{+}\mathrm{\infty }}{\mathit{X}}_{\mathit{k}}^{\mathrm{-}\mathrm{3}\mathit{,}\mathrm{2}}{\mathit{X}}_{\mathit{l}}^{\mathrm{-}\mathrm{3}\mathit{,}\mathrm{2}}\mathrm{\Im }\mathrm{(}{\mathit{z}}_{\mathit{k}}{\mathrm{e}}^{\mathrm{i}\mathrm{(}\mathit{k}\mathrm{-}\mathit{l}\mathrm{)}\mathit{M}}\mathrm{)}\mathit{,}& & \end{array}$(73)with $\begin{array}{ccc}{\mathit{z}}_{\mathit{k}}\mathrm{=}\frac{\mathrm{i}\mathit{\tau }{\mathit{\omega }}_{\mathit{k}}}{\mathrm{1}\mathrm{+}\mathrm{i}\mathit{\tau }{\mathit{\omega }}_{\mathit{k}}}\mathrm{+}\left(\frac{\mathit{\tau }{\mathit{\omega }}_{\mathit{k}}}{\mathrm{1}\mathrm{+}\mathrm{i}\mathit{\tau }{\mathit{\omega }}_{\mathit{k}}}\mathrm{-}\frac{\mathit{\tau }\mathrm{(}{\mathit{\omega }}_{\mathit{k}}\mathrm{+}\mathit{\nu }\mathrm{)}}{\mathrm{1}\mathrm{+}\mathrm{i}\mathit{\tau }\mathrm{(}{\mathit{\omega }}_{\mathit{k}}\mathrm{+}\mathit{\nu }\mathrm{)}}\right)\mathit{\rho }{\mathrm{e}}^{\mathrm{i}\mathit{\nu t}}& & \mathrm{+}\left(\frac{\mathit{\tau }{\mathit{\omega }}_{\mathit{k}}}{\mathrm{1}\mathrm{+}\mathrm{i}\mathit{\tau }{\mathit{\omega }}_{\mathit{k}}}\mathrm{-}\frac{\mathit{\tau }\mathrm{(}{\mathit{\omega }}_{\mathit{k}}\mathrm{-}\mathit{\nu }\mathrm{)}}{\mathrm{1}\mathrm{+}\mathrm{i}\mathit{\tau }\mathrm{(}{\mathit{\omega }}_{\mathit{k}}\mathrm{-}\mathit{\nu }\mathrm{)}}\right){\mathit{\rho }}^{\mathit{\star }}{\mathrm{e}}^{\mathrm{-}\mathrm{i}\mathit{\nu t}}\mathit{.}\\ & & \end{array}$(74)We now average over the mean anomaly M, but keeping the secular libration νt unaveraged, that is $\begin{array}{ccc}\mathrm{⟨}\begin{array}{c}\mathrm{¨}\\ \mathit{\theta }\end{array}{\mathrm{⟩}}_{\mathit{M}}\mathrm{=}\mathrm{-}\mathrm{𝒦}\sum _{\mathit{k}\mathrm{=}\mathrm{-}\mathrm{\infty }}^{\mathrm{+}\mathrm{\infty }}{\mathrm{(}{{\mathit{X}}_{\mathit{k}}^{\mathrm{-}\mathrm{3}\mathit{,}\mathrm{2}}}^{\mathrm{)}}}^{\mathrm{2}}\mathrm{\Im }\mathrm{\left(}{\mathit{z}}_{\mathit{k}}\mathrm{\right)}\mathit{.}& & \end{array}$(75)

6.1. Low libration frequency regime (τν≪ 1)

When τν ≪ 1, the expression (74) can be approximated by $\begin{array}{ccc}{\mathit{z}}_{\mathit{k}}\mathrm{\approx }\frac{\mathrm{i}\mathit{\tau }{\mathit{\omega }}_{\mathit{k}}}{\mathrm{1}\mathrm{+}\mathrm{i}\mathit{\tau }{\mathit{\omega }}_{\mathit{k}}}& & \end{array}$(76)which is exactly the expression without libration. The torque and the evolution of the system are thus identical to those described in Sect. 4.

6.2. High libration frequency regime (τν≫ 1)

Within this regime, it is necessary to distinguish the resonant term from the others. For k ≠ k₀, we have | ω_k | ≫ ν, and thus expression (74) is again $\begin{array}{ccc}{\mathit{z}}_{\mathit{k}}\mathrm{\approx }\frac{\mathrm{i}\mathit{\tau }{\mathit{\omega }}_{\mathit{k}}}{\mathrm{1}\mathrm{+}\mathrm{i}\mathit{\tau }{\mathit{\omega }}_{\mathit{k}}}\mathrm{·}& & \end{array}$(77)Conversely, for k = k₀, we have | ω_k0 | ≪ ν, and expression (74) becomes $\begin{array}{ccc}{\mathit{z}}_{{\mathit{k}}_{\mathrm{0}}}\mathrm{\approx }\frac{\mathrm{i}\mathit{\tau }{\mathit{\omega }}_{{\mathit{k}}_{\mathrm{0}}}}{\mathrm{1}\mathrm{+}\mathrm{i}\mathit{\tau }{\mathit{\omega }}_{{\mathit{k}}_{\mathrm{0}}}}\mathrm{+}\frac{\mathrm{2}\mathrm{i}{\mathit{\theta }}_{\mathrm{1}}}{\mathrm{1}\mathrm{+}\mathrm{i}\mathit{\tau }{\mathit{\omega }}_{{\mathit{k}}_{\mathrm{0}}}}\mathrm{·}& & \end{array}$(78)As a result, combining expressions (75), (77), and (78), we get $\begin{array}{ccc}\mathrm{⟨}\begin{array}{c}\mathrm{¨}\\ \mathit{\theta }\end{array}{\mathrm{⟩}}_{\mathit{M}}\mathrm{=}\mathit{T}\mathrm{-}{\mathit{\nu }}^{\mathrm{2}}{\mathit{\theta }}_{\mathrm{1}}\mathit{,}& & \end{array}$(79)where T is the non-resonant torque (Eq. (41)), and $\begin{array}{ccc}{\mathit{\nu }}^{\mathrm{2}}\mathrm{=}\mathrm{2}\mathrm{𝒦}\hspace{0.17em}{\mathrm{(}{{\mathit{X}}_{{\mathit{k}}_{\mathrm{0}}}^{\mathrm{-}\mathrm{3}\mathit{,}\mathrm{2}}}^{\mathrm{)}}}^{\mathrm{2}}\frac{\mathrm{1}}{\mathrm{1}\mathrm{+}{\mathit{\tau }}^{\mathrm{2}}{\mathit{\omega }}_{{\mathit{k}}_{\mathrm{0}}}^{\mathrm{2}}}\mathrm{·}& & \end{array}$(80)As in the previous section, the non-resonant torque is responsible for the slow variation of Ω, i.e., , and it remains $\begin{array}{ccc}\mathrm{⟨}\begin{array}{c}\mathrm{¨}\\ {\mathit{\theta }}_{\mathrm{1}}\end{array}{\mathrm{⟩}}_{\mathit{M}}\mathrm{+}{\mathit{\nu }}^{\mathrm{2}}{\mathit{\theta }}_{\mathrm{1}}\mathrm{=}\mathrm{0}& & \end{array}$(81)whose solution corresponds to our hypothesis (Eq. (68)).

7. Application to close-in planets

In this section we apply the model described in Sect. 2 to two distinct situations of exoplanets: hot Jupiters and hot super-Earths (Table 1). For that purpose, we numerically integrate the set of Eqs. (15)−(25) for the deformation of the planet, that simultaneously account for orbital and spin evolution due to the rotational flattening and to tides raised by the parent star.

In Table 1 we list the observational parameters for the planets considered in the present study. All these planets were detected through the transiting method combined with radial velocity measurements. Therefore, all orbital parameters are known with good precision, as well as the planetary masses and radius. Thus, the only unknown parameters in the model are the Love numbers (k_e and k_f) and the relaxation times (τ_e and τ).

The fluid Love numbers depend on the internal differentiation of the bodies¹. In the Solar System we have k_f ≈ 0.9 for the terrestrial planets, and 0.3 <k_f< 0.5 for the gaseous planets (Yoder 1995). The elastic Love numbers were measured with accuracy only for the Earth and Mars, respectively, k_e ≈ 0.3 (Mathews et al. 1995), and k_e ≈ 0.15 (Konopliv et al. 2006). For the gaseous planets only the ratio k₂/Q is known for Jupiter and Saturn, respectively, k₂/Q ≈ 1.1 × 10^-5 (Lainey et al. 2009), and k₂/Q ≈ 2.3 × 10^-4 (Lainey et al. 2012). The relaxation times are totally unknown for both terrestrial and gaseous planets, but if one is able to estimate Q, then τ and τ_e can be estimated through expressions (46) and (56).

7.1. Hot Jupiters

The elastic response of giant gaseous planets is unknown. Since most of the mass is in their huge atmospheres, for simplicity we assume that only the fluid deformation is important, that is, we adopt a purely viscous rheology. We note, however, that the core of these planets also experiences tidal effects, that in some cases can be as strong as those present in the convective envelope (Remus et al. 2012b; Guenel et al. 2014). In addition, other tidal mechanisms, like the excitation of inertial waves, are expected to take place in fluid planets, that also enhance the tidal dissipation (e.g., Ogilvie & Lin 2004; Favier et al. 2014).

In a viscous rheology model, we have τ_e = 0 and k_e = 0. This model is equivalent to the one proposed by Ferraz-Mello (2013), who adopts Q ~ 10⁵ and η ~ 10¹¹ Pa s, and then obtains from expression (7) that τ ~ 1 s. This timescale is intimately connected with the viscous time-lag Δτ = τ_v between the perturbation and high-tide observed (τ_v ~ δ/ω, Eq. (44)). Nevertheless, Zahn (1977) has shown that for stars possessing a convective envelope, most of the energy flux is transported by turbulent convection, leading to a strong delay of the equilibrium tide, a “convective viscous time” τ_c ~ 1 yr, that can be related with the viscous time through (Zahn 1977, Eq. (4.2)) $\begin{array}{ccc}{\mathit{\tau }}_{\mathrm{c}}\mathrm{\approx }\frac{{\mathit{R}}^{\mathrm{3}}}{\mathit{Gm}}\frac{\mathrm{1}}{{\mathit{\tau }}_{\mathit{v}}}\mathrm{·}& & \end{array}$(82)The “convective time” can also be considered for gaseous planets (that is, adopting τ = τ_c), since they can admit regions of turbulent convection (e.g., Guillot 1999; Ogilvie & Lin 2004). For planets, the effective viscosity of turbulent convection falls rapidly as the oscillation frequency is increased, and in the limit of high tidal frequencies, there is a viscoelastic response to deformation (Ogilvie & Lesur 2012). However, for low tidal frequencies (for instance, a spin-orbit commensurability) some different results should be expected. Therefore, the full deformation of gaseous planets may be on the order of a few years or even decades (for a review see Socrates et al. 2012).

Hot Jupiters are close-in Jupiter-mass planets with orbital periods less than 10 days. The in-situ formation of these planets is unlikely, because of the insufficient disk mass close to the star. The Solar System formation theories suggest that giant planets preferentially form close to the ice-line where water is condensed into ice and the necessary building blocks for the formation of planets could be found in large quantities (e.g., Ida & Lin 2008). Their presence at a close-in orbit can be explained by considering core-accretion simultaneously with disk driven migration (e.g., Lin et al. 1996; Alibert et al. 2005; Mordasini et al. 2009). An alternative explanation is obtained through gravitational interactions with a distant inclined binary companion (e.g., Wu & Murray 2003; Correia et al. 2011) or planet-planet scattering (e.g. Rasio & Ford 1996; Beaugé & Nesvorný 2012) combined with tidal effects. Contrarily to the migration scenario, where the eccentricities are quickly damped to zero, the gravitational interaction mechanism can also explain why some proto-hot Jupiters are still observed with high eccentricity values (e.g., Dawson & Murray-Clay 2013). Since for circular orbits the only possibility for the spin is synchronous rotation (e.g., MacDonald 1964; Hut 1980), in this paper we focus on eccentric planets. The solution for circular orbits is also obtained at the end of the evolution for initially eccentric planets.

Here, we apply our model to the planets HAT-P-34 b (Bakos et al. 2012) and HD 80606 b (Hébrard et al. 2010), both belonging to single-planet systems. The main difference between these two planets is their present orbital period, only 5.45 day for HAT-P-34 b, and about 111 day for HD 80606 b (Table 1). However, due to the angular momentum conservation, when the eccentricity becomes fully damped, the final semi-major axis and orbital period are constrained (e.g., Correia 2009) $\begin{array}{ccc}{\mathit{a}}_{\mathrm{f}}\mathrm{=}\mathit{a}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{e}}^{\mathrm{2}}\mathrm{)} \mathrm{\iff } {\mathit{P}}_{\mathrm{f}}\mathrm{=}\mathit{P}{\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{e}}^{\mathrm{2}}\mathrm{)}}^{\mathrm{3}\mathit{/}\mathrm{2}}\mathit{.}& & \end{array}$(83)This gives 3.94 day and 5.19 day for HAT-P-34 b and HD 80606 b, respectively, that is, both planets will evolve to “regular” hot Jupiter configurations. For HAT-P-34 b we begin our simulations with P = 10.82 day and e = 0.7, while for HD 80606 b we start with the present observed values (Table 1). The initial rotation period is 0.5 day for both planets, and the planet is assumed perfectly spherical (J₂ = C₂₂ = S₂₂ = 0). Because the relaxation time is unknown, we adopt different values for τ ranging from 10^-5−10⁰ yr, in order to cover all kinds of behavior. For the fluid Love number we use Jupiter’s value k_f = 0.5 (Yoder 1995).

In Fig. 6 we show the future evolution of the HD 80606 b orbit (top), rotation (middle) and shape (bottom) using different values of the relaxation time. The shape of the planet is described by its gravity field coefficients. However, since C₂₂ and S₂₂ values constantly change for a non-synchronous planet, we plot the prolateness of the planet, $\begin{array}{ccc}\mathit{ϵ}\mathrm{=}\mathrm{\left|}\mathit{Z}\mathrm{\right|}\mathrm{=}\sqrt{{\mathit{C}}_{\mathrm{22}}^{\mathrm{2}}\mathrm{+}{\mathit{S}}_{\mathrm{22}}^{\mathrm{2}}}\mathit{,}& & \end{array}$(84)which is equivalent to the value of C₂₂ along the instantaneous principal long axis of inertia. In all situations, the planet attempts to point the instantaneous long axis to the star, increasing its prolateness until it reaches the pericenter of the orbit. At this stage, the prolateness decreases again, reaching a minima at the apocenter.

Fig. 6 Evolution of HD 80606 b with time for different τ values. We plot the semi-major axis (in AU) and the eccentricity (top), the ratio between the planet rotation rate and the orbital mean motion (middle), and the planet J2 and ϵ (bottom). The green line gives the pseudo-synchronous equilibrium (Eq. (60)) (middle), and the average equilibrium values for J2 and ϵ, respectively (Eqs. (85) and (86)) (bottom).

For τ = 10^-4 yr (Fig. 6, left), the relaxation time is much shorter than the orbital period, so between two consecutive passages at the apocenter, the prolateness of the planet is back to its minimal value. The only effect on the rotation rate is then to decrease it from its initial value, until it reaches the pseudo-synchronous equilibrium (Eq. (60)), since we have τω_k ≪ 1 (Sect. 5.1). Because of the conservation of the angular momentum, the eccentricity of the orbit is also progressively damped to zero. When the eccentricity is close to zero, Ω /n ≈ 1 + 6e² (Eq. (60)), so the equilibrium rotation tends to the synchronous rotation. At this stage the planet is able to acquire a permanent prolateness, always pointing the long axis to the star.

For τ = 10^-2 yr (Fig. 6, center), the relaxation time is also shorter than the orbital period during the initial stages of the spin evolution, so the rotation behaves similarly to the previous case (τ = 10^-4 yr). However, as the orbit shrinks, the relaxation time and the orbital period become equivalent τn ~ 1. Then, when the planet is back at the pericenter, its prolateness is still excited by the previous pericenter passage. As a consequence, the planet increases its global prolateness at each pericenter passage. This effect is maximized near a spin-orbit commensurability, since the rotation rate is close to a harmonic of the orbital mean motion, Ω /n = p = k₀/ 2, and thus there exist one tidal frequency ω_2p that is near zero (Sect. 5.2). The gravitational torque exerted by the star on the prolate planet dominates the tidal toque and tends to keep the rotation near the spin-orbit commensurability (Fig. 4). Therefore, in this regime the rotation rate closely follows the spin-orbit equilibria.

Finally, for τ = 1 yr (Fig. 6, right), the relaxation time is permanently larger than the orbital period. Therefore, the prolateness of the planet is continuously excited by a nearby spin-orbit resonance, as explained for the final stages in the previous case. The spin remains trapped in a spin-orbit equilibrium until the eccentricity decreases below a given threshold (Fig. 5). Then, the spin decreases one more step until it is captured in the following lower-order spin-orbit resonance. The process continues until the synchronous equilibrium is reached for small values of the eccentricity, since this is the only spin-orbit resonance that is always stable (Eq. (65)).

For lower and higher τ values the system asymptotically behaves as in the τ = 10^-5 yr and 10⁰ yr cases, respectively, the major change is in the evolution timescale, that is longer.

7.2. Evolution versus eccentricity

In Fig. 6, the evolution timescale is different for each τ value, since tidal dissipation is controled by this parameter (Fig. 2). However, because of the conservation of the angular momentum (Eq. (83)), the final semi-major axis is always the same. In addition, there is a perfect correspondence between the semi-major axis and the eccentricity (Fig. 7), that is, the orbital evolution is always the same, only the timescale changes. Therefore, in order to better compare the behavior of the spin with different τ values, instead of using the time, it is preferable to use the semi-major axis or the eccentricity in the evolution axis.

Fig. 7 Evolution of the semi-major axis as a function of the eccentricity for different systems and τ values. The dotted lines correspond to the theoretical prediction given by expression (83). Solid lines correspond to the orbital evolution obtained using numerical simulations with different τ values. For each system, the simulations exactly match the theoretical prediction.

In Fig. 8, we plot the evolution of the spin of HAT-P-34 b and HD 80606 b as a function of the eccentricity for different τ values. The initial rotation period is arbitrarily set at 0.5 day for both planets in all simulations. This value is not critical, because for all τ values the tidal torque on the spin is so strong that the rotation quickly evolves into an equilibrium configuration. We observe that the spin evolution of both planets is identical, although the initial eccentricity of HAT-P-34 is 0.7, while the initial eccentricity of HD 80606 is 0.933. We can also clearly observe the two tidal regimes described in Sect. 5. For low eccentricity values, the orbital period of both planets is around 10^-2 yr. Thus, for τ ≪ 10^-2 yr, the spin closely follows the pseudo-synchronous equilibrium (Eq. (60)), since τn ≪ 1. For τ ~ 10^-2 yr, we observe the transition between regimes, while for τ ≫ 10^-2 yr, the spin is captured in half-integer spin-orbit resonances, since τn ≫ 1. These resonances become destabilized as the eccentricity decays, in agreement with expression (65). The main consequence of these metastable equilibria is to allow faster rotation rates for an identical eccentricity.