Issue 
A&A
Volume 629, September 2019



Article Number  A132  
Number of page(s)  27  
Section  Planets and planetary systems  
DOI  https://doi.org/10.1051/00046361/201935905  
Published online  17 September 2019 
Final spin states of eccentric ocean planets
^{1}
Laboratoire d’Astrophysique de Bordeaux, Université Bordeaux, CNRS,
B18N, allée Geoffroy SaintHilaire,
33615
Pessac,
France
^{2}
University of Bern, Center for Space and Habitability,
Gesellschaftsstrasse 6,
3012
Bern,
Switzerland
email: pierre.auclairdesrotour@csh.unibe.ch
^{3}
Observatoire de Genève, Université de Genève,
51 Chemin des Maillettes,
1290
Sauvergny, Switzerland
^{4}
AIM, CEA, CNRS, Université ParisSaclay, Université Paris Diderot, Sorbonne Paris Cité,
91191
GifsurYvette Cedex,
France
Received:
16
May
2019
Accepted:
12
July
2019
Context. Eccentricity tides generate a torque that can drive an ocean planet towards asynchronous rotation states of equilibrium when enhanced by resonances associated with the oceanic tidal modes.
Aims. We investigate the impact of eccentricity tides on the rotation of rocky planets hosting a thin uniform ocean and orbiting cool dwarf stars such as TRAPPIST1, with orbital periods ~1−10 days.
Methods. Combining the linear theory of oceanic tides in the shallow water approximation with the Andrade model for the solid part of the planet, we developed a global model including the coupling effects of ocean loading, selfattraction, and deformation of the solid regions. From this model we derive analytic solutions for the tidal Love numbers and torque exerted on the planet. These solutions are used with realistic values of parameters provided by advanced models of the internal structure and tidal oscillations of solid bodies to explore the parameter space both analytically and numerically.
Results. Our model allows us to fully characterise the frequencyresonant tidal response of the planet, and particularly the features of resonances associated with the oceanic tidal modes (eigenfrequencies, resulting maxima of the tidal torque, and Love numbers) as functions of the planet parameters (mass, radius, Andrade parameters, ocean depth, and Rayleigh drag frequency). Resonances associated with the oceanic tide decrease the critical eccentricity beyond which asynchronous rotation states distinct from the usual spinorbit resonances can exist. We provide an estimation and scaling laws for this critical eccentricity, which is found to be lowered by roughly one order of magnitude, switching from ~0.3 to ~0.06 in typical cases and to ~0.01 in extremal ones.
Key words: hydrodynamics / planet / star interactions / planets and satellites: oceans / planets and satellites: terrestrial planets
© P. AuclairDesrotour et al. 2019
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
1 Introduction
One of the most essential questions raised by the discovery of rocky exoplanets is the nature of their climate and surface conditions. Particularly, this question motivated a major part of studies dealing with the TRAPPIST1 system (Gillon et al. 2017), where an ultracool Mdwarf star harbours seven Earthsized planets among which four – namely planets d, e, f, and g – may be potentially habitable (e.g. Bolmont et al. 2017; Bourrier et al. 2017; Grimm et al. 2018; Papaloizou et al. 2018; Unterborn et al. 2018; Barr et al. 2018; Turbet et al. 2018; Dobos et al. 2019). To constrain the climate of such planets, it is crucial to preliminarily constrain their rotation. As tides drive the long term evolution of planetary systems, this requires the characterisation of the possible asynchronous states of equilibrium where rocky planets may be tidally locked.
Rocky planets orbiting cool dwarf stars in tightly packed systems combine temperate surfaces conditions with small orbital periods, which makes them privileged objects of study. The TRAPPIST1 system illustrates this configuration well since its seven rocky planets orbit within a radius disk of 0.065 au, and two of them – planets d and e – are located in the habitable zone of the host star with orbital periods of 4.05 and 6.10 days, respectively (Gillon et al. 2017; Grimm et al. 2018). Under these conditions, the presence of a considerable amount of liquid water on the planet makes the existence of global oceans likely (e.g. Bolmont et al. 2017; Bourrier et al. 2017; Turbet et al. 2018).
The tidal response of a freesurface oceanic layer strongly differs from that of solid bodies. Dry rocky planets undergoing the tidal gravitational potential of their host star, or perturber, are subject to small distortions due to their elasticity (e.g. Henning et al. 2009; Efroimsky 2012; Remus et al. 2012). The resulting tidal elongation corresponds to the hydrostatic adjustment between gravity and elasticity, which is the socalled equilibrium tide (Zahn 1966; Ogilvie & Lin 2004). Because of the lag induced by internal dissipative processes, the tidal bulge generates a torque acting on the rotation of the body. In the case of a circular and coplanar planet, this torque is created solely by the semidiurnal tide and drives the planet towards spinorbit synchronous rotation. Unlike solid tides, which are not strongly dependent on the tidal frequency, oceanic tides exhibit a frequencyresonant behaviour that enables large variations of the tidal torque (e.g. Webb 1980; Tyler 2011; Chen et al. 2014; Matsuyama 2014; AuclairDesrotour et al. 2018).
In thin fluid shells, such behaviour results from the propagation of gravitoinertial surface modes (or barotropic modes) forced by the tidal gravitational potential, which corresponds to the solutions of Laplace’s tidal equation (e.g. LonguetHiggins 1968; Webb 1980). In deeper oceans, stable stratification leads to the propagation of internal gravitoinertial waves, which are restored by the Archimedean force and provide a supplementary resonant contribution to the oceanic tidal torque (e.g. Tyler 2011; AuclairDesrotour et al. 2018). These modes form the baroclinic tide. The component of the tidal response associated with the propagation of waves is named “dynamical tide” (e.g. Zahn 1975; Ogilvie & Lin 2004) in opposition with the equilibrium tide.
The tidal response of the ocean is tightly coupled with that of the solid part by selfattraction and loading, meaning that the contributions of the two layers cannot be separated but merge into an effective tidal response of the planet. Although these couplings are often ignored for the sake of simplicity (e.g. Tyler 2011; Chen et al. 2014; AuclairDesrotour et al. 2018), they significantly modify the global tidal response of the planet (e.g. Matsuyama 2014; Matsuyama et al. 2018) and should be taken into accountin the calculation of the effective tidal torque and Love numbers (Love 1911), which quantify the impact of tidal dissipation on the evolution of the planetperturber system (see e.g. Munk & MacDonald 1960, for solid bodies).
The specific action of tides on the planetary spin can be characterised in the idealised coplanar and circular configuration, that is in the absence of obliquity and eccentricity (e.g. AuclairDesrotour et al. 2018). In this configuration, the semidiurnal tidal component predominates and drives the planet towards the spinorbit synchronous rotation, except in the case of thermal atmospherictides (e.g. Lindzen & Chapman 1969), where the body is torqued away from synchronisation (e.g. Gold & Soter 1969; Ingersoll & Dobrovolskis 1978; Dobrovolskis & Ingersoll 1980; Correia & Laskar 2001, 2003; Leconte et al. 2015; AuclairDesrotour et al. 2017a,b, 2019).
Eccentricity affects this evolution by introducing supplementary forcing terms, weighted by the Hansen coefficients (e.g. Hughes 1981), which can torque the body away from synchronisation (e.g. Greenberg 2009; Makarov & Efroimsky 2013; Correia et al. 2014). A dry rocky planet with an eccentric orbit may thus be tidally locked into one of the spinorbit resonances induced by eccentricity tides (3:2, 2: 1, 5: 2, 3: 1, etc.; see Makarov 2012; Makarov & Efroimsky 2013; Correia et al. 2014), which correspond to asynchronous rotation rates. This question has been investigated mainly for the rocky bodies of the Solar system, such as Mercury, by means of twolayer tidal models composed of a molten core and a solid crust (e.g. Peale & Boss 1977; Correia & Laskar 2009; Henning & Hurford 2014; Noyelles et al. 2014). The spinorbit rotation equilibria were also retrieved in the case of giant planets by studies using a viscoelastic model of the tidally dissipated energy (e.g. Storch & Lai 2014).
As shown by many studies examining the case of icy satellites in the Solar system (e.g. Tyler 2008, 2009, 2011; Chen et al. 2014; Beuthe 2016; Matsuyama 2014; Matsuyama et al. 2018), the resonances of oceanic modes may enhance eccentricity tidal components by several orders of magnitude and increase as well the resulting tidal heating. As a consequence, they are also likely to generate asynchronous rotation states of equilibrium distinct from the spinorbit resonances identified by early works for solid bodies, including for low eccentricities. This means that the combination of eccentricity tides with resonances enables the existence of such states for rocky planets exhibiting quasicircular orbits.
In the present work, we investigate this mechanism by considering the case of an idealised ocean planet with an eccentric orbit and no obliquity, with the central object of the system assumed to be a TRAPPIST1like dwarf star. Treating the ocean as a spherical thin shell in the shallow water approximation, we follow Matsuyama (2014) and take the effects of ocean loading, selfattraction, and deformation of solid regions selfconsistently into account. We shall emphasise here that we do not include the coupling between the internal tidal heating of the planet and its structure and surface conditions, these latter conditions being fixed in the model. Therefore, the planetary system is essentially parametrised by the mass of the host star, which is set to the TRAPPIST1 value. The implications of tidal heating on the possible existence of water oceans are discussed in the conclusion.
To describe the tidal response of the solid part, we used the Andrade model (Andrade 1910; CastilloRogez et al. 2011; Efroimsky 2012) with values of parameters provided by a spectral code that computes the viscoelastic tidal oscillations of the body expanded in spherical harmonics by taking as input its density, rigidity, and viscosity profiles (Takeuchi & Saito 1972; Tobie et al. 2005). For a given mass and composition for the planet, the internal structure used for these calculations is derived from the model detailed in Sotin et al. (2007). For the derived internal structure of the planet, we then used the model of Tobie et al. (2005), which solves an equation of state selfconsistently across the radial direction.
The choice of the Andrade model is mainly motivated by the abundance of experimental data supporting it for telluric planets, the model having been reported to match over a wide range of experimental conditions (e.g. Andrade 1910, 1914; Cottrell & Aytekin 1947; Duval 1978; Jackson 1993). Particularly, the Andrade model is believed to better describe the behaviour of terrestrial bodies than the Maxwell rheology in the highfrequency range (e.g. Efroimsky & Lainey 2007; Efroimsky 2012). However, the methodology applied in this article is general and can thus be easily applied to any rheology, including more complex ones such as the SundbergCooper rheology (e.g. Renaud & Henning 2018).
The ability of eccentricity tides to drive the planet away from the spinorbit synchronous rotation is quantified by the minimal eccentricity for which asynchronous states may exist. All along the article, this eccentricity is called the “critical eccentricity” and denoted by e_{AR} (“AR” refers to asynchronous rotation). As we focus on the frequency interval bounded by the 1: 1 (synchronisation) and 3:2 spinorbit resonances, e < e_{AR} means that the planet is driven towards the spinorbit synchronous rotation, while e > e_{AR} corresponds to asynchronous final rotation states of equilibrium.
In Sect. 2 we introduce the physical setup, the main parameters, and the reference frames of the study. In Sect. 3, we expand the perturbing tidal gravitational potential in eccentricity series, and establish the equations and quantities governing the tidal response of the solid part and ocean. Particularly, we detail the features of oceanic tidal modes derived from Laplace’s tidal equation. In Sect. 4, we establish the analytic expressions of the tidal torque and Love numbers in the general case, and in the asymptotic cases corresponding to pure solid and oceanic tidal responses. This allows us to analyse the effect of resonances on eccentricity terms, and to characterise the frequency behaviour of the tidal torque. In Sect. 5, we derive analytic estimations of the critical eccentricity in the quasiadiabatic asymptotic regime.
In Sect. 6, the final rotation states of Earth and superEarthsized planets are calculated numerically as functions of the eccentricity and ocean depth for various orbital periods and Rayleigh drag timescales, which highlights that the critical eccentricity may be decreased by one order of magnitude owing to the action of resonances associated with oceanic tidal modes. The existence of these final spin states is discussed in terms of possible capture in the 1: 1 spinorbit resonance using the theory of Goldreich & Peale (1966). In Sect. 7, the critical eccentricity is calculated as a function of the orbital period and ocean depth in order to unravel the regions of the parameter space where asynchronous spin equilibria are compatible with low eccentricities. In Sect. 8, we integrate the evolution of the planet spin over time for various eccentricities, and thus quantify the evolution timescale associated with the tidal torque. Finally, we give our conclusions and discuss the limitations of the model in Sect. 9.
2 Physical setup
We examined the simplified case of a single terrestrial planet orbiting its host star with an eccentric orbit. The planet, of mass M_{p} and radius R_{p}, is treated as a bilayered body basically composed of an internal solid part and an external thin incompressible ocean of uniform thickness H_{oc} ≪ R_{p} and density ρ_{oc}. Its orbital motion is described by its mean motion n_{⋆} and eccentricity e. The orbital angular momentum vector is denoted n_{⋆}. The planetocentric referential is , where Z_{G} ≡n_{⋆}∕n_{⋆} (the symbol ≡ meaning “defined by”), X_{G} and Y_{G} designate the directions of two distant stars defining the orbital plane, and O is the planet’s gravity center.
We assumed the absence of obliquity, which allows us to define the reference frame corotating with the planet as with Z_{p} = Z_{G}. The vectors X_{p} and Y_{p} define the equatorial plane of the planet, which is also its orbital plane in the present case. The rotating motion of with respect to is defined by the spin vector Ω = ΩZ_{p}, Ω being the rotation rate of the planet. Hence, denoting t the time, and . Finally the spherical coordinates are radius r, colatitude θ, and longitude φ, and the associated unitvector basis .
The surface planet gravity is denoted as g. The planet mean density is and the mass of the ocean . In addition, the mass and radius of the solid part are denoted M_{c} ≡ M_{p} − M_{oc} ≈ M_{p} and R_{c} ≡ R_{p} − H_{oc} ≈ R_{p}, respectively.To simplify calculations, we assumed in the following that this layer can be treated as a homogeneous body of density .
The ocean is considered as perfectly coupled with the solid part of the planet by frictional forces, meaning that the whole planet rotates as a solid body. This assumption holds provided that the ocean is thin with respect to the planet radius. In the case of a deep ocean, the solid part of the planet and the upper layers of the ocean decouple, which allow zonal flows to develop and potentially give rise to differential rotation depending on the strength of friction. Owing to this decoupling, one may envision that a sufficiently deep ocean and the solid part could be rotating in different spin states, similar to the core and the mantle in rocky planets (e.g. Correia & Laskar 2009).
In the framework of the Rayleigh drag approximation chosen for this study (e.g. Vallis 2006), the typical timescale of the coreocean coupling by viscous friction is , where σ_{R} is the effective Rayleigh drag frequency parametrising the dissipative term in the momentum equation (see Eq. (18)). In the case of the Earth, with σ_{R} = 10^{−5} s^{−1} (Webb 1980), τ_{R} ≈ 30 h, which is far smaller than the evolution timescale of the planet spin.
3 Tidal dynamics
Owing to its eccentric orbit, the planet is subject to eccentricity tides, which act on its spin rotation in a different way from the standard semidiurnal tide. While the semidiurnal tide drives the body towards the spinorbit synchronous rotation, eccentricity tides tend to desynchronise it, and thereby enable the existence of nonsynchronised states of equilibrium. In this section, we establish the components of the tidal potential and discuss the desynchronising mechanism of eccentricity tides. We then use the classic tidal theory to analytically describe the dynamics of the planet tidal response by successively treating the solid part and the oceanic layer.
Fig. 1 Tidal distortion of terrestrial planet resulting from gravitational forcing of host star. The planet structure is composed of a solid part of radius R_{c} and a thin uniform ocean of external radius R_{p} rotating as a solid body at the spin rotation rate Ω. The tidal responses of the two layers are coupled together and characterised by different frequency behaviours. The corresponding angular lags are designated by δ_{c} and δ_{oc}, and the surface displacements by ξ_{c} and ξ_{oc}. The ocean depth modified by solid and oceanic tides is denoted by h. Reference frames introduced in Sect. 2 are drawn in purple. 
3.1 Perturbing tidal gravitational potential
The whole planet is tidally forced by the gravitational potential of the host star. For large star–planet distances r_{⋆}, that is r_{⋆} ≫ R_{p} typically, this gravitational potential is expressed in the accelerated frame of the planet as (1)
the first term of the righthand side of the equation corresponding to the attraction by the host star, of mass M_{⋆}, and the second term to the centrifugal force due to the orbital motion. We note that the symbol ^{^} is employed here and throughout the article to highlight real quantities with respect to complex ones, the latter being preferentially used in analytical developments in the general case. Conversely, a few specific complex quantities that are usually real will be highlighted by the symbol ^{˜} in order to avoid confusion.
In the thin layer approximation, the tidal gravitational potential is approximated by its value at the planet surface (r = R_{p}), (2)
where we have removed the constant component as it does not contribute to the tidal force. The associated complex gravitational potential U_{T}, such that , is expanded in Fourier series of the time and spherical harmonics, following Kaula’s theory (e.g. Kaula 1966). In the absence of obliquity, it is thus written (3)
where l and m designate the latitudinal and longitudinal degrees, s an integer, σ_{m,s} = mΩ − sn_{⋆} the forcing tidal frequency of the mode associated with the doublet , the normalised associated Legendre function associated with the doublet (see Appendix A), and the surface tidal gravitational potential associated with the triplet . This later isgiven by (4)
where we have introduced the semimajor axis a, and the dimensionless coefficients A_{l,m,s}, following the notation by Ogilvie (2014) (Eq. (3) of the review; we note that we use here the normalised associated Legendre functions, which explains how we get different numerical factors). The Kronecker symbol is δ_{l,k}, such that δ_{l,k} = 1 if l = k and δ_{l,k} = 0 otherwise. By analogy with the Kronecker symbol, we define the coefficient by δ_{s<0}, such that δ_{s<0} = 1 if s < 0, and 0 otherwise. The A_{l,m,s} coefficients are thus expressed as (5)
In the aboveexpression, the are the unnormalised associated Legendre functions, and the eccentricity functions are the socalledHansen coefficients (Hughes 1981; Polfliet & Smeyers 1990; Laskar 2005), which are calculated numerically in the study using a fast Fourier transform (see Appendix B for details).
Owing to the radial behaviour , terms of degrees l > 2 are negligible with respect to second order components when the starplanet distance far exceeds the radius of the planet (Mathis & Le PoncinLafitte 2009). They can thus be ignored in the framework of this study. Quadrupolar terms (l = m = 2) are associated with the eccentricity frequencies (6)
with . In the case of a circular orbit (e = 0), the tidal potential reduces to the semidiurnal component (8)
At smalleccentricities (e ≪ 1), eccentricity tides are predominantly forced by terms associated with s = 1 (the westward propagating potential) and s = 3 (the eastward propagating potential), given by (e.g. Ogilvie 2014, Table 1)
which are associated with the tidal frequencies σ_{2,1} = 2Ω − n_{⋆} and σ_{2,3} = 2Ω − 3n_{⋆}, respectively. At high eccentricities (e ≲ 1), the order s of the predominant component increases with e, while the spectrum of forcing terms widens (see Ogilvie 2014, Fig. 3). Besides, we note that terms associated with negative s are always negligible compared to those associated with positive s.
The transition between low and higheccentricity regimes occurs for , that is at the transition eccentricity e_{trans} ≈ 0.243 in the second order approximation in e (see Table C.1 in Correia et al. 2014), independently from the system parameters. In the absence of resonances, the semidiurnal tide remains predominant as long as e ≲ e_{trans}, meaning that the critical eccentricity beyond which asynchronous final rotation states of equilibrium can exist is e_{AR} = e_{trans} in this case. As shown in the following, e_{AR} can be strongly decreased by the resonances associated with the oceanic tidal response.
Figure 2 illustrates the desynchronising mechanism of eccentricity tides, through the action on rotation of the components of the tidal gravitational potential given by Eqs. (8)–(10). The semidiurnal component, given by Eq. (8), is locked on the average direction of the satellite. The two other components, given by Eqs. (9) and (10) and induced by eccentricity, travel clockwise and counterclockwise, respectively (see Greenberg 2009, Fig. 3). Dissipative processes, such as viscous friction, induce an angular lag between the tidal bulge and the direction of the star.
The configuration shown by Fig. 2 is supersynchronous and corresponds to . In this case, the planet is torqued towards the spinorbit synchronous rotation by the semidiurnal component and the eccentricity component of degree s = 1, which results from the westward displacement of the star in the reference frame of the planet in the vicinity of the apoapsis. Conversely, its rotation is accelerated by the component of degree s = 3 resulting from the westward displacement of the star in the vicinity of the periastron. In the figure, this component is assumed to be predominant. Hence the spin rotation of the planet is accelerated in the vicinity of the periastron.
In the case of subsynchronous rotation (Ω ≤ n_{⋆}), the eccentricity component of degree s = 3 torques the planet towards synchronisation, as is the case for the semidiurnal component. The desynchronising component is that of degree s = 1 in this case, from the moment that . The role played by the two predominating eccentricity components in the approximation of small eccentricities may be generalised to other ones (). The degrees eccentricity component tends to drive the planet away from synchronisation provided that Ω satisfies the conditions (11)
Otherwise, the component pushes the planet towards synchronisation. The global effect resulting from the combinations of tidal components depends on the behaviour of the solid part and the oceanic layer coupled together. Particularly, it is sensitive to resonances proper to the oceanic tidal response (Webb 1980; AuclairDesrotour et al. 2018). These resonances induce a strong dependence of the torque on the tidal frequency, which is not the case of the behaviour of the solid part (e.g. Efroimsky & Lainey 2007). We thus characterise the global tidal response of the planet in the following section.
Fig. 2 Desynchronising mechanism of eccentricity tides. Tidal bulge raised on a planet by the host star in the heliocentric reference frame (left) and in a planetocentric reference frame rotating with the planet mean motion (right). The tidal gravitational forcing is assumed to be dominated by its eccentricity components, given by Eqs. (9) and (10). Because of the angular lag generated by dissipative processes, the rotation of the planet is accelerated by the component of degree s = 3 (the eastward propagating potential) in the vicinity of the periastron. Conversely, it is decelerated by the component of degree s = 1 (the westward propagating potential) to a lesser extent in the vicinity of the apoapsis. 
3.2 Tidal response of the solid part
As a first step, we focus on the tidal response of the solid part of the planet. In the general case, the distortion associated with the gravitational tidal forcing of the host star can be considered a small perturbation, meaning that variations of physical quantities characterising the planet are small with respect to the hydrostatic background. This is the framework where the linear tidal theory used in this work can be applied.
Computing the response of the solid part requires the preliminary assumption of a rheology for the material that composes it. Hence, following Efroimsky (2012), we consider that the solid part behaves as an isotropic viscoelastic body of unrelaxed shear modulus μ. This body undergoes a timeperiodic tidal force offrequency σ. Thus, in the linear approximation, any quantity describing the distortion can be written as , where q^{σ} stands for the associated complex quantity, which depends on the forcing frequency. We introduce the stress tensor and the strain tensor in the linear approximation (12)
the notation ξ designating the displacement vector. The rheology of the solid is determined by the relationship between and .
To model this relationship, we make use of a generalised version of Hooke’s law, which is a linear constitutive law that governs the distortion of isotropic elastic materials as long as this distortion does not exceed the elastic limit of the material. Neglecting compressibility, we only take the deviatoric stresses and strain into consideration. Therefore, Hooke’s law reduces in the present case to (13)
where is the complex, shear modulus accounting for the material rheology. We note that the imaginary part of characterises the anelasticity of the material due to the action of internal dissipative processes, such as friction. In the case of a purely elastic body, would be real.
In this study, we opt for an Andrade rheology (see Andrade 1910; CastilloRogez et al. 2011; Efroimsky 2012, for details), which allows us to write as an explicit function of the forcing frequency, (14)
where we have introduced the gamma function Γ (Abramowitz & Stegun 1972); the dimensionless rheological exponent α (the values of α are determined experimentally and typically fallwithin the interval 0.2–0.4 for olivinerich rocks, see Castelnau et al. 2008); the Maxwell relaxation time τ_{M} associated with the viscous response of the material (τ_{M} is defined by the ratio of viscosity η to unrelaxedrigidity, τ_{M} ≡ η∕μ, see e.g. CastilloRogez et al. 2011); and the Andrade – or anelastic – time τ_{A} related to anelasticity. This last parameter is the characteristic timescale associated with the material creep.
By setting τ_{A} = +∞, one retrieves the standard Maxwell viscoelastic rheology (Greenberg 2009; Efroimsky 2012; Correia et al. 2014). In this case, the behaviour of the solid part follows two asymptotic regimes defined by the hierarchy between the tidal period P_{T} and the Maxwell time. In the zerofrequency limit (P_{T} ≫ τ_{M}), the response is dominated by viscous friction, leading to the purely imaginary complex shear modulus (again σ designates the forcing frequency of the tidal mode). In the highfrequency regime (P_{T} ≪ τ_{M}), the behaviour of the solid part tends to be purely elastic, with . In the case of telluric planets, the order of magnitude of the Maxwell time is of several centuries (typically τ_{M} ~ 10^{4}–10^{5} days; for the Earth τ_{M} is about 500yr, see e.g. Efroimsky 2012), which is far greater than usual tidal periods (P_{T} ~ 10^{0}−10^{2} days). As a consequence, the imaginary part of drops rapidly as σ increases in the highfrequency regime, which leads to underestimating the energy that is tidally dissipated in the planet interior (this effect has repercussions on the frequencydependence of the imaginary part of Love numbers, as shown by Fig. 3 in the following).
The main interest of the Andrade model is to improve this behaviour by providing more realistic orders of magnitude of the energy that is tidally dissipated within planetary interiors with a minimal number of additional parameters, namely α and τ_{A}. In the Andrade model, the decay of the anelastic component with the forcing frequency is attenuated with respect to that described by the Maxwell model, the slope of the decay being determined by the rheological parameter α (this is illustratedby Fig. 3, which will be discussed below, where the imaginary part of the quadrupole Love number is plotted as a function of the tidal frequency). This description of the tidal behaviour may be refined using more sophisticated rheological models like the SundbergCooper rheology, which implies a larger tidal dissipation around a critical temperature and frequency (Renaud & Henning 2018).
As discussed by early studies (CastilloRogez et al. 2011; Efroimsky 2012), estimating the Andrade time is far from being trivial. Particularly, there is no reason for τ_{A} and τ_{M} to be comparable in the general case since they are related to two different physical mechanisms. The order of magnitude of τ_{A} can be derived from microphysics through the study of the propagation of seismic waves, as done for instance in the case of the Earth mantle (e.g. Karato & Spetzler 1990; Tan et al. 2001; Jackson et al. 2002). For more details, we refer the reader to CastilloRogez et al. (2011), where different sets of data are compared in order to determine the Andrade parameter .
For the purpose of this work, we will use in calculations the effective values of parameters derived from the numerical integration of the equations of the elastogravitational theory (Love 1911; Takeuchi & Saito 1972) with realistic radial profiles of background quantities (see Table 1). The tidal model used to obtain these values is described in Tobie et al. (2005), and treats the deformation of a spherically symmetric, nonrotating, and elastic body, for which the momentum, Poisson, and mass conservation equations are solved by means of the classical propagator matrix method (e.g. Sabadini & Vermeersen 2004) over a sampling of the tidal frequency.
The radial profiles of background quantities used as inputs of the tidal model are computed with a code (as detailed in Sotin et al. 2007)that numerically solves the hydrostatic balance across the radial direction assuming an equation of state and a composition for the planet. In these calculations of the internal structure, the planet is considered to be a dry rocky body owing to the negligible impact of the thin oceanic shell on the parameters of the solid viscoelastic tidal response. In the case of a thick oceanic layer, these parameters should be derived by taking the ocean into account as the latter can no longer be neglected (Dermott 1979; Remus et al. 2012).
Values of Andrade parameters are obtained by fitting the Andrade model to the frequencyspectra derived from the calculations of tidal oscillations using the method of Tobie et al. (2005). These values may be regarded as the effective parameters of an equivalent body that is homogeneous in density and presents the same frequency behaviour as the studied body of realistic internal structure. An application of this method to the TRAPPIST1 system may be found in Breton et al. (2018), and will be detailed in a forthcoming article (Bolmont et al. in prep.).
The variation of mass distribution associated with the tidal distortion is described by Love numbers. In the linear approach, each degreel mode of theexpansion in spherical harmonics (see Eq. (3)) has its associated Love numbers, denoted by , , , and in the present work. The first two parameters, and , stand for the tidal gravitational and displacement Love numbers, respectively. They characterise the variations of the selfgravitational potential and vertical displacement induced by the internal gravitational tidal forcing at the surface of the solid body.
The notations and designate the load Love numbers for the gravitational potential and vertical displacement of the surface, respectively. These numbers characterise the response of the planet to the gravitational force and surface pressure induced by the variation of mass distribution of a thin external fluid layer. In the case of a uniform interior, the degreel Love numbers (e.g. Munk & MacDonald 1960) are given by (15)
where we have introduced the dimensionless effective rigidity (16)
The contribution of the degreel component to the tidal torque exerted on the solid part is proportional to (e.g. Makarov 2012). For comparison, the imaginary part of the quadrupole Love number is plotted in Fig. 3 for the Andrade and Maxwell rheology (e.g. Makarov & Efroimsky 2013; Correia et al. 2014) as a function of the normalised tidal frequency στ_{M}. The behaviour of the two models does not differ in the lowfrequency regime, where increases with σ linearly. A maximum is reached for P_{T} ~ τ_{M}. In the highfrequency regime, the torque decays as σ increases, but the Andrade model is characterised by a weaker slope than the Maxwell model, as discussed above. Basically, the slope of the imaginary part of the complex degree2 Love number is − 0.25 per decade for the Andrade model (α = 0.25), and − 1 per decade for the Maxwell model (Fig. 3).
Fig. 3 Andrade (solid black line) and Maxwell (dotted blue line) models in frequency domain.The negative imaginary part of the complex quadrupole Love number (see Eq. (15)) characterising a solid planet is plotted as a function of the normalised tidal frequency στ_{M} in logarithmic scales. We used the values of the Andrade parameters given by Table 1 for the Earth. 
Parameters of Andrade model derived from elastogravitational theory (Tobie et al. 2005; Breton et al. 2018).
3.3 Tidal response of the ocean
We now consider the tidal response of the ocean. Similar to the solid part, the oceanic layer undergoes the gravitational tidal forcing of the perturber. It is also affected by the distortion of the solid part, which induces both a variation of the planet selfgravitational potential and a displacement of the oceanic floor, denoted by ξ_{c}. We take this solidocean coupling into account by following along the line by Matsuyama (2014). In this earlier work, the equations of dynamics are written in the shallow water approximation, the fluid being supposedly incompressible. Mean flows are ignored (the whole planet rotates as a solid, the ocean included), and the fluctuations of quantities associated with the perturbation are assumed to be small with respect to background quantities, which is the postulate required by the linear approach.
The perturbation is thus described by an horizontal velocity field V = V_{θ}e_{θ} + V_{φ}e_{φ} and the variation of the oceanic depth η = ξ_{oc} − ξ_{c}, the notation ξ_{oc} designating the oceanic surface displacement (see Fig. 1). In the frame rotating with the planet (), the dynamics of the oceanic tidal response is governed by the momentum equation (Matsuyama 2014), (18)
and the equation of mass conservation, (19)
This should be taken into account with the horizontal gradient operator in spherical coordinates, (20)
the horizontal component of the velocity divergence, (21)
the Rayleigh drag frequency σ_{R} characterising the friction with the oceanic floor and the conversion of barotropic tidal flows into internal gravity waves (Wunsch 1975); and a perturbed potential, denoted by and defined further. This denotation encompasses both the tidal gravitational forcing and coupling with the solid part. We note that the radial component of the Coriolis acceleration − 2ΩV_{φ}e_{r} is neglected in the shallow water approximation since fluid motions are supposed to be dominated by horizontal flows. As a consequence, the Coriolis acceleration in Eq. (18) is simply given by .
In the linear approximation, the response of a given mode is proportional to the corresponding component of the forcing. It follows that Ψ, V, and η can be expanded in Fourier series of t and φ similarly as the tidal gravitational potential (Eq. (3)), (22)
The Fourier coefficients Ψ^{m,σ}, V^{m,σ}, and η^{m,σ} may themselves be expanded as series of the normalised associated Legendre functions , (23)
where , , and designate the degreel components of the expansion.
By introducing the complex tidal frequency and spin parameter , defined by (24)
and substituting perturbed quantities by their Fourier expansions (Eq. (22)) in Eq. (18), we may express the components of the horizontal velocity field as functions of Ψ,
Then, substituting and by the above expressions in Eq. (19), we end up with (27)
the notation referring to Laplace’s tidal operator, and defined by (e.g. Lee & Saio 1997) (28)
We recognise Laplace’s tidal equation in Eq. (27), which determines the horizontal structure of the fluid tidal response (e.g. Lindzen & Chapman 1969). If we assume that the Fourier components Ψ^{m,σ} and η^{m,σ} may be written as (29)
(the parameter n being an integer) and consider that the latitudinal functions are bounded at the poles, Eq. (27) defines an eigenfunctioneigenvalue problem. Hence, the set of eigenfunctions , called Hough functions after Hough’s work (Hough 1898), is associated with a set of eigenvalues through the relationship (30)
We also introduce here the corresponding complex equivalent depths , which are defined by analogy with real equivalent depths (e.g. Taylor 1936; Lindzen & Chapman 1969), by (31)
The Hough function may be written as a combination of the normalised associated Legendre functions, (32)
Conversely, the normalised associated Legendre function may be written as (33)
In the preceding expressions, the and are complex overlap coefficients. They are computed in the meantime as the eigenvalues using the standard method detailed by Wang et al. (2016), which is based on a series of the normalised associated Legendre functions. Additionally, we introduce here the overlap coefficients (34)
which are used to weight the degreen components of the oceanic tidal torque.
In the absence of friction (σ_{R} = 0), corresponds to the resonant configuration where the phase velocity of the forced degreen mode equalises the characteristic propagation velocity of largewavelength surface gravity waves. The nonfrictional case has been thoroughly discussed in early studies (e.g. LonguetHiggins 1968; Lindzen & Chapman 1969; Lee & Saio 1997) so that we do not need to enter into detail here. Some main aspects of this regime are revisited here.
When σ_{R} = 0, the complex quantities defined by Eq. (24) reduce to and , which are the real forcing frequency and spin parameter usually found in the literature. The set of Hough functions divides into two families. The first family, referred to as gravity modes or gmodes (see e.g. Lee & Saio 1997), is defined for . This family corresponds to the ordinary spherical modes modulated by rotation and naturally reduces to the set of the associated Legendre functions in the nonrotating case. For , another family of functions appear. These functions are generally called rotational modes or rmodes, and develop outside of the equatorial band where gmodes are confinedby Coriolis effects.
When a Rayleigh drag is introduced, the sets of Hough functions and associated eigenvalues become complex in the general case (Volland 1974a,b; AuclairDesrotour et al. 2018). The real case discussed above corresponds to the asymptotic regime where . For , the friction affects the behaviour of tidal modes. While the ratios and decay, gmodes and rmodes tend to merge together and converge towards the functions of the nonrotating case, namely the associated Legendre functions. Hence, in the asymptotic limit (), the friction is strong enough to annihilate the distortion caused by the Coriolis effects. From a mathematical point of view, the introduction of friction regularises the solution in the zerofrequency limit. In the absence of dissipation, the number of Hough modes necessary to approximate the solution given by Eq. (29) diverges as (i.e. Ω → 0 or σ → +∞). This is no longer the case when a drag is taken into account.
We now come back to Ψ, which is still to be defined. Following Matsuyama (2014), we introduce the tilt factors associated with the tidal gravitational forcing of the perturber and the distortion of the oceanic layer, denoted by and , respectively. In the framework of the thin shell approximation^{1}, these factors are defined by
They characterise the effective forcing of the oceanic layer including the effects of the tidal distortion of the solid part, which undergoes both the tidal gravitational forcing and the ocean loading. The degreel component of the potential Ψ^{m,σ} is thus defined as (37)
By expanding Ψ^{m,σ} and η^{m,σ} in series of Hough functions, we transform Eq. (27) into (38)
In practice, the series of Hough functions and associated Legendre functions given by Eqs. (23) and (29) are truncated in numerical calculations. By N we denote the number of functions used to approximate series, which is arbitrarily chosen to be large enough to be associated with negligible overlap coefficients. It follows that the of Eq. (29) are the solutions of an algebraic system of the form (40)
In Eq. (40), the σ_{n,k} are the complex characteristic frequencies defined by (41)
where designates the horizontal wavenumber of the degreen mode. The components of the force vector in Eq. (40) are defined by (42)
We note that we have kept components of degrees greater than 2 in the preceding expression for the sake of generality. In reality, these components can be neglected since we assumed that R_{p} ≪ a. Finally, the components of the tidal displacement of the ocean surface in the basis of the normalised associated Legendre functions are simply deduced from the by using the overlap coefficients introduced in Eq. (32), (43)
4 Tidal torque created by eccentricity tides
The modelling of the solid and oceanic tides detailed in the preceding section allows us to determine the tidal torque exerted on the planet. The obtained analytic formulae are used to further compute the evolution of nonsynchronised rotation states of equilibrium with the planet eccentricity and ocean depth.
4.1 General case
In the general case, the solid and oceanic tidal responses are coupled to each other by the gravitational force and ocean loading. In the quadrupolar approximation (R_{p} ≪ a), the tidal torque is thus given by (44)
where is the quadrupolar component of the effective gravitational Love number of the planet, , defined in the thin shell approximation by (45)
In this expression, the degreel component of the tidal gravitational potential and variation of ocean thickness are given by Eqs. (4) and (43), respectively.
4.2 Pure solid tidal response
In the absence of oceanic layer, R_{c} = R_{p} and the tidal torque exerted on the planet with respect to the spinaxis is simply expressed as (Efroimsky & Williams 2009; Makarov 2012; Makarov & Efroimsky 2013; Correia et al. 2014) (46)
In this case, the tidal response of the planet is characterised by the solid gravitational Love number given by Eq. (15). The expression of this parameter in the framework of the Andrade model may be found in Eq. (6) of Makarov (2012) for instance. As discussed in Sect. 3.2, the Andrade and Maxwell times generally far exceed typical tidal periods. As a consequence the imaginary part of may be approximated by (47)
except in the zerofrequency limit, where tidal periods become comparable with the Andrade and Maxwell times in order of magnitude.
4.3 Pure oceanic tidal response
The tidal response of the oceanic layer is more complex than that of the solid part. To examine it, we have to consider the particular case where the coupling with the solid part vanishes. In this simplified case, the rigidity of the solid body is supposed to be infinite (μ → +∞), so that the Love numbers defined by Eq. (15) are annihilated. Moreover, the Cowling approximation (e.g. Cowling 1941; Unno et al. 1989) is assumed, which means that the term in ρ_{oc} ∕ρ_{c} resulting from the selfattraction of the ocean in Ψ is ignored. In this case, the oceanic tidal torque may be written similarly as that of the solid part (Eq. (46)), (49)
where designates the degree2 tidal gravitational Love number of the ocean.
Owing to the Cowling approximation, σ_{n,k} = 0 for n≠k. The matrix of the algebraic system given by Eq. (40) is thus diagonal, which leads to the solutions obtained in early studies (e.g. Webb 1980; AuclairDesrotour et al. 2018), (50)
the notation and designating the complex eigenfrequencies of the degreen Hough mode, defined by (see e.g. Webb 1980, Eq. (2.12) therein) (51)
In this case, the quadrupolar Love number introduced in Eq. (49) is simply expressed as (52)
Using the mean density of the planet (), it may also be put into the form (53)
which highlights the fact that a resonance occurs when the equivalent depth of a given mode equalises the ocean depth in the absence of friction.
The friction of tidal flows against the oceanic floor affects both the overlap coefficients and the equivalent depths of Hough modes. Since these parameters are complex in the general case, characterising the dependence of the tidal torque on σ_{R} is not straightforward. However, this may be done in the quasiadiabatic asymptotic regime (), where the imaginary part of is negligible, and for gmodes, given that the latter are associated with positive . By introducing the characteristic frequency of the degreen surface gravity wave, (54)
we express the imaginary part of the oceanic Love number as (55)
The preceding expression characterises the shape of a resonant peak in the frequencyspectra of the tidal torque. The frequency atwhich the resonant peak of a mode reaches a maximum is expressed as (56)
In the quasiadiabatic regime (), it reduces to (57)
The corresponding maximum of for the resonance associated with the degreen Hough mode is denoted by and obtained by substituting σ with Eq. (57) in Eq. (55). If the contributions of nonresonant components are neglected, we thus obtain (58)
which shows that, as a first approximation, the maximum of the peak associated with the degreen Hough mode varies with the ocean depth as and is inversely proportional to the typical frequency characterising the Rayleigh drag, in agreement with scaling laws derived using simplified Cartesian models (e.g. Auclair Desrotour et al. 2015).
When the frequency associated with an eccentricity term becomes equal to the eigenfrequency of the resonant mode, this eccentricity term is enhanced by a factor corresponding to the ratio of the peak maximal value over the level of the nonresonant background.It can thus generate a torque strong enough to compensate the predominating contribution of the semidiurnal component in the loweccentricity regime.
This mechanism is illustrated by Fig. 4, which represents the main features of the tidal torque as functions of the semidiurnal frequency, , in the supersynchronous regime (Ω > n_{⋆}). In this figure, the tidal response of the solid part is ignored in order to focus on the oceanic tidal torque, which is reduced to the contributionof the degree0 Hough mode for simplification. For pedagogical purposes, we ignore the effect of rotation, which induces Rossby modes through the Coriolis terms in the momentum equation. This allows us to assume that the overlap coefficients and eigenvalues associated with the considered Hough mode do not vary with the tidal frequency, for a fixed ocean depth.
The simplified oceanic tidal response is thus solely due to the prograde and retrograde degree0 Hough modes. The corresponding frequencies of these resonant modes are ±σ_{0} and the resulting tidal torque contributions, which both drive the system towards synchronism, are represented schematically by blue cones in Fig. 4. As Eq. (54) shows, the frequencies of these intrinsic modes of the ocean increase with the ocean depth H_{oc}. Indeed, Fig. 4 shows that the blue cones’ positions shift towards the higher frequencies when H_{oc} increases from panels a to f.
One of these peaks is located in the supersynchronous frequency range (σ_{2,2} > 0) and the other one in the subsynchronous frequency range (σ_{2,2} < 0). Synchronisation is symmetrical and is due to the nonrotating approximation assumed in Fig. 4. We emphasise the fact that the peak of positive torque appearing in the supersynchronous regime results from the forcing of the resonance located in the subsynchronous frequency range by a prograde eccentricity tidal potential (σ_{2,s} = −σ_{0}). This configuration corresponds for instance to case 1 of Fig. 2, where the degree3 eccentricity potential propagates eastward.
In this framework, Fig. 4 simultaneously shows the forcing frequencies σ_{2,s} = 2Ω − sn_{⋆} of the semidiurnal mode (s = 2) and of two eccentricity modes (s = 3 and s = 4) in the cases where a peak of positive torque is generated by one of these eccentricity modes in the interval 0 < σ_{2,2} < n_{⋆}. These particular frequencies are symbolised by the superscript ▲, and the created peak is designated by a red cone. Hence, the contribution of one of the eccentricity components potentially leads to asynchronous equilibrium states located at . This configuration occurs when the peak is sufficiently important to counterbalance the semidiurnal tidal torque, whose resonances are represented by blue cones as we plot the schematic torques as a function of σ_{2,2}.
For the cases a–c of Fig. 4, the assumption made about (i.e. ) leads to the fact that only the degree3 eccentric term can coincide with the prograde resonant frequency of the ocean (subsynchronous blue cone), that is . When the depth of the ocean increases (from d to f), keeping in the range 0n_{⋆}, this leads to an excitation of the ocean resonant mode only with the degree4 eccentricity mode (for ). This illustrates how increasing the depth of the ocean leads to the excitation of higher eccentricity modes in the tidal potential.
Configuration (a) represents the case where σ_{0} = n_{⋆}∕2. In this case, the peaks created by the resonance of the degree3 and degree2 (semidiurnal) component are exactly superposed. The semidiurnal one is the strongest in the low eccentricity regime since . As a consequence, the oceanic tide leads the planet toward synchronisation (Ω = n_{⋆}) in this configuration, and cannot generate an asynchronous rotation state of equilibrium in the range − n_{⋆} < σ_{2,2} < n_{⋆}.
While the ocean depth increases (configurations (b) and (c)), the resonance moves away from the synchronous rotation and becomes stronger, its eigenfrequency and maximum value scaling as (see Eqs. (54) and (58)). The peak generated by the degree3 eccentricity component increases as well in intensity and gets closer to synchronisation in the meantime. This means that the resulting asynchronous state of equilibrium – if it exists – also gets closer to synchronisation until being annihilated by the solid tidal torque, which outweighs the oceanic torque when P_{T} ~ τ_{M}.
Configuration (d) shows the switch from the degree3 term to the degree4 term that occurs while the ocean depth keeps increasing. The peak of positive tidal torque now results from the amplification of the degree4 eccentricity term by the resonance associated with the tidal oceanic mode. It is smaller than that generated by the degree3 term since it scales quadratically with the forcing tidal potential (see Eqs. (7) and (49)), which decays while the degrees component increases in the loweccentricity regime (for e ≳ e_{trans}, nonlinearities lead eccentricity components of higher degrees to predominate; see e.g. Ogilvie 2014, Fig. 3).
Similar to the peak created by the degree3 component, the peak associated with the degree4 eccentricity term moves towards synchronisation while the ocean depth increases, as shown by configurations (e) and (f). When σ_{0} becomes greater than2n_{⋆}, the term generating the peak of positive tidal torque switches from the degree4 to the degree5 eccentricity term, and so on, as discussed at the end of Sect. 3.1 (see Eq. (11)).
The mechanism highlighted by Fig. 4 for the degree0 Hough mode can be generalised to other modes, each of them being able to amplify an eccentricity term provided that the associated resonance dominates the nonresonant background level. This picture is also completed by the action of solid tides, which enable the existence of tidallylocked asynchronous rotation states of equilibrium for the rotation rates , as discussed in Sect. 5.
In reality, the simplified oceanic tidal response is not solely due to the degree0 Hough mode of frequency σ_{0} since in the frequency range of interest. The real tidal response is composed of g and rmodes (i.e. the gravity modes modified by rotation and the Rossby modes, respectively), the latter modes being restored by the spin rotation of the planet. Therefore, the resulting behaviour is much more complex than that described in Fig. 4 although it exhibits the main highlighted features.
In the quasiadiabatic regime, the maximum of the tidal torque given by Eq. (58) can be arbitrarily high depending on the value of the Rayleigh drag frequency. When σ_{R} → 0, this maximum tends to infinity. This is an artefact of the linear theory. In reality, the large amplitudes of tidal fields associated with a resonance in the quasiadiabatic regime violate the small perturbation approximation. As a consequence, the tidal response becomes nonlinear and the maximum is attenuated with respect to that predicted by the linear theory.
In the case of the Earth, the Rayleigh drag frequency used to model the friction with the oceanic floor was estimated at σ_{R} ≈ 10^{−5} s^{−1} (e.g. Webb 1980), which is comparable with typical tidal frequencies. In this case, the enhanced oceanic tidal response resulting from a resonance remains in the framework of the linear theory, and the associated amplification of the tidal torque generally does not exceed a decade in logarithmic scale (see AuclairDesrotour et al. 2018, Fig. 5).
As the effective Rayleigh drag frequency accounts here for the amount of energy dissipated through the interaction of tidal waves with the Earth complex topography, we assume that 10^{−5} s^{−1} is an upper estimation for σ_{R}. It seems likely that σ_{R} takes smaller values in the case of a planet hosting a uniform ocean. To take these discrepancies into account, we consider the values 10^{−7} and 10^{−6} s^{−1} in addition to 10^{−5} s^{−1} in numerical calculations.
In spite of the limitations mentioned above, studying the quasiadiabatic regime is a useful step to understand how the resonances characterising the oceanic tidal response can enhance the desynchronising effect of eccentricity tides. This is thus the framework that we adopt to derive a theoretical estimation of the critical eccentricity in Sect. 5.
Fig. 4 Positive tidal torque generated by oceanic eccentricity tides in supersynchronous frequency range (Ω > n_{⋆}) for increasing oceanic depths (a–f). An oceanic mode is characterised by a resonance of eigenfrequency σ_{n} (blue peaks). We focus here on the degree0 mode, of frequency σ_{0}, in the nonrotating approximation (Coriolis effects are ignored). An eccentricity term of frequency σ_{2,s} = 2Ω − sn_{⋆} generates a positive torque, which is enhanced by the resonance when (red peak). The semidiurnal (σ_{2,2}) and first eccentricity frequencies (σ_{2,3} and σ_{2,4}) associated with the peak thus created are symbolised by the superscript ▲. The peak raised by the eccentricity tide occurs when , which also corresponds to and . The direction of peaks indicates the sign of the tidal torque. In the supersynchronous regime (σ_{2,2} > 0), downwardblue peaks tend to drive the planet towards tidal locking in spinorbit synchronous (σ_{2,2} = 0), while upward red peaks tend to drive it away from this state of equilibrium. 
4.4 Frequency behaviour of tidal torque
The tidal torque resulting from the planet tidal response is characterised by complex frequency behaviour, which slightly diverges from the idealised picture shown by Fig. 4. This behaviour results both from the gravitational and surface coupling between the solid and oceanic layers, and from the dependence of the oceanic tidal response on Coriolis effects in the frequency range of interest.
It seems helpful to focus on a given case to visualise how the tidal torque exerted on the planet varies with the tidal frequency before going further in the analysis. As shown previously, the amplification of eccentricity terms by resonances associated with the oceanic tide most likely takes place at short orbital periods (P_{⋆} ~ 1–10 days) since this regime corresponds to H_{oc} ~ 1–10 km. This regime is typically encountered in tightlypacked systems such as that hosted by the TRAPPIST1 ultracool dwarf star (Gillon et al. 2017; Grimm et al. 2018), where planet e exhibits a 6.10day orbital period (Gillon et al. 2017).
We thus consider throughout the whole study the case of an Earthsized planet orbiting TRAPPIST1 (M_{⋆} = 0.09 M_{⊙}, see Van Grootel et al. 2018) with a sixday orbital period. We note that the mass of the host star is the only parameter of the system here, since the internal structure and surface conditions of the planet are fixed. For illustrative purposes, the planet is assumed to host a 0.2 kmdeep uniform ocean, which corresponds to the resonance of one of the dominating oceanic modes for the degree3 eccentricity term in the range 0 < σ_{2,2} < n_{⋆}. We set the typical frequency of the Rayleigh friction to σ_{R} = 1.0 × 10^{−6} s^{−1}. Finally, the response of the solid part is described by the values of Andrade parameters given by Table 1, as mentioned in Sect. 3.2.
The tidal torque exerted on the planet is plotted in Fig. 5 as a function of the normalised semidiurnal tidal frequency , for circular (e = 0) and eccentric (e = 0.2) orbits. Panels from top to bottom correspond to the case of a pure solid tide (H_{oc} = 0), the case of a pure oceanic tide (μ = +∞), and the general case, respectively. In the three cases, the tidal torque is calculated using the expression given by Eq. (44), with H_{oc} = 0 or μ = +∞ to obtain the pure solid and oceanic response, respectively. To highlight the effect of eccentricity components, regimes where the planet is driven towards or away from the spinorbit synchronous rotation state of equilibrium are designated by solid or dashed lines.
First, we consider the torque generated by a pure solid tide (Fig. 5, top panel). We retrieve here the behaviour investigated by Makarov & Efroimsky (2013) and Correia et al. (2014). In the absence of eccentricity, the semidiurnal component solely generates a torque, which varies with the frequency following the Andrade frequencydependence illustrated by Fig. 3. For tidal periods smaller than the Maxwell and Andrade times, we recover the scaling law given by Eq. (47). This torque drives the planet towards the spinorbit synchronisation. A sufficiently large eccentricity induces additional asynchronous rotation states of equilibrium associated with the rotation rates (Makarov &Efroimsky 2013; Correia et al. 2014), which are discussed in Sect. 5.
We then move to the pure oceanic tidal response (Fig. 5, middle panel), where the rigidity of the solid part is supposed to be infinite. In this case, the tidal torque is composed of a batch of resonances associated with the main Hough modes coupled with the tidal gravitational potential. In the circular configuration, two dominating peaks may be observed. For e = 0.2, additional peaks appear. Each of these peaks is generated by the enhanced contribution of an eccentricity term. For instance, the peak noticed for clearly corresponds to the configuration illustrated by Fig. 4b, where the frequency associated with the degree3 eccentricity term meets the eigenfrequency of one of the dominating modes.
When the typical tidal periods of predominating tidal components exceed the characteristic timescale of inertial waves and are less than or comparable with the energy decay timescale of the ocean in the meantime (), the tidal forcing tends to couple with an infinite number of Hough modes as increases (this corresponds to the regime of subinertial waves, defined by , and discussed by, e.g. Lee & Saio 1997). The method adopted to solve Laplace’s tidal equation (see Eq. (30)) in this study cannot correctly treat this phenomenon because the number of computed modes is fixed by the dimension of the truncated matrix used in calculations (Wang et al. 2016). The effect of truncation becomes more important as increases.
As a result, the oceanic tidal response is poorly described by the theory for , which degrades the estimation of the induced tidal torque and leads to the unrealistic discontinuous change of sign observed for the circular case in the vicinity of the synchronisation (Fig. 5, middle panel). This feature of the model prevents us from exploring the region of the parameter space characterised by very small drag frequencies. Fortunately, the action of the drag annihilates the distortion caused by the spin rotation by making Hough functions converge towards the associated Legendre functions, as discussed in Sect. 3.3. The change of sign thus does not occur from the moment that , which is the case in this study if σ_{R} ≳ 10^{−6} s^{−1}.
We finally consider the tidal torque in the general case, where the tidal responses of the solid and oceanic layers are coupled together by gravitational and pressure forces (Fig. 5, bottom panel). We mainly retrieve the frequencyresonant behaviour of the ocean, but the tidally dissipated energy is attenuated by one order of magnitude. This is due to the elastic adjustment of the deformable solid part, which tends to compensate the horizontal gradient of mass distribution associated with the elevation of the ocean surface level. As we can now better visualise the dependence of the tidal torque on the forcing frequency, we can proceed to an analytical estimation of the critical eccentricity beyond which asynchronous rotation states of equilibrium may exist.
Fig. 5 Logarithm of tidal torque exerted on planet. In the cases of the pure solid (top panel) and oceanic (middle panel) tidal responses, and in the general case (bottom panel), these are shown as functions of the normalised semidiurnal frequency for a circularorbit (red line) and an orbit of eccentricity e = 0.2. The torques , , and are computed using the expression given by Eq. (44), and the values given by Table 1 for an Earthsized planet orbiting TRAPPIST1 (M_{⋆} = 0.09 M_{⊙}) with a sixday orbital period. The Rayleigh drag frequency is set to σ_{R} = 1.0 × 10^{−6} s^{−1}. Solid (dashed) lines designate the regimes where the tidal torque drives the planet towards (away from) spinorbit synchronousrotation. 
5 Analytical estimation of critical eccentricity for an ocean planet
In the previous sections, we detailed the formalism that describes the tidal response of an ocean planet by including the interactions between the solid part and the oceanic layer, and derived the expressions of the planetary tidal torque and second order Love number. First, we use here these solutions to establish the conditions that have to be satisfied to end up with tidallylocked asynchronous rotation states of equilibrium. Second, we derive an analytical estimation of the critical eccentricity beyond which such states can exist in the low eccentricity regime. The obtained results are used in Sects. 6 and 7.
In the nonresonant case, eccentricity tides may drive the planet away from synchronisation only if the associated components of the tidal gravitational potential are comparable with or greater than the component associated with the semidiurnal tide, that is for e ≳ e_{trans} ~ 0.25 typically, as shown in Sect. 3. The critical eccentricity e_{AR} may be strongly lowered when a resonance occurs, since this later generates an important increase of the tidally dissipated energy. This mechanism enables the existence of asynchronous rotation states of equilibrium at smaller eccentricities. In this section, the planet is considered an idealised rotationally symmetric body, and the triaxial torque due to triaxiality is thus ignored. The action of triaxiality regarding the trapping of the planet in spinorbit resonances is investigated in Sect. 6.2.
If we assume that e ≪ 1, three conditions must be satisfied at the same time to lock the planet into a nonsynchronised rotation rate: (i) the forcing frequency associated with one of the eccentricity components of the response (see Eq. (49)) must be equal to one of the eigenfrequencies of the oceanic surface gravity modes, given by Eq. (54); (ii) the resonant eccentricity component must drive the planet away from synchronisation as discussed in Sect. 3; and (iii) the resonance must dominate the nonresonant background, which is the sum of nonresonant terms. We examine these conditions for eccentricity components of degrees s ≥ 3, given that they predominate.
Condition (ii) is expressed by Eq. (11). In this configuration, the degrees forcing frequencies are negative while the semidiurnal forcing frequency is positive. Assuming condition (i), we suppose that the degreen mode is resonant. This leads to (59)
where σ_{2,s} = 2Ω − sn_{⋆} is the degrees eccentricity frequency defined by Eq. (6) and σ_{n} the eigenfrequency of the degreen Hough mode given by Eq. (54). By combining Eqs. (11) and (59), and introducing the notation , we hence derive from (i) and (ii) the first condition on the ocean depth (H_{oc}), which must satisfy the inequality (60)
This inequality determines the range of depths for which a resonance may occur. We observe that the range of H_{oc} widens with s, meaning that the smaller the ocean depth, the larger the number of eccentricity components that may be excited by the resonance. Similarly, the upper limit of H_{oc} increases with n_{⋆}, which favours resonances in the case of closein terrestrial planets. Conversely, this limit is inversely proportional to the eigenvalue of the mode. In the quasiadiabatic regime, the are sorted in ascending order (e.g. Lee & Saio 1997). Thus, the upper limit of H_{oc} is lower for high degrees n than for low degrees. Basically, this means that the resonance of the degree0 mode will be expressed preferentially with respect to those associated with other modes. A highdegree mode cannot enter into resonance except in the case of a very thin ocean.
Establishing an analytical formulation of condition (iii) is challenging – if not impossible – in the general case because of the coupling between the solid and fluid layers. We thus choose to ignore this coupling as well as the effect of the variation of planet selfattraction (this refers to the Cowling approximation mentioned above). Moreover, as we focus on the case of small eccentricities (e ≪ 1), the major part of the nonresonant background level results from the semidiurnal component (see Eq. (8)), which tends to drive the planet towards synchronisation. We thus neglect other terms and assume that the nonresonant part of the response reduces to the semidiurnal tide.
The case of a purely solid planet has been studied by Makarov & Efroimsky (2013) and Correia et al. (2014), who used Maxwell (e.g. Greenberg 2009) and Andrade (e.g. Efroimsky 2012) rheologies to model the frequency behaviour of the body, respectively.These early works show that eccentricity tides drive the planet towards final rotation rates that are multiples of n_{⋆} ∕2 in the regime where P_{T} ≪ τ_{M} and for sufficiently high eccentricities. In this case the spin equilibria occur when the torque associated with one of the eccentricity components becomes resonant (see Fig. 3). These positions approximately correspond to the rotation rates for which the forcing frequencies annihilate, that is . For small parameters α of the Andrade rheology (typically α = 0.25), the order of magnitude of the tidal torque does not vary much with the tidal frequency. As a consequence, asynchronous spin equilibria may not exist at small eccentricities, if the departure between the frequencies of the semidiurnal and degree3 components (that is n_{⋆}) is not large enough.
The sharp variations of the tidal torque associated with the resonances characterising the ocean tidal response enable asynchronous rotation state of equilibrium to exist at small eccentricities. If we neglect the coreocean coupling resulting from gravitational and surface forces, condition (iii) may be simply written as (61)
where σ_{2,s} = 2Ω − sn_{⋆} = −σ_{n} (see Eq. (59)) designates the resonant forcing frequency of the degrees quadrupole component (the resonance being associated with the degreen Hough mode), and the semidiurnal frequency.
First, we consider the case where the oceanic tide always predominates over the solid tide, which amounts to assuming that the solid part of the planet is of infinite rigidity (μ → +∞). In this case, . The associated Legendre function is generally well coupled with the degree0 Hough mode, which means that except for n = 0. As a consequence, the total tidal torque is dominated by the contribution of the degree0 Hough mode outside of the resonances associated with other modes (AuclairDesrotour et al. 2018), and may be reduced to this component as a first approximation.
It follows that the condition given by Eq. (61) is simply expressed as (62)
To obtain this expression, we have assumed that the degrees eccentricity component (i) is resonant, that is σ_{2,s} = −σ_{0} and (Eq. (59)); (ii) made use of Eqs. (55) and (58) for the degree0 Hough mode; and (iii) introduced the notation .
The condition on the eigenfrequency σ_{0} resulting from the preceding expression is implicit in the general case since the overlap coefficients ( and ) and eigenvalue () associated with the degree0 Hough mode both depend on the forcing frequencies σ_{2,s} and σ_{2,2}. This dependence can be neglected in the regime of rapid rotation () since spin parameters determining the solutions of Laplace’s tidal equation (Eq. (30)) hardly vary in this case. In the general case (), the dependence of the spin parameters on the forcing frequencies is both stronger and much more complex, as shown for example by Lee & Saio (1997) and Townsend (2003). It is fully taken into account in the numerical calculations performed in Sects. 6–8.
Nevertheless, the fact that analytic solutions of Laplace’s tidal equation do not exist for the general case leads us to ignore this dependence here, and we obtain thereby that condition (iii) is satisfied for (63)
where we have introduced the supposed constant parameter (65)
Equations (63) and (64) determine the interval of H_{oc} in which the semidiurnal torque dominates the torque induced by the degrees eccentricity component, and thus leads to synchronisation. The width of this interval depends on the parameter γ, which compares the Rayleigh friction frequency to a characteristic frequency proportional to the degrees eccentricity component of the tidal gravitational potential. By recalling Eq. (60), we hence observe that the degrees eccentricity term cannot drive the planet away from spinorbit synchronous rotation if γ_{s} ≫ 1 since the torque induced by the semidiurnal component is systematically stronger.
Asynchronous rotation rates may exist only if γ_{s} ≲ 1. In the asymptotic limit (γ_{s} → 0), the eccentricity component always predominates when its forcing frequency is equal to the eigenfrequency of the degree0 Hough mode (σ_{0}). In this case, the only condition that has to be satisfied to make possible the existence of asynchronous final rotation rates is the condition given by Eq. (60). This highlights the criterion for which resonances may enable eccentricity terms to predominate at small eccentricities, (66)
When this criterion is satisfied, the amplifying effect of the resonance compensates the smallness of the eccentricity forcing component with respect to the semidiurnal one.
By focusing on the degree3 component, which is the stronger in the low frequency limit (e.g. Ogilvie 2014), and using Eq. (10), we immediately derive from Eq. (66) an analytical estimation of the critical frequency beyond which tidallylocked asynchronous states of equilibrium can exist, (67)
This expression shows that the critical eccentricity is actually determined by only two parameters when the tidal response of the planet is reduced to the oceanic tide: the orbital and Rayleigh drag frequencies. For a typical orbital period P_{⋆} = 6 days, n_{⋆} = 1.21 × 10^{−5} s^{−1}. This leads to e_{AR} ≈ 0.1 for σ_{R} = 10^{−5} s^{−1} and e_{AR} ≈ 0.01 for σ_{R} = 10^{−6} s^{−1}. We retrieve these orders of magnitude in the numerical results detailed in Sect. 6.
We investigated the above case where the oceanic tidal torque always predominates over the torque exerted on the solid part. We now consider the case where the solid tidal torque predominates outside of the resonances of the oceanic tidal response, the nonresonant component of the oceanic tidal torque being assumed negligible compared to the solid tidal torque. Condition (iii) is thus expressed as (68)
By assuming that , considering that the degrees is amplified by the resonance associated with the degree0 Hough mode, and combining together Eqs. (47) and (58), we put the preceding equation into the form (69)
where we have introduced the normalised frequency X = σ_{0}∕n_{⋆} and the dimensionless constant (70)
Similarly as what we did to establish Eq. (67), we suppose that the oceanic tidal response behaves similarly to the nonrotating case, and ignore the implicit dependence of the overlap coefficient on the spin rotation (Ω) and forcing (σ_{2,s}) frequencies through the spin parameter (). We thus assume that . Moreover, as we are interested in the loweccentricity regime, we focus on the effect of the degree3 eccentricity term.
In the general case, the two roots defined by the inequality given by Eq. (69) cannot be derived analytically since α is real. However, the calculation of the maximum of the polynomial function corresponding to the lefthand member of the equation is straightforward. Denoting this maximum as u_{α} we obtain, for s = 3, (71)
Since the inequality given by Eq. (69) can be satisfied only if , it follows that an estimation of the critical eccentricity can be derived by substituting in this latter inequality. We end up with (72)
The comparison between the two obtained expressions of e_{AR}, given by Eqs. (67) and (72), highlight the role played by the tidal response of the solid part, which significantly attenuates the sensitivity of the studied amplification mechanism to the strength of the drag. This sensitivity is important in the regime where the oceanic tide is the predominating source of the tidally dissipated energy. In this regime, the critical eccentricity for asynchronous states of equilibrium scales as e_{AR} ∝ σ_{R}. When the solid tide determines the nonresonant background level, this scaling law switches to . As the nonresonant background level of the oceanic component of the tidal torque scales as ∝ σ_{R} (see e.g. AuclairDesrotour et al. 2018, Fig. 5), this means that the critical eccentricity ceases to decay proportionally with σ_{R} from the moment that σ_{R} becomes smaller than a critical value.
The above analysis provides a reference basis to interpret numerical results. It enables us to go deeper into the exploration of the parameter space with the calculation of the planet finalrotation state as a function of its eccentricity and ocean depth.
6 Evolution of planet final rotation with eccentricity and ocean depth
In this section, we investigate how the finalrotation rate of a planet evolves with its eccentricity and ocean depth. First, we focus on the equilibria determined by the tidal torque. Second, we discuss the existence of pseudosynchronous equilibria close to the synchronisation by considering the role played by the triaxial torque – that is the torque due to the planet inherent and permanent triaxiality – and using the theory of capture in spinorbit resonances (Goldreich 1966; Goldreich & Peale 1966).
6.1 Tidal locking in asynchronous states of equilibrium
We treat the cases of an Earthsized planet and of a 10M_{⊕} superEarth orbiting the TRAPPIST1 dwarf star (M_{⋆} = 0.09 M_{⊙}). In the first case, we consider two orbital periods: P_{⋆} = 1 day and P_{⋆} = 6 days. In the second case, the superEarth orbits the star with a sixday orbital period. As the Rayleigh drag frequency cannot be accurately specified for an exoplanet of unknown topography, we perform the calculations for σ_{R} = 10^{−7}, 10^{−6}, 10^{−5} s^{−1}, the lower and upper values corresponding to weak and strong drags, respectively. The values of Andrade parameters used for the two planets are those specified in Table 1.
For given eccentricity and ocean depth, the final rotation of a planet is determined by successive iterations from an initial value, specified thereafter, with a frequency step scaling as ∝ σ_{R} in order to take into account the dependence of resonance peak widths on σ_{R}. At each step of the research, the sign of the tidal torque exerted on the planet is computed with TRIP (Gastineau & Laskar 2011) using the expression derived in the general case, and given by Eq. (44). As the research evolves following the direction defined by the sign of the tidal torque, the final rotation state of equilibrium corresponds to a change of sign. When this change of sign is encountered, the finalrotation state is determined with an arbitrary precision using a combination of the dichotomy and secant methods. As mentioned in the introduction, we ignore for the moment the torque due to the planet triaxiality and consider that the final state of equilibrium is solely determined by the annihilation of the tidal torque exerted on the planet.
We focus on the existence of asynchronous rotation states in the interval 0 < σ_{2,2} < n_{⋆} (see Fig. 4), which corresponds to the supersynchronous frequency range (Ω > n_{⋆}). In this interval, the tidal torque generated by the solid part in the absence of rotation reaches its maximal values for and , that is, in the vicinity of the interval bounds. As a consequence, if the initial value of the semidiurnal tidal frequency is too close to 0 or n_{⋆}, the solid tide predominates and the final state of equilibrium is either the spinorbit synchronous rotation or the asynchronous rotation state induced by the degree3 eccentricity term, (see discussion in Sect. 5).
As a consequence, the initial semidiurnal frequency has to be set sufficiently far from the bounds of the studied frequency interval to observe the effect of oceanic tides on the final rotation state of equilibrium. In this work, we arbitrarily start the research at the middle of the interval, that is with the initial semidiurnal frequency σ_{2,2} = n_{⋆}∕2. This corresponds to the initial rotation rate .
Results are plotted in Fig. 6. Each panel of the figure represents the normalised finalrotation period P_{rot}∕P_{⋆} as a function of the eccentricity and ocean depth in logarithmic scales for a given case defined by the type of planet (Earth or superEarth), the orbital period (P_{⋆}), and the Rayleigh drag frequency (σ_{R}). We note that, although global oceans with depths of ≲0.1 km are very unlikely owing to the extremely flat topography it would imply, the lower bound of is set to in order to emphasise the asymptotic regime of dry planets. The yellow colour designates the region of the parameter space where the final state is the spinorbit synchronous rotation, while shades of blues indicate asynchronous rotation states of equilibrium such that . The six topand middle panels correspond to the Earth, the three bottom ones to the 10M_{⊕} superEarth.
Numerical calculations highlight the mechanism discussed in the previous sections: a resonance associated with an oceanic mode can decreasethe lower eccentricity of the region where asynchronous final states exist. This region approximately corresponds to log_{10}e ≳−0.6 (i.e. e ≈ 0.25) in absence of ocean, that is in the asymptotic limit H_{oc} → 0. We note that this value is approximately the transition eccentricity e_{trans} at which the degree3 eccentricity tidal potential becomes equal to the semidiurnal tidal potential (see Sect. 3.1).
As the depth of the ocean increases, the eigenfrequency of the predominant oceanic mode first equalises the frequency associated with the degree3 eccentricity term, thus leading to the most important decay of the critical eccentricity. It then meets the frequency of the degree4 eccentricity term, and so on, as illustrated by Fig. 4. This generates a series of peaks located at the corresponding resonant depths, denoted by H_{oc ;s}. Moreover, the asynchronous rotation state of equilibrium resulting from a resonance tends to get closer to the synchronisation while H_{oc} increases, since it follows the decay of the resonant eigenfrequency (see Fig. 4).
Owing to the scaling law given by Eq. (60), the resonant ocean depths increase when the orbital period switches from six days to one. The smaller the orbital period, the deeper the ocean likely to be subject to resonances. A decay of the Rayleigh drag frequency (from left to right) accentuates the observed features, and particularly the decay of the critical eccentricity in resonant configurations as discussed in Sect. 5 (see Eqs. (67) and (72)). In these cases, the value of the critical eccentricity drops from e_{AR} ≈ 0.25 to e_{AR} ≈ 0.015 (see e.g. the top right and middle panels of Fig. 6). In the most extremal treated case (P_{⋆} = 1 day and σ_{R} = 10^{−7} s^{−1}, middle right panel), the lowest critical eccentricity reached is e_{AR} ≈ 0.006. However, this case is unlikely since it clearly goes beyond the scope of the linear theory.
By considering the cases defined by a sixday orbital period and the Rayleigh frequencies 10^{−5} and 10^{−6} s^{−1} (top left and middle panels), we observe the transition from the frictional regime () to the quasiadiabatic regime (). In the first configuration, the ocean depth minimising the critical eccentricity is affected by the drag, while it is not the case for σ_{R} ≲ 10^{−6} s^{−1}.
Finally, we note that the results obtained for the superEarth are practically identical to those obtained for the Earth. This is mainly due to the fact that the oceanic tidal response generally predominates over the tidal response of the solid part, except in a very small vicinity of the lower and upper bounds of the studied frequency interval. Thus, the final rotation of the planet is not very sensitive to the properties of the solid part, although it is strongly affected by its tidal response (see Fig. 5). Besides, the scale factors and associated respectively with the eigenfrequency of a mode – given by Eq. (54) – and the corresponding maximum of the tidal torque – given by Eq. (58) – are approximately the same for the two planets. As a consequence, variations related to the change of mass and radius from a case to another are small.
In the dry asymptotic limit (H_{oc} → 0), one might expect that the critical eccentricity depends on the Andrade parameters characterising the solid part of the planet, which would lead to two different values for the Earth and the superEarth. This is not the case however, given that e_{AR} ≈ 0.25 for both planets. To explain this feature, we come back to the values given by Table 1. The Maxwell and Andrade timescales of both planets far exceed typical tidal periods. As a consequence, the imaginary part of the degree2 Love number is well approximated by its highfrequency asymptotic functional form, given by Eq. (47). This allows us to factorise all eccentricity terms (Eq. (46)) by (73)
and write the tidal torque exerted on the solid body as (74)
As may be noticed, the sign of the torque in this equation only depends on the rheological exponent α, set to 0.25 for both planets, and the eccentricity tidal frequencies and Hansen coefficients. Thus, the critical eccentricity is not affected at all by the other parameters of the solid body, except in the close vicinity of the resonances associated with eccentricity terms. This explains why no difference can be detected between the Earth and the superEarth in the dry asymptotic regime (H_{oc} → 0) in Fig. 4 (top and bottom panels).
The above observations suggest that final rotation states of equilibrium distinct from spinorbit resonances () are not very sensitive to most of the parameters of the solid part (size, mass, bulk rigidity, Maxwell and Andrade times), although they are affected by the rheological exponent (α) in the higheccentricity regime. The important parameters for these states are the ocean parameters (depth, density, and Rayleigh drag frequency) and the orbital period.
Fig. 6 Normalised finalrotation period P_{rot}∕P_{⋆} as a function of eccentricity (horizontal axis) and ocean depth (vertical axis) in logarithmic scales. Two planets are considered: an Earthsized planet with sixday (top panels) and oneday (middle panels) orbital periods, and a 10 M_{⊕} superEarth with a sixday orbital period (bottom panels). From left to right: configurations associated with the values of the Rayleigh drag frequency σ_{R} = 10^{−5}, 10^{−6}, 10^{−7} s^{−1}. Values of parameters used for the solid part are given by Table 1. In all cases, the stellar mass is M_{⋆} = 0.09 M_{⊙}, the initial rotation rate is (middle of the interval 0 < σ_{2,2} < n_{⋆}, see Fig. 4), and 10 Hough modes are taken into account in the calculation of the tidal response. The degrees of eccentricity modes causing resonances are indicated in the case of the sixdayperiod Earthsized planet with a weak drag(top right panel), consistent with the discussion on Fig. 4. 
6.2 Capture in 1:1 spinorbit resonance
In reality, the final rotation state of equilibrium of the planet is not determined by the tidal torque solely, but also depends on its triaxiality, that is, its permanent, nonaxisymmetric deformation. As highlighted by Goldreich (1966) and Goldreich & Peale (1966), the triaxiality of the planet creates an additional gravitational torque, which tends to trap it in spinorbit resonances.
In the framework of the coplanar case, the torque due to triaxiality (noted as ) adds to the tidal torque (noted as the superscript “tide” here to avoid confusion with the total torque) in the righthand member of the equation of the spin evolution, (75)
This torque is expressed in the general case as (see Appendix C) (76)
In the preceding equations, the symbol ^{∙} is the time derivative, A, B, and C the principal moments of inertia of the planet (in increasing magnitude), and the angles associated with the eccentricity frequencies σ_{2,s}, with the notation ϑ and M designating the planet rotation angle and mean anomaly, respectively. The triaxiality, , is the dimensionless parameter that quantifies the degree of permanent nonrotationally symmetric deformation of the planet. It is annihilated for a spherically symmetric body, where A = B = C.
In the vicinity of the spinorbit resonance s∕2, the torque due to triaxiality is dominated by the degrees term in Eq. (76) when averaged over an orbital period (Goldreich 1966; Goldreich & Peale 1966; Ribas et al. 2016) and can thus be approximated by^{2} (77)
which is oscillatory and leads to libration when the planet is in spinorbit resonance. This approximation allows us to rewrite the equation of the spin evolution as (78)
where we have made use of the relationship and denoted the second order time derivative by ¨ . Assuming that the triaxial torque predominates in the vicinity of the resonance, we neglect the tidal torque (). It follows that the maximum absolute value that can reach inside the resonance is given by (Goldreich & Peale 1966, Eq. (17)) (79)
Following Ribas et al. (2016), Δ is called the width of the resonance in the following, since the separatrix between the librating (trapped) and circulating states is defined as (80)
The width of the resonance defines the frequency range where a capture may occur, that is the range where the planet may be driven towards the exact spinorbit ratio of the resonance by the triaxial torque. In this range, the existence of rotation equilibria distinct from the spinorbit resonance depends on the combination of the triaxial and tidal torques. As a consequence, asynchronous equilibria defined by the annihilation of the tidal torque and plotted in Fig. 6 may not exist if .
In addition to this criterion, Goldreich & Peale (1966) derived a probability of capture in the case of weak tidal torques. Their theory, generalised by Makarov (2012), established that this probability only depends on the ratio between the even and odd components of the tidal torque with respect to the s∕2 resonance, defined as
The even part of the torque tends to make the planet traverse the resonance, while the odd part drives it back towards this configuration. Thus the capture probability of the s∕2 spinorbit resonance is defined for any tidal torque as (Goldreich & Peale 1966; Makarov 2012; Ribas et al. 2016) (83)
where the integrals should be performed over the separatrix between the librating and circulating states, given by Eq. (80). We note that when the integral of the odd component is greater than that of the even component, meaning that the capture always occurs in this case.
We focus on the 1:1 spinorbit resonance. In this case, the mechanism described in this section has strong repercussions on the climate and surface conditions of the planet since it determines whether the body is tidally locked in synchronous rotation or in a nonsynchronous state. For comparisons with results obtained in previous sections, we introduce the normalised width of the 1: 1 resonance, (84)
where designates the rotation period at the boundary of the resonance, that is such that Ω −n_{⋆} = Δ.
This normalised width is plotted as a function of the logarithm of the triaxiality of the planet (Fig. 7, left panel), with indicative levels corresponding to the triaxialities of several rocky planets of the Solar system (in decreasing orders of magnitude): ~ 1.4 × 10^{−4} for Mercury, as derived by Ribas et al. (2016) from the gravity moments measured by Smith et al. (2012), ~ 2 × 10^{−5} for the Earth, and ~6 × 10^{−6} for Venus (Yoder 1995). Similarly, the length of the Solar day at the boundary of the resonance () is plotted as a function of the orbital period and triaxiality in logarithmic scales (Fig. 7, middle panel). In these two plots, we assume , since this Hansen coefficient hardly varies for e ≲ 0.3.
The first plot shows that typical values of the normalised width of the 1: 1 spinorbit resonance fall within the interval 0.5–3% of the orbital period, which corresponds to the yellow areas of the parameter space in Fig. 6. The Solar day associated with these values, varies from ~50 to ~ 3000 Earth days in the range P_{⋆} ~ 1−10 Earth days. As a first approximation, such day lengths may be considered as upper estimations of those that can be reached by asynchronous rotation equilibria.
Using this criterion, we calculate the triaxiality necessary to make the width of the synchronisation coincide with the final state determined by the tidal torque, (85)
where Ω_{eq} designates the rotation rate of the final state derived from the tidal torque. This triaxiality is plotted in logarithmic scale in the case of an Earthsized planet of sixday orbital period with σ_{R} = 10^{−6} s^{−1} (Fig. 6, top middle panel) as a function of the logarithms of the eccentricity and ocean depth (Fig. 7, right panel).
The very small triaxialities obtained outside of oceanic tidal resonances for e ≲ 0.25 reveal that the planet is likely to end in tidally locked spinorbit synchronous rotation in this region of the parameter space, instead of being driven towards slightly asynchronous states of equilibrium. Conversely, the high triaxialities observed for asynchronous states induced by the oceanic tide (red areas) suggest that the corresponding final rotation rates are not affected by the triaxial torque unless triaxiality is very high.
In this section, we characterised the evolution of final rotations with the eccentricity and the ocean depth in light of the theory of capture in spinorbit resonance, and recovered the features described by the analytical theory detailed in Sect. 5. We now focus on the critical eccentricity below which the planet is driven towards spinorbit synchronous rotation. We aim to identify the regions of the parameter space where this eccentricity is the smallest.
Fig. 7 Left panel: normalised width of 1:1 spinorbit resonance as a function of planet triaxiality. The blue zone shows the area covered by librations inside the resonance, which is the resonance width. Dotted and dashed brown, green and orange lines correspond to the triaxialities of Mercury (~ 1.4 × 10^{−4}; Smith et al. 2012; Ribas et al. 2016), Earth (~2 × 10^{−5}), and Venus (~6 × 10^{−6}; Yoder 1995), respectively.Middle panel: solar day logarithm of planet at 1:1 resonance boundary as a function of its orbital period in Earth days (horizontal axis) and triaxiality (vertical axis) in logarithmic scales. Right panel: triaxiality logarithm. This is given such that the final asynchronous state of equilibrium of the planet corresponds to the width of the 1:1 spinorbit resonance as a function of the eccentricity (horizontal axis) and ocean depth (vertical axis, km) of the planet in logarithmic scales. The treated case corresponds to the top middle panel of Fig. 6 (M_{p} = M_{⊕}, P_{⋆} = 6 days, and σ_{R} = 10^{−6} s^{−1}). 
7 Evolution of critical eccentricity with orbital period and ocean depth
The critical eccentricity is the parameter that determines whether or not asynchronous rotation states of equilibrium may exist outside of spinorbit resonances (Ω_{eq}:n_{⋆} = 3:2, 2:1, 5:2, …; e.g. Makarov & Efroimsky 2013; Correia et al. 2014). It provides qualitative information on the rotation state at which a planet can be found. Thus, we examine in this section how the critical eccentricity depends on the orbital period of the planet and its ocean depth. Since we obtained similar results for the Earth and superEarth studied in the preceding section, we consider here the Earth case solely.
As highlighted by Fig. 6, oceanic tides can generate asynchronous states in the vicinity of the synchronisation (Ω_{eq} ≈ n_{⋆}). In these states,the planet may be tidally locked in the 1:1 spinorbit resonance by the torque due to triaxiality. Hence, we are interested in the asynchronous states that are outside of the resonance. For this reason, we define the critical eccentricity as the minimal eccentricity such that the finalrotation period satisfies the condition P_{rot}∕P_{⋆} ≤ 0.95. This means that the departure between the finalrotation period and the orbital period of the planet has to be greater than 5% of the orbital period, which is the upper estimation of the width of the 1: 1 resonance obtained for a triaxiality ~10^{−3} (see Fig. 7, left panel).
In Fig. 8, we plot the logarithm of the critical eccentricity calculated using the above definition (top panels) and the associated normalised finalrotation period P_{rot}∕P_{⋆} (bottom panels) as functions of the logarithms of the orbital period (horizontal axis) and ocean depth (vertical axis), for a strong drag (σ_{R} = 10^{−5} s^{−1}, left), and a weak drag (σ_{R} = 10^{−6} s^{−1}, right). To obtain these maps, we performed the calculations detailed in Sect. 6 for each of the sampled orbital periods with the same parameterisation.
We first consider the maps of the critical eccentricity (top panels). In these maps, yellowgreen areas designate the nonresonant regime, where the critical eccentricity is high (log_{10}e_{AR} ≈−0.6). This corresponds to the asymptotic configuration of a dry rocky planet. Resonances of the oceanic tidal response induce blue diagonal bands, where the critical eccentricity is decreased with respect to the reference value. Each of these bands, from left to right, is associated with an eccentricity degree in ascending order, the main one being due to the degree3 eccentricity term. The diagonal pattern follows the scaling law derived analytically for the resonant ocean depths (see Eq. (60)), that is (86)
Moreover, the contrast between the nonresonant and resonant regimes becomes stronger as the Rayleigh drag frequency and orbital period decay, in agreement with the scaling laws given by Eqs. (67) and (72), that is (87)
respectively. It should be remembered that these scaling laws were derived in the framework of approximations such as the quasiadiabatic () and nonrotating (no Coriolis effects) approximations, which leads to differences with the results of numerical calculations, and particularly in the strong drag configuration (σ_{R} = 10^{−5} s^{−1}).
We now move to the maps representing the finalrotation rates associated with the computed critical eccentricity (bottom panels). Yellow areas designate states of equilibrium close to synchronisation, and blue areas the 3: 2 spinorbit resonance. In order to avoid confusion when looking at these maps, one should bear in mind that they show both the loweccentricity regime, which is the subject of the analytical theory detailed in Sects. 4 and 5, and the higheccentricity regime, where the derived results do not apply since it is beyond the scope of this work.
We shall also stress here a somehow counterintuitive feature of the behaviour of the final spin rate that we already emphasised in Fig. 4 and discussed in Sect. 5: as the ocean depth increases, the peak generated by a resonant eccentricity term moves towards the spinorbit synchronisation, and so does the associated asynchronous state of equilibrium. Thus, the asynchronous state generated by a strong eccentricity term is not necessarily far from the synchronisation. We retrieve this feature in Fig. 7 (bottom panels). While no correlation can be observed between the final rotation and the critical eccentricity in the loweccentricity regime, the displacement of the state of equilibrium induces a colour gradation from blue to yellow along the direction of ascending ocean depths.
Up to now, we examined the final rotation state of equilibrium of the planet without considering how much time is necessary to reach this state. The evolution timescale of the planet rotation rate is related to the tidally dissipated energy, which strongly depends on the forcing frequencies associated with tidal components (e.g. AuclairDesrotour et al. 2014). It is thus important to complete the present work with a study of the evolution of the planet spin, which is the objective of the next section.
Fig. 8 Logarithm of critical eccentricity e_{AR} (top panels) and associated normalised finalrotation period P_{rot}∕P_{⋆} (bottom panels) as functions of logarithms of orbital period (horizontal axis) and ocean depth (vertical axis), for σ_{R} = 10^{−5} s^{−1} (left) and σ_{R} = 10^{−6} s^{−1} (right). The calculated critical eccentricity is such that the finalrotation period satisfies the condition P_{rot} ∕P_{⋆} ≤ 0.95 to retain solely asynchronous states clearly separated from the synchronisation. In the eccentricity maps (top panels), the yellow colour designates the region of the nonresonant regime, where the critical eccentricity is high (log_{10} e_{AR} ≈−0.6). Conversely, the blue shades designate the regions where e_{AR} is low owing to the resonances of the oceanic tidal response. In the finalrotation maps, the yellow colour designates states that are close to synchronisation, while the dark blue colour corresponds to the 3/2 spinorbit resonance, induced by the degree3 eccentricity term. Values of parameters used for the solid part are given by Table 1. In both cases, the stellar mass is M_{⋆} = 0.09 M_{⊙}, the initial rotation rate is , and 10 Hough modes are taken into account in the calculation of the tidal response. 
8 Evolution timescale
We return to the case of the Earthsized planet orbiting the TRAPPIST1 star with a sixday orbital period. The features of the planet and of its solid part in particular are determined by the parameters given by Table 1. As done in Sect. 6, we consider the strong and weak drag configurations characterised by σ_{R} = 10^{−5} s^{−1} and σ_{R} = 10^{−6} s^{−1}, respectively. Three orbits of various eccentricity are studied. In the first one, e = 0.063, which is the value of the critical eccentricity predicted by the theory in the weak drag configuration (see Fig. 6, top middle panel). In the second orbit, e is set to 0.2,that is just below the value of the critical eccentricity obtained in the absence of ocean. In the third orbit, e is set to 0.3, which is slightly greater than this critical eccentricity.
In each configuration we compute the evolution of the spin rotation rate of a set of planets characterised by ocean depths sampled in the interval − 2 ≤ γ ≤ 0, where with H_{oc} given in kilometres. As mentioned in Sect. 6, the range of values H_{oc} ≲ 0.1 km does not correspond to realistic cases since the global ocean approximation is no longer valid for such small depths. The chosen lower boundary γ = −2 is intended to emphasise the continuous transition between the asymptotic regime of dry planets (H_{oc} → 0) in our ocean planet tidal model and the case of a dry solid body described by the Andrade model.
The initial value of the spin rotation rate is set to , which corresponds to the middle of the interval of forcing frequencies defined by the 1: 1 and 3: 2 spinorbit resonances. If we ignore the librating motion of the planet and the associated triaxial torque, the evolution of the spin angular velocity is determined by the equation (88)
which is integrated over time using a variable step size BulirschStoer method (Press et al. 2007). The normalised moment of inertia is set to Č = 0.3308 (Bullen 2012) and the tidal torque exerted on the planet is calculated using the expression given by Eq. (44).
We plot in Fig. 9 the normalised rotation period P_{rot}∕P_{⋆} (vertical axis) as a function of time in logarithmic scale (horizontal axis) for each of the above defined configurations: the strong (top panels) and weak (bottom panels) drag regimes with eccentricities e = 0.063, 0.2, and 0.3 (from left to right). Colours of solid lines varies from orange to blue as the ocean depth increases.
The order of magnitude of the evolution timescale τ_{ev} can be estimated by considering the contribution of the solid part in the absence of ocean in the circular case. By assuming that σ_{2,2} ≈ n_{⋆}∕2 = π∕P_{⋆}, we thus obtain (89)
where the expression of as a function of the Andrade parameters is given by Eq. (47). The substitution of parameters by their numerical values leads to τ_{ev} ~ 10^{3} yr, and we hence recover the order of magnitude given by numerical calculations.
The plots of Fig. 9 show the three possible regimes of the evolution of the planetary spin. In the first regime (e = 0.063, left panels), the planet is driven towards the spinorbit synchronous rotation in most of cases. Particularly, in the strong drag configuration, all of the planets are locked in this state of equilibrium after ~ 5000 yr whatever the depth of their oceanic layer. For σ_{R} = 10^{−6} s^{−1}, two of the studied cases avoid the synchronisation owing to the action of oceanic tides and evolve towards asynchronous rotation states located between the 1:1 and 3:2 spinorbit resonances.
The second regime (e = 0.2, middle panels) corresponds to the transition between low and high eccentricities. In this regime, a solid planet is driven towards the synchronisation while most of those hosting an ocean converge towards the nonsynchronised rotation states of equilibrium identified in our exploration of the parameter space (see Fig. 6, top left and middle panels). We note that trajectories of planets are affected by the variations of the tidal torque with the spin angular velocity. A resonant peak associated with an oceanic tidal mode generates a rapid evolution, which slows down when the tidal torques returns to its nonresonant level.
In the third regime (e = 0.3, right panels) the contribution of nonresonant eccentricity terms is so strong that it counterbalances the action of the semidiurnal component. As a consequence, intermediate rotation sates of equilibrium cannot exist and all of the planets converge towards the 3: 2 spinorbit resonances where they are locked by solid tides. The ocean just slightly modifies the time necessary to reach the state of equilibrium in this case.
9 Conclusions
By combining the linear theory of oceanic tides with realistic models of the solid tide, we derived a selfconsistent analytic model to characterise the response of an ocean planet undergoing eccentricity tides. Following early works (e.g. Tyler 2011, 2014; Chen et al. 2014; Matsuyama 2014; Matsuyama et al. 2018), this model assumes the shallow water approximation and includes the effects of ocean loading, selfattraction, and deformation of the solid regions, which tightly couple the oceanic tidal response with that of the solid part. The energy dissipation induced by the interactions of tidal flows with the planet topography is taken into account with a Rayleigh drag, which enabled us to extend the usual theoretical treatments applied in adiabatic regimes to the dissipative case.
The tidal response of the solid part is described in our approach by the Andrade model, with parameters computed from advanced models of the internal structure and tidal oscillations of solid bodies. The Andrade model solves the inability of commonly used models – such as the Maxwell model for instance – to properly describe how the tidally dissipated energy scales with the tidal frequency, which usually leads to underestimate it by several orders of magnitude in the highfrequency regime.
Firstly, we derived an analytic solution for the tidal torque exerted on the planet, which allowed us to characterise the main features of this torque in the quasiadiabatic regime. Particularly, we identified the dependence of eigenfrequencies associated with oceanic modes on the planet parameters and quantified the maximum of resonant peaks. This led us to derive an analytic estimation of the critical eccentricity beyond which asynchronous rotation states of equilibrium can be induced by eccentricity tides.
Secondly, we used the obtained expression of the tidal torque to numerically explore the parameter space in the case of an Earth and a superEarth. We thus showed that resonances of oceanic modes are likely to decrease the critical eccentricity by one order of magnitude, enabling thereby the existence of asynchronous rotation states of equilibrium distinct from the spinorbit resonances induced by tides in the case of solid bodies (Makarov & Efroimsky 2013; Correia et al. 2014). Particularly, for typical configurations, the critical eccentricity can switch from ~ 0.3 to ~ 0.06 owing to the resonant amplification of the degree3 (or eastward propagating) eccentricity term, which is the largest one in the loweccentricity regime. Calculations allowed us to determine the region of the parameter space where the impact of resonances on the critical eccentricity is the most significant, and to establish scaling laws characterising its dependence on the system features (orbital period of the planet, ocean depth, and Rayleigh drag frequency).
Thirdly, we highlighted the action of resonances associated with oceanic modes on the history of the planetary rotation by coupling the time evolution equation of the spin with our tidal model. Results shows the three possible regimes of evolution: (i) the loweccentricity limit (e ≪ 0.2); where planets are driven towards synchronisation except in case of strong resonance of the oceanic tidal response; (ii) the transition regime (e ≈ 0.2), where a large fraction of ocean planets converge towards asynchronous states while solid ones are driven towards the spinorbit synchronous rotation; and (iii) the higheccentricity regime (e ≳ 0.3), where the eccentricity is strong enough to make planets converge towards the 3: 2 spinorbit resonance induced by solid tides whatever the depth of their oceans.
These results can be used to better constrain the rotation of discovered rocky planets. In the case of TRAPPIST1 planets, eccentricities should be limited to low values (Luger et al. 2017; Grimm et al. 2018; Turbet et al. 2018). Planets b and c are in a runaway greenhouse state so that they are too hot to sustain surface liquid water. Planet d, depending on the habitable zone model, could be able to sustain surface liquid water or it could be in a runaway greenhouse state. However, being relatively closein, planet d could experience a strong tidal heating (e.g. Makarov et al. 2018), which could trigger a tidal runaway greenhouse (Barnes et al. 2013) and prevent the presence of a global surface ocean. Planet f to h are too cold to sustain surface liquid water but could still host interior oceans below a layer of ice (Turbet et al. 2018). Our formalism does not allow us to model such planets at this point.
However, planet e has the highest probability of hosting a global ocean of liquid water. It could have retained its initial water content (Bourrier et al. 2017), and if there is water, it could be liquid (Turbet et al. 2018). And due to the small eccentricity of its orbit, it is likely to be tidally locked in synchronous rotation in our framework, unless the ocean has a very specific depth and tidal flows undergo a weak drag. Note that the strong tidal heating obtained by Makarov et al. (2018) could trigger a runaway greenhouse state for planet e as well. However they provide an upper estimate of the tidal heating and other studies have found much lower tidal heat fluxes for planet e (Barr et al. 2018; Bolmont et al., in prep.), which are compatible with surface liquid water.
In the case of Proxima, the eccentricity could be as high as 0.1 (AngladaEscudé et al. 2016; Ribas et al. 2016), meaning that the planet could be locked in a spin state different from the 1: 1 and 3: 2 spinorbit resonances, depending on the depth of a supposed global ocean. Thus the synchronous rotation state is less probable in this case than in that of TRAPPIST1 planets.
Our approach based on the linear theory of tides can help to better understand the action of eccentricity tides on the rotation of ocean planets, and provides solutions that can be implemented in evolutionary codes to compute the tidal Love numbers and torque in a realistic and nevertheless efficient way. Particularly, the linear analysis leads to a full characterisation of the frequency behaviour of these parameters, which is crucial in the case of oceanic tides. However, it presents limitations that make it complementary with general circulation models.
The main limitation of our model is that it cannot be applied to large tidal distortions since these imply nonlinear mechanisms. Thus, resonances violating the small perturbation approximation are beyond the scope of the linear theory and shall be studied by integrating the full primitive equations of the ocean (e.g. Vallis 2006). This limitation prevented us from quantifying the critical eccentricity in the quasiadiabatic regime, below σ_{R} ~ 10^{−6} s^{−1}. A second limitation is related to the shallow water approximation, which eludes the role played by the oceaninternal structure on the tidally dissipated energy. Although the formalism can easily be extended to stably stratified deep oceans if the effects of ocean loading, selfattraction, and deformation of the solid regions are ignored (AuclairDesrotour et al. 2018), taking these ingredients into accountcomplicates the analytic treatment of the problem significantly.
Finally, we would like to draw the attention of the reader to the importance of the Rayleigh drag frequency in this approach. The Rayleigh drag frequency accounts for the global effect of topography and ocean streams, which converts barotropic tidal flows into internal gravity waves. For Earth, this parameter has been estimated from the energy dissipated by the semidiurnal Lunar tide (e.g. Webb 1980) and satellite measurements of the ocean surface elevation (e.g. Egbert & Ray 2001, 2003). Given the impact it has on the planet tidal response, it should be characterised as a function of the planet properties (spin rotation, topography, ocean depth) in order to improve the predictions of the linear theory from a quantitative point of view. Moreover, in future studies one should investigate how the presence of a partial or global ice cap may affect the results obtained here for a freesurface ocean planet, following the methodology adopted for icy satellites (e.g. Kamata et al. 2015; Matsuyama et al. 2018).
Fig. 9 Evolution of normalised rotation period P_{rot}∕P_{⋆} as function of time (yr) in logarithmic scale. This is given in the case of an Earthsized planet orbiting the TRAPPIST1 star with a sixday orbital period. The Rayleigh drag frequency of the ocean is set to σ_{R} = 10^{−5} s^{−1} for cases of strong drag (top panels), and to σ_{R} = 10^{−6} s^{−1} for cases of weak drag (bottom panels). The planet eccentricity takes the values e = 0.063, 0.2, and 0.3 from left to right. In each panel, a set of planets with various ocean depths is studied. These configurations are parametrised by , where the ocean depth is given in kilometers. The parameter γ is sampled from −2 (the driest case of the sample, solid orange line) to 0 (the most humid case, solid blue line). Values used for the planet parameters are given by Table 1. 
Acknowledgements
The authors thank the anonymous referee for constructive comments that improved the manuscript. This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreements No. 679 030/WHIPLASH and 771 620/EXOKLEIN), and has been carried out within the framework of the NCCR PlanetS supported by the Swiss National Science Foundation. S.M. acknowledges funding by the European Research Council through the ERC grant No. 647 383/SPIRE and by the PLATO CNES grant at CEA Saclay. E.B. and P.A.D. are grateful to O. Grasset and G. Tobie, whose models were used in this work to calculate the internal structure of the studied rocky planets and the tidal viscoelastic response of their solid part, respectively. This research has made use of NASA’s Astrophysics Data System.
Appendix A Associated Legendre functions
Associated Legendre functions are defined for and as (Abramowitz & Stegun 1972; Arfken & Weber 2005) (A.1)
the designating the Legendre polynomials, given by (A.2)
The normalised associated Legendre functions are characterised by (A.3)
where the Kronecker symbol δ_{l,k} is, such that δ_{l,k} = 1 if l = k, and otherwise δ_{l,k} = 0. Hence they are expressed as (A.4)
Appendix B Hansen coefficients
Hansen coefficients are derived from the Fourier coefficients of the eccentricity (e) and mean anomaly (M) function (e.g. Hughes 1981; Polfliet & Smeyers 1990; Laskar 2005) (B.1)
where v designate the true anomaly and r_{⋆} the star planet distance, these two quantities being themselves expressed as
Hansen coefficients are thus defined by (e.g. Hughes 1981) (B.4)
We note that the are real since is an even function of M.
By changing the sign of s and using the symmetry property , we obtain (B.5)
which is convenient to put the perturbing tidal potential into the form given by Eq. (3).
In practice, considering the fact that Hansen coefficients are the Fourier coefficients of the function , one can efficiently get a numerical estimate of 2^{N} + 1 coefficients for − K ≤ s ≤ K with K = 2^{N−1} by means ofa fast Fourier transform (FFT; Adel Sharaf & Hassan Selim 2010; Correia et al. 2014). The integer N shall be chosen appropriately so that , and increases with the eccentricity as the spectrum of forcing terms widens (see Fig. 3 in Ogilvie 2014).
Appendix C Torque due to triaxiality
In this appendix, we detail the derivation of the expression of the triaxial torque in the general case, given by Eq. (76). This torque is expressed as a function of the planet spin angle (ϑ) and true anomaly (v) as (e.g. Goldreich & Peale 1966, Eq. (2)) (C.1)
which can be rewritten as , where (C.2)
By using Eq. (B.5) and the notation , we thus obtain (C.3)
Appendix D Nomenclature
Notations used along this work and their designations in order of appearance in the text.
References
 Abramowitz, M., & Stegun, I. A. 1972, Handbook of Mathematical Functions (Mineola, New York: Dover Publications) [Google Scholar]
 Adel Sharaf, M., & Hassan Selim, H. 2010, Res. Astron. Astrophys., 10, 1298 [NASA ADS] [CrossRef] [Google Scholar]
 Andrade, E. N. D. C. 1910, Proc. Royal Soc. Lond. Ser. A, 84, 1 [Google Scholar]
 Andrade, E. N. D. C. 1914, Proc. Royal Soc. Lond. Ser. A, 90, 329 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 AngladaEscudé, G., Amado, P. J., Barnes, J., et al. 2016, Nature, 536, 437 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Arfken, G. B.,& Weber, H. J. 2005, Mathematical Methods for Physicists, 6th edn. (Cambridge: Academic Press) [Google Scholar]
 AuclairDesrotour, P., Le PoncinLafitte, C., & Mathis, S. 2014, A&A, 561, L7 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Auclair Desrotour, P., Mathis, S., & Le PoncinLafitte, C. 2015, A&A, 581, A118 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 AuclairDesrotour, P., Laskar, J., & Mathis, S. 2017a, A&A, 603, A107 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 AuclairDesrotour, P., Laskar, J., Mathis, S., & Correia, A. C. M. 2017b, A&A, 603, A108 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 AuclairDesrotour, P., Mathis, S., Laskar, J., & Leconte, J. 2018, A&A, 615, A23 [EDP Sciences] [Google Scholar]
 AuclairDesrotour, P., Leconte, J., & Mergny, C. 2019, A&A, 624, A17 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Barnes, R., Mullins, K., Goldblatt, C., et al. 2013, Astrobiology, 13, 225 [NASA ADS] [CrossRef] [Google Scholar]
 Barr, A. C., Dobos, V., & Kiss, L. L. 2018, A&A, 613, A37 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Beuthe, M. 2016, Icarus, 280, 278 [CrossRef] [Google Scholar]
 Bolmont, E., Selsis, F., Owen, J. E., et al. 2017, MNRAS, 464, 3728 [Google Scholar]
 Bourrier, V., de Wit, J., Bolmont, E., et al. 2017, AJ, 154, 121 [NASA ADS] [CrossRef] [Google Scholar]
 Breton, S., Bolmont, E., Tobie, G., Mathis, S., & Grasset, O. 2018, in SF2A2018: Proceedings of the Annual meeting of the French Society of Astronomy and Astrophysics, 195 [Google Scholar]
 Bullen, K. E. 2012, The Earth’s Density (Berlin: Springer Science & Business Media) [Google Scholar]
 Castelnau, O., Duval, P., Montagnat, M., & Brenner, R. 2008, J. Geophys. Res. Solid Earth, 113, B11203 [NASA ADS] [CrossRef] [Google Scholar]
 CastilloRogez, J. C., Efroimsky, M., & Lainey, V. 2011, J. Geophys. Res. Planets, 116, E09008 [NASA ADS] [CrossRef] [Google Scholar]
 Chen, E. M. A., Nimmo, F., & Glatzmaier, G. A. 2014, Icarus, 229, 11 [NASA ADS] [CrossRef] [Google Scholar]
 Correia, A. C. M., & Laskar, J. 2001, Nature, 411, 767 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Correia, A. C. M., & Laskar, J. 2003, J. Geophys. Res. Planets, 108, 5123 [NASA ADS] [CrossRef] [Google Scholar]
 Correia, A. C. M., & Laskar, J. 2009, Icarus, 201, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Correia, A. C. M., Boué, G., Laskar, J., & Rodríguez, A. 2014, A&A, 571, A50 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Cottrell, A. H., & Aytekin, V. 1947, Nature, 160, 328 [CrossRef] [Google Scholar]
 Cowling, T. G. 1941, MNRAS, 101, 367 [NASA ADS] [Google Scholar]
 Dermott, S. F. 1979, Icarus, 37, 310 [NASA ADS] [CrossRef] [Google Scholar]
 Dobos, V., Barr, A. C., & Kiss, L. L. 2019, A&A, 624, A2 [Google Scholar]
 Dobrovolskis, A. R., & Ingersoll, A. P. 1980, Icarus, 41, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Duval, P. 1978, J. Glaciol., 21, 621 [NASA ADS] [CrossRef] [Google Scholar]
 Efroimsky, M. 2012, ApJ, 746, 150 [NASA ADS] [CrossRef] [Google Scholar]
 Efroimsky, M., & Lainey, V. 2007, J. Geophys. Res. Planets, 112, E12003 [NASA ADS] [CrossRef] [Google Scholar]
 Efroimsky, M., & Williams, J. G. 2009, Celest. Mech. Dyn. Astron., 104, 257 [Google Scholar]
 Egbert, G. D., & Ray, R. D. 2001, J. Geophys. Res., 106, 22 [Google Scholar]
 Egbert, G. D., & Ray, R. D. 2003, Geophys. Res. Lett., 30, 1907 [NASA ADS] [CrossRef] [Google Scholar]
 Gastineau, M., & Laskar, J. 2011, ACM Commun. Comput. Algebra, 44, 194 [Google Scholar]
 Gillon, M., Triaud, A. H. M. J., Demory, B.O., et al. 2017, Nature, 542, 456 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Gold, T., & Soter, S. 1969, Icarus, 11, 356 [NASA ADS] [CrossRef] [Google Scholar]
 Goldreich, P. 1966, AJ, 71, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Goldreich, P., & Peale, S. 1966, AJ, 71, 425 [NASA ADS] [CrossRef] [Google Scholar]
 Greenberg, R. 2009, ApJ, 698, L42 [Google Scholar]
 Grimm, S. L., Demory, B.O., Gillon, M., et al. 2018, A&A, 613, A68 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Henning, W. G., & Hurford, T. 2014, ApJ, 789, 30 [NASA ADS] [CrossRef] [Google Scholar]
 Henning, W. G., O’Connell, R. J., & Sasselov, D. D. 2009, ApJ, 707, 1000 [NASA ADS] [CrossRef] [Google Scholar]
 Hough, S. S. 1898, Roy. Soc. London Philos. Trans. Ser. A, 191, 139 [Google Scholar]
 Hughes, S. 1981, Celest. Mech., 25, 101 [NASA ADS] [CrossRef] [Google Scholar]
 Ingersoll, A. P., & Dobrovolskis, A. R. 1978, Nature, 275, 37 [NASA ADS] [CrossRef] [Google Scholar]
 Jackson, I. 1993, Geophys. Res. Lett., 20, 2115 [NASA ADS] [CrossRef] [Google Scholar]
 Jackson, I., Fitz Gerald, J. D., Faul, U. H., & Tan, B. H. 2002, J. Geophys. Res. Solid Earth, 107, 2360 [Google Scholar]
 Kamata, S., Matsuyama, I., & Nimmo, F. 2015, J. Geophys. Res. Planets, 120, 1528 [NASA ADS] [CrossRef] [Google Scholar]
 Karato, S., & Spetzler, H. A. 1990, Rev. Geophys., 28, 399 [CrossRef] [Google Scholar]
 Kaula, W. M. 1966, Theory of Satellite Geodesy. Applications of Satellites to Geodesy (North Chelmsford: Courier Corporation) [Google Scholar]
 Laskar, J. 2005, Celest. Mech. Dyn. Astron., 91, 351 [NASA ADS] [CrossRef] [Google Scholar]
 Leconte, J., Wu, H., Menou, K., & Murray, N. 2015, Science, 347, 632 [NASA ADS] [CrossRef] [Google Scholar]
 Lee, U., & Saio, H. 1997, ApJ, 491, 839 [NASA ADS] [CrossRef] [Google Scholar]
 Lindzen, R. S., & Chapman, S. 1969, Space Sci. Rev., 10, 3 [NASA ADS] [CrossRef] [Google Scholar]
 LonguetHiggins, M. S. 1968, Phil. Trans. R. Soc. London, Ser. A, 262, 511 [Google Scholar]
 Love, A. E. H. 1911, Some Problems of Geodynamics (Cambridge: Cambridge University Press) [Google Scholar]
 Luger, R., Sestovic, M., Kruse, E., et al. 2017, Nat. Astron., 1, 0129 [Google Scholar]
 Makarov, V. V. 2012, ApJ, 752, 73 [NASA ADS] [CrossRef] [Google Scholar]
 Makarov, V. V., & Efroimsky, M. 2013, ApJ, 764, 27 [NASA ADS] [CrossRef] [Google Scholar]
 Makarov, V. V., Berghea, C. T., & Efroimsky, M. 2018, ApJ, 857, 142 [NASA ADS] [CrossRef] [Google Scholar]
 Mathis, S., & Le PoncinLafitte, C. 2009, A&A, 497, 889 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Matsuyama, I. 2014, Icarus, 242, 11 [NASA ADS] [CrossRef] [Google Scholar]
 Matsuyama, I., Beuthe, M., Hay, H. C. F. C., Nimmo, F., & Kamata, S. 2018, Icarus, 312, 208 [NASA ADS] [CrossRef] [Google Scholar]
 Munk, W. H.,& MacDonald, G. J. F. 1960, J. Geophys. Res., 65, 2169 [NASA ADS] [CrossRef] [Google Scholar]
 Noyelles,B., Frouard, J., Makarov, V. V., & Efroimsky, M. 2014, Icarus, 241, 26 [NASA ADS] [CrossRef] [Google Scholar]
 Ogilvie, G. I. 2014, ARA&A, 52, 171 [Google Scholar]
 Ogilvie, G. I., & Lin, D. N. C. 2004, ApJ, 610, 477 [NASA ADS] [CrossRef] [Google Scholar]
 Papaloizou, J. C. B., Szuszkiewicz, E., & Terquem, C. 2018, MNRAS, 476, 5032 [NASA ADS] [CrossRef] [Google Scholar]
 Peale, S. J., & Boss, A. P. 1977, J. Geophys. Res., 82, 3423 [NASA ADS] [CrossRef] [Google Scholar]
 Polfliet, R., & Smeyers, P. 1990, A&A, 237, 110 [Google Scholar]
 Press, W. H., Teukolsky, S. A., Vetterling, W. T., & Flannery, B. P. 2007, Numerical Recipes: The Art of Scientific Computing, 3rd edn. (Cambridge: Cambridge University Press) [Google Scholar]
 Remus, F., Mathis, S., Zahn, J.P., & Lainey, V. 2012, A&A, 541, A165 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Renaud, J. P., & Henning, W. G. 2018, ApJ, 857, 98 [NASA ADS] [CrossRef] [Google Scholar]
 Ribas, I., Bolmont, E., Selsis, F., et al. 2016, A&A, 596, A111 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Sabadini, R., & Vermeersen, B. 2004, Global Dynamics of the Earth: Applications of Normal Mode Relaxation Theory to SolidEarth Geophysics (Dordrecht, The Netherlands: Kluwer Academic Publishers) [Google Scholar]
 Smith, D. E., Zuber, M. T., Phillips, R. J., et al. 2012, Science, 336, 214 [NASA ADS] [CrossRef] [Google Scholar]
 Sotin, C., Grasset, O., & Mocquet, A. 2007, Icarus, 191, 337 [NASA ADS] [CrossRef] [Google Scholar]
 Storch, N. I., & Lai, D. 2014, MNRAS, 438, 1526 [CrossRef] [Google Scholar]
 Takeuchi, H., & Saito, M. 1972, Methods in Computational Physics (Cambridge: Academic Press), 11, 217 [Google Scholar]
 Tan, B. H., Jackson, I., & Fitz Gerald, J. D. 2001, Phys. Chem. Miner., 28, 641 [Google Scholar]
 Taylor, G. I. 1936, Proc. Royal Soc. Lond. Ser. A, 156, 318 [Google Scholar]
 Tobie, G., Mocquet, A., & Sotin, C. 2005, Icarus, 177, 534 [NASA ADS] [CrossRef] [Google Scholar]
 Townsend, R. H. D. 2003, MNRAS, 340, 1020 [NASA ADS] [CrossRef] [Google Scholar]
 Turbet, M., Bolmont, E., Leconte, J., et al. 2018, A&A, 612, A86 [CrossRef] [EDP Sciences] [Google Scholar]
 Tyler, R. H. 2008, Nature, 456, 770 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Tyler, R. H. 2009, Geophys. Res. Lett., 36, L15205 [Google Scholar]
 Tyler, R. 2011, Icarus, 211, 770 [NASA ADS] [CrossRef] [Google Scholar]
 Tyler, R. 2014, Icarus, 243, 358 [NASA ADS] [CrossRef] [Google Scholar]
 Unno, W., Osaki, Y., Ando, H., Saio, H., & Shibahashi, H. 1989, Nonradial Oscillations of Stars (Tokyo: University of Tokyo Press) [Google Scholar]
 Unterborn, C. T., Hinkel, N. R., & Desch, S. J. 2018, Res. Notes Amer. Astron. Soc., 2, 116 [NASA ADS] [CrossRef] [Google Scholar]
 Vallis, G. K. 2006, Atmospheric and Oceanic Fluid Dynamics (Cambridge: Cambridge University Press), 770 [Google Scholar]
 Van Grootel, V., Fernandes, C. S., Gillon, M., et al. 2018, ApJ, 853, 30 [NASA ADS] [CrossRef] [Google Scholar]
 Volland, H. 1974a, J. Atm. Terr. Phys., 36, 445 [NASA ADS] [CrossRef] [Google Scholar]
 Volland, H. 1974b, J. Atm. Terr. Phys., 36, 1975 [NASA ADS] [CrossRef] [Google Scholar]
 Wang, H., Boyd, J. P., & Akmaev, R. A. 2016, Geosci. Model Dev., 9, 1477 [NASA ADS] [CrossRef] [Google Scholar]
 Webb, D. J. 1980, Geophys. J., 61, 573 [NASA ADS] [CrossRef] [Google Scholar]
 Wunsch, C. 1975, Rev. Geophys. Space Phys., 13, 167 [NASA ADS] [CrossRef] [Google Scholar]
 Yoder, C. F. 1995, Icarus, 117, 250 [NASA ADS] [CrossRef] [Google Scholar]
 Zahn, J. P. 1966, Ann. Astrophys., 29, 313 [NASA ADS] [Google Scholar]
 Zahn, J.P. 1975, A&A, 41, 329 [Google Scholar]
In the case of a thick oceanic layer, the ocean depth intervenes in the expressions of tilt factors as well as in those of the solid Love numbers (Dermott 1979; Remus et al. 2012).
In the vicinity of the spinorbit resonance s∕2, the angle γ_{2,s} is small. Therefore, integrating Eq. (76) over an orbital period results in all eccentricity terms vanishing except the degrees term, which reduces Eq. (75) to the equation of motion of a simple pendulum in the absence of tides (see e.g. Eq. (5) in Goldreich 1966). However, we note that vanishes at the resonance (ϑ∕M = s∕2).
All Tables
Parameters of Andrade model derived from elastogravitational theory (Tobie et al. 2005; Breton et al. 2018).
Notations used along this work and their designations in order of appearance in the text.
All Figures
Fig. 1 Tidal distortion of terrestrial planet resulting from gravitational forcing of host star. The planet structure is composed of a solid part of radius R_{c} and a thin uniform ocean of external radius R_{p} rotating as a solid body at the spin rotation rate Ω. The tidal responses of the two layers are coupled together and characterised by different frequency behaviours. The corresponding angular lags are designated by δ_{c} and δ_{oc}, and the surface displacements by ξ_{c} and ξ_{oc}. The ocean depth modified by solid and oceanic tides is denoted by h. Reference frames introduced in Sect. 2 are drawn in purple. 

In the text 
Fig. 2 Desynchronising mechanism of eccentricity tides. Tidal bulge raised on a planet by the host star in the heliocentric reference frame (left) and in a planetocentric reference frame rotating with the planet mean motion (right). The tidal gravitational forcing is assumed to be dominated by its eccentricity components, given by Eqs. (9) and (10). Because of the angular lag generated by dissipative processes, the rotation of the planet is accelerated by the component of degree s = 3 (the eastward propagating potential) in the vicinity of the periastron. Conversely, it is decelerated by the component of degree s = 1 (the westward propagating potential) to a lesser extent in the vicinity of the apoapsis. 

In the text 
Fig. 3 Andrade (solid black line) and Maxwell (dotted blue line) models in frequency domain.The negative imaginary part of the complex quadrupole Love number (see Eq. (15)) characterising a solid planet is plotted as a function of the normalised tidal frequency στ_{M} in logarithmic scales. We used the values of the Andrade parameters given by Table 1 for the Earth. 

In the text 
Fig. 4 Positive tidal torque generated by oceanic eccentricity tides in supersynchronous frequency range (Ω > n_{⋆}) for increasing oceanic depths (a–f). An oceanic mode is characterised by a resonance of eigenfrequency σ_{n} (blue peaks). We focus here on the degree0 mode, of frequency σ_{0}, in the nonrotating approximation (Coriolis effects are ignored). An eccentricity term of frequency σ_{2,s} = 2Ω − sn_{⋆} generates a positive torque, which is enhanced by the resonance when (red peak). The semidiurnal (σ_{2,2}) and first eccentricity frequencies (σ_{2,3} and σ_{2,4}) associated with the peak thus created are symbolised by the superscript ▲. The peak raised by the eccentricity tide occurs when , which also corresponds to and . The direction of peaks indicates the sign of the tidal torque. In the supersynchronous regime (σ_{2,2} > 0), downwardblue peaks tend to drive the planet towards tidal locking in spinorbit synchronous (σ_{2,2} = 0), while upward red peaks tend to drive it away from this state of equilibrium. 

In the text 
Fig. 5 Logarithm of tidal torque exerted on planet. In the cases of the pure solid (top panel) and oceanic (middle panel) tidal responses, and in the general case (bottom panel), these are shown as functions of the normalised semidiurnal frequency for a circularorbit (red line) and an orbit of eccentricity e = 0.2. The torques , , and are computed using the expression given by Eq. (44), and the values given by Table 1 for an Earthsized planet orbiting TRAPPIST1 (M_{⋆} = 0.09 M_{⊙}) with a sixday orbital period. The Rayleigh drag frequency is set to σ_{R} = 1.0 × 10^{−6} s^{−1}. Solid (dashed) lines designate the regimes where the tidal torque drives the planet towards (away from) spinorbit synchronousrotation. 

In the text 
Fig. 6 Normalised finalrotation period P_{rot}∕P_{⋆} as a function of eccentricity (horizontal axis) and ocean depth (vertical axis) in logarithmic scales. Two planets are considered: an Earthsized planet with sixday (top panels) and oneday (middle panels) orbital periods, and a 10 M_{⊕} superEarth with a sixday orbital period (bottom panels). From left to right: configurations associated with the values of the Rayleigh drag frequency σ_{R} = 10^{−5}, 10^{−6}, 10^{−7} s^{−1}. Values of parameters used for the solid part are given by Table 1. In all cases, the stellar mass is M_{⋆} = 0.09 M_{⊙}, the initial rotation rate is (middle of the interval 0 < σ_{2,2} < n_{⋆}, see Fig. 4), and 10 Hough modes are taken into account in the calculation of the tidal response. The degrees of eccentricity modes causing resonances are indicated in the case of the sixdayperiod Earthsized planet with a weak drag(top right panel), consistent with the discussion on Fig. 4. 

In the text 
Fig. 7 Left panel: normalised width of 1:1 spinorbit resonance as a function of planet triaxiality. The blue zone shows the area covered by librations inside the resonance, which is the resonance width. Dotted and dashed brown, green and orange lines correspond to the triaxialities of Mercury (~ 1.4 × 10^{−4}; Smith et al. 2012; Ribas et al. 2016), Earth (~2 × 10^{−5}), and Venus (~6 × 10^{−6}; Yoder 1995), respectively.Middle panel: solar day logarithm of planet at 1:1 resonance boundary as a function of its orbital period in Earth days (horizontal axis) and triaxiality (vertical axis) in logarithmic scales. Right panel: triaxiality logarithm. This is given such that the final asynchronous state of equilibrium of the planet corresponds to the width of the 1:1 spinorbit resonance as a function of the eccentricity (horizontal axis) and ocean depth (vertical axis, km) of the planet in logarithmic scales. The treated case corresponds to the top middle panel of Fig. 6 (M_{p} = M_{⊕}, P_{⋆} = 6 days, and σ_{R} = 10^{−6} s^{−1}). 

In the text 
Fig. 8 Logarithm of critical eccentricity e_{AR} (top panels) and associated normalised finalrotation period P_{rot}∕P_{⋆} (bottom panels) as functions of logarithms of orbital period (horizontal axis) and ocean depth (vertical axis), for σ_{R} = 10^{−5} s^{−1} (left) and σ_{R} = 10^{−6} s^{−1} (right). The calculated critical eccentricity is such that the finalrotation period satisfies the condition P_{rot} ∕P_{⋆} ≤ 0.95 to retain solely asynchronous states clearly separated from the synchronisation. In the eccentricity maps (top panels), the yellow colour designates the region of the nonresonant regime, where the critical eccentricity is high (log_{10} e_{AR} ≈−0.6). Conversely, the blue shades designate the regions where e_{AR} is low owing to the resonances of the oceanic tidal response. In the finalrotation maps, the yellow colour designates states that are close to synchronisation, while the dark blue colour corresponds to the 3/2 spinorbit resonance, induced by the degree3 eccentricity term. Values of parameters used for the solid part are given by Table 1. In both cases, the stellar mass is M_{⋆} = 0.09 M_{⊙}, the initial rotation rate is , and 10 Hough modes are taken into account in the calculation of the tidal response. 

In the text 
Fig. 9 Evolution of normalised rotation period P_{rot}∕P_{⋆} as function of time (yr) in logarithmic scale. This is given in the case of an Earthsized planet orbiting the TRAPPIST1 star with a sixday orbital period. The Rayleigh drag frequency of the ocean is set to σ_{R} = 10^{−5} s^{−1} for cases of strong drag (top panels), and to σ_{R} = 10^{−6} s^{−1} for cases of weak drag (bottom panels). The planet eccentricity takes the values e = 0.063, 0.2, and 0.3 from left to right. In each panel, a set of planets with various ocean depths is studied. These configurations are parametrised by , where the ocean depth is given in kilometers. The parameter γ is sampled from −2 (the driest case of the sample, solid orange line) to 0 (the most humid case, solid blue line). Values used for the planet parameters are given by Table 1. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.