© P. Auclair-Desrotour et al. 2019

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 TRAPPIST-1 system (Gillon et al. 2017), where an ultra-cool M-dwarf star harbours seven Earth-sized 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 TRAPPIST-1 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 free-surface 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 so-called 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 spin-orbit synchronous rotation. Unlike solid tides, which are not strongly dependent on the tidal frequency, oceanic tides exhibit a frequency-resonant behaviour that enables large variations of the tidal torque (e.g. Webb 1980; Tyler 2011; Chen et al. 2014; Matsuyama 2014; Auclair-Desrotour et al. 2018).

In thin fluid shells, such behaviour results from the propagation of gravito-inertial surface modes (or barotropic modes) forced by the tidal gravitational potential, which corresponds to the solutions of Laplace’s tidal equation (e.g. Longuet-Higgins 1968; Webb 1980). In deeper oceans, stable stratification leads to the propagation of internal gravito-inertial waves, which are restored by the Archimedean force and provide a supplementary resonant contribution to the oceanic tidal torque (e.g. Tyler 2011; Auclair-Desrotour 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 self-attraction 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; Auclair-Desrotour 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 planet-perturber 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. Auclair-Desrotour et al. 2018). In this configuration, the semidiurnal tidal component predominates and drives the planet towards the spin-orbit 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; Auclair-Desrotour 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 spin-orbit 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 two-layer 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 spin-orbit rotation equilibria were also retrieved in the case of giant planets by studies using a visco-elastic 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 spin-orbit 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 quasi-circular 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 TRAPPIST-1-like 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, self-attraction, and deformation of solid regions self-consistently 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 TRAPPIST-1 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; Castillo-Rogez et al. 2011; Efroimsky 2012) with values of parameters provided by a spectral code that computes the visco-elastic 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 self-consistently 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 high-frequency 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 Sundberg-Cooper rheology (e.g. Renaud & Henning 2018).

The ability of eccentricity tides to drive the planet away from the spin-orbit 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 spin-orbit resonances, e < e_AR means that the planet is driven towards the spin-orbit 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 quasi-adiabatic asymptotic regime.

In Sect. 6, the final rotation states of Earth and super-Earth-sized 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 spin-orbit 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 bi-layered 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 planeto-centric 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 co-rotating 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 unit-vector 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 core-ocean 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⁻⁵ s⁻¹ (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 spin-orbit synchronous rotation, eccentricity tides tend to desynchronise it, and thereby enable the existence of non-synchronised 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 Rc and a thin uniform ocean of external radius Rp 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 right-hand 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 semi-major 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 so-calledHansen 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 star-planet distance far exceeds the radius of the planet (Mathis & Le Poncin-Lafitte 2009). They can thus be ignored in the framework of this study. Quadrupolar terms (l = m = 2) are associated with the eccentricity frequencies (6)

and expressed as (7)

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 high-eccentricity 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 super-synchronous and corresponds to . In this case, the planet is torqued towards the spin-orbit 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 sub-synchronous 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 degree-s 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; Auclair-Desrotour 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 visco-elastic body of unrelaxed shear modulus μ. This body undergoes a time-periodic 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; Castillo-Rogez 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 olivine-rich 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. Castillo-Rogez 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 zero-frequency 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 high-frequency 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⁴–10⁵ days; for the Earth τ_M is about 500yr, see e.g. Efroimsky 2012), which is far greater than usual tidal periods (P_T ~ 10⁰−10² days). As a consequence, the imaginary part of drops rapidly as σ increases in the high-frequency regime, which leads to underestimating the energy that is tidally dissipated in the planet interior (this effect has repercussions on the frequency-dependence 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 Sundberg-Cooper rheology, which implies a larger tidal dissipation around a critical temperature and frequency (Renaud & Henning 2018).

As discussed by early studies (Castillo-Rogez 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 Castillo-Rogez 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 elasto-gravitational 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, non-rotating, 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 visco-elastic 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 frequency-spectra 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 TRAPPIST-1 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 degree-l 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 self-gravitational 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 degree-l Love numbers (e.g. Munk & MacDonald 1960) are given by (15)

where we have introduced the dimensionless effective rigidity (16)

with (17)

The contribution of the degree-l 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 low-frequency regime, where increases with σ linearly. A maximum is reached for P_T ~ τ_M. In the high-frequency 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 degree-2 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.

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 self-gravitational potential and a displacement of the oceanic floor, denoted by ξ_c. We take this solid-ocean 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 degree-l 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 eigenfunction-eigenvalue 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 degree-n 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 degree-n mode equalises the characteristic propagation velocity of large-wavelength surface gravity waves. The non-frictional case has been thoroughly discussed in early studies (e.g. Longuet-Higgins 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 g-modes (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 non-rotating case. For , another family of functions appear. These functions are generally called rotational modes or r-modes, and develop outside of the equatorial band where g-modes 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; Auclair-Desrotour 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, g-modes and r-modes tend to merge together and converge towards the functions of the non-rotating 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 zero-frequency 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¹, 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 degree-l 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)

the being given by (39)

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 degree-n 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 non-synchronised 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 degree-l 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 spin-axis 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)

with (48)

except in the zero-frequency 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 self-attraction 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 degree-2 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; Auclair-Desrotour et al. 2018), (50)

the notation and designating the complex eigenfrequencies of the degree-n 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 quasi-adiabatic asymptotic regime (), where the imaginary part of is negligible, and for g-modes, given that the latter are associated with positive . By introducing the characteristic frequency of the degree-n 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 frequency-spectra of the tidal torque. The frequency atwhich the resonant peak of a mode reaches a maximum is expressed as (56)

In the quasi-adiabatic regime (), it reduces to (57)

The corresponding maximum of for the resonance associated with the degree-n Hough mode is denoted by and obtained by substituting σ with Eq. (57) in Eq. (55). If the contributions of non-resonant components are neglected, we thus obtain (58)

which shows that, as a first approximation, the maximum of the peak associated with the degree-n 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 non-resonant background.It can thus generate a torque strong enough to compensate the predominating contribution of the semidiurnal component in the low-eccentricity 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 super-synchronous 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 degree-0 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 degree-0 Hough modes. The corresponding frequencies of these resonant modes are ±σ₀ 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 super-synchronous frequency range (σ_2,2 > 0) and the other one in the sub-synchronous frequency range (σ_2,2 < 0). Synchronisation is symmetrical and is due to the non-rotating approximation assumed in Fig. 4. We emphasise the fact that the peak of positive torque appearing in the super-synchronous regime results from the forcing of the resonance located in the sub-synchronous frequency range by a prograde eccentricity tidal potential (σ_2,s = −σ₀). This configuration corresponds for instance to case 1 of Fig. 2, where the degree-3 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 degree-3 eccentric term can coincide with the prograde resonant frequency of the ocean (sub-synchronous blue cone), that is . When the depth of the ocean increases (from d to f), keeping in the range 0-n_⋆, this leads to an excitation of the ocean resonant mode only with the degree-4 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 σ₀ = n_⋆∕2. In this case, the peaks created by the resonance of the degree-3 and degree-2 (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 degree-3 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 degree-3 term to the degree-4 term that occurs while the ocean depth keeps increasing. The peak of positive tidal torque now results from the amplification of the degree-4 eccentricity term by the resonance associated with the tidal oceanic mode. It is smaller than that generated by the degree-3 term since it scales quadratically with the forcing tidal potential (see Eqs. (7) and (49)), which decays while the degree-s component increases in the low-eccentricity regime (for e ≳ e_trans, non-linearities lead eccentricity components of higher degrees to predominate; see e.g. Ogilvie 2014, Fig. 3).

Similar to the peak created by the degree-3 component, the peak associated with the degree-4 eccentricity term moves towards synchronisation while the ocean depth increases, as shown by configurations (e) and (f). When σ₀ becomes greater than2n_⋆, the term generating the peak of positive tidal torque switches from the degree-4 to the degree-5 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 degree-0 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 non-resonant background level. This picture is also completed by the action of solid tides, which enable the existence of tidally-locked 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 degree-0 Hough mode of frequency σ₀ since in the frequency range of interest. The real tidal response is composed of g- and r-modes (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 quasi-adiabatic 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 quasi-adiabatic regime violate the small perturbation approximation. As a consequence, the tidal response becomes non-linear 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⁻⁵ s⁻¹ (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 Auclair-Desrotour 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⁻⁵ s⁻¹ 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⁻⁷ and 10⁻⁶ s⁻¹ in addition to 10⁻⁵ s⁻¹ in numerical calculations.

In spite of the limitations mentioned above, studying the quasi-adiabatic 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 super-synchronous 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 degree-0 mode, of frequency σ0, in the non-rotating 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 super-synchronous regime (σ2,2 > 0), downwardblue peaks tend to drive the planet towards tidal locking in spin-orbit 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 tightly-packed systems such as that hosted by the TRAPPIST-1 ultra-cool dwarf star (Gillon et al. 2017; Grimm et al. 2018), where planet e exhibits a 6.10-day orbital period (Gillon et al. 2017).

We thus consider throughout the whole study the case of an Earth-sized planet orbiting TRAPPIST-1 (M_⋆ = 0.09 M_⊙, see Van Grootel et al. 2018) with a six-day 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 km-deep uniform ocean, which corresponds to the resonance of one of the dominating oceanic modes for the degree-3 eccentricity term in the range 0 < σ_2,2 < n_⋆. We set the typical frequency of the Rayleigh friction to σ_R = 1.0 × 10⁻⁶ s⁻¹. 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 spin-orbit 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 frequency-dependence 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 spin-orbit 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 degree-3 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 sub-inertial 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⁻⁶ s⁻¹.

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 frequency-resonant 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 Earth-sized planet orbiting TRAPPIST-1 (M⋆ = 0.09 M⊙) with a six-day 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) spin-orbit 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 tidally-locked 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 non-resonant 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 spin-orbit 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 non-synchronised 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 non-resonant background, which is the sum of non-resonant 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 degree-s forcing frequencies are negative while the semidiurnal forcing frequency is positive. Assuming condition (i), we suppose that the degree-n mode is resonant. This leads to (59)

where σ_2,s = 2Ω − sn_⋆ is the degree-s eccentricity frequency defined by Eq. (6) and σ_n the eigenfrequency of the degree-n 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 close-in terrestrial planets. Conversely, this limit is inversely proportional to the eigenvalue of the mode. In the quasi-adiabatic 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 degree-0 mode will be expressed preferentially with respect to those associated with other modes. A high-degree 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 self-attraction (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 non-resonant 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 non-resonant 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 degree-3 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 core-ocean 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 degree-s quadrupole component (the resonance being associated with the degree-n 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 degree-0 Hough mode, which means that except for n = 0. As a consequence, the total tidal torque is dominated by the contribution of the degree-0 Hough mode outside of the resonances associated with other modes (Auclair-Desrotour 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 degree-s eccentricity component (i) is resonant, that is σ_2,s = −σ₀ and (Eq. (59)); (ii) made use of Eqs. (55) and (58) for the degree-0 Hough mode; and (iii) introduced the notation .

The condition on the eigenfrequency σ₀ resulting from the preceding expression is implicit in the general case since the overlap coefficients ( and ) and eigenvalue () associated with the degree-0 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)

or (64)

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 degree-s 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 degree-s eccentricity component of the tidal gravitational potential. By recalling Eq. (60), we hence observe that the degree-s eccentricity term cannot drive the planet away from spin-orbit 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 degree-0 Hough mode (σ₀). 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 degree-3 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 tidally-locked 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⁻⁵ s⁻¹. This leads to e_AR ≈ 0.1 for σ_R = 10⁻⁵ s⁻¹ and e_AR ≈ 0.01 for σ_R = 10⁻⁶ s⁻¹. 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 non-resonant 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 degree-s is amplified by the resonance associated with the degree-0 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 = σ₀∕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 non-rotating 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 low-eccentricity regime, we focus on the effect of the degree-3 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 left-hand 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 non-resonant background level, this scaling law switches to . As the non-resonant background level of the oceanic component of the tidal torque scales as ∝ σ_R (see e.g. Auclair-Desrotour 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 final-rotation 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 final-rotation 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 pseudo-synchronous 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 spin-orbit resonances (Goldreich 1966; Goldreich & Peale 1966).

6.1 Tidal locking in asynchronous states of equilibrium

We treat the cases of an Earth-sized planet and of a 10-M_⊕ super-Earth orbiting the TRAPPIST-1 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 super-Earth orbits the star with a six-day orbital period. As the Rayleigh drag frequency cannot be accurately specified for an exoplanet of unknown topography, we perform the calculations for σ_R = 10⁻⁷, 10⁻⁶, 10⁻⁵ s⁻¹, 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 final-rotation 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 super-synchronous 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 spin-orbit synchronous rotation or the asynchronous rotation state induced by the degree-3 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 final-rotation 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 super-Earth), 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 spin-orbit 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 10-M_⊕ super-Earth.

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₁₀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 degree-3 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 degree-3 eccentricity term, thus leading to the most important decay of the critical eccentricity. It then meets the frequency of the degree-4 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⁻⁷ s⁻¹, 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 six-day orbital period and the Rayleigh frequencies 10⁻⁵ and 10⁻⁶ s⁻¹ (top left and middle panels), we observe the transition from the frictional regime () to the quasi-adiabatic 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⁻⁶ s⁻¹.

Finally, we note that the results obtained for the super-Earth 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 super-Earth. 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 degree-2 Love number is well approximated by its high-frequency 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 super-Earth 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 spin-orbit 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 high-eccentricity regime. The important parameters for these states are the ocean parameters (depth, density, and Rayleigh drag frequency) and the orbital period.

Fig. 6 Normalised final-rotation period Prot∕P⋆ as a function of eccentricity (horizontal axis) and ocean depth (vertical axis) in logarithmic scales. Two planets are considered: an Earth-sized planet with six-day (top panels) and one-day (middle panels) orbital periods, and a 10 M⊕ -super-Earth with a six-day 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 six-day-period Earth-sized planet with a weak drag(top right panel), consistent with the discussion on Fig. 4.

6.2 Capture in 1:1 spin-orbit 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, non-axisymmetric 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 spin-orbit 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 right-hand 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 non-rotationally symmetric deformation of the planet. It is annihilated for a spherically symmetric body, where A = B = C.

In the vicinity of the spin-orbit resonance s∕2, the torque due to triaxiality is dominated by the degree-s 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² (77)

which is oscillatory and leads to libration when the planet is in spin-orbit 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 spin-orbit ratio of the resonance by the triaxial torque. In this range, the existence of rotation equilibria distinct from the spin-orbit 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 spin-orbit 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 spin-orbit 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 non-synchronous 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⁻⁴ for Mercury, as derived by Ribas et al. (2016) from the gravity moments measured by Smith et al. (2012), ~ 2 × 10⁻⁵ for the Earth, and ~6 × 10⁻⁶ 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 spin-orbit 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 Earth-sized planet of six-day orbital period with σ_R = 10⁻⁶ s⁻¹ (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 spin-orbit 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 spin-orbit 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 spin-orbit 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 spin-orbit 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 spin-orbit 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 (Mp = 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 spin-orbit 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 super-Earth 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 spin-orbit 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 final-rotation period satisfies the condition P_rot∕P_⋆ ≤ 0.95. This means that the departure between the final-rotation 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⁻³ (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 final-rotation 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⁻⁵ s⁻¹, left), and a weak drag (σ_R = 10⁻⁶ s⁻¹, 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, yellow-green areas designate the non-resonant regime, where the critical eccentricity is high (log₁₀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 degree-3 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 non-resonant 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 quasi-adiabatic () and non-rotating (no Coriolis effects) approximations, which leads to differences with the results of numerical calculations, and particularly in the strong drag configuration (σ_R = 10⁻⁵ s⁻¹).

We now move to the maps representing the final-rotation rates associated with the computed critical eccentricity (bottom panels). Yellow areas designate states of equilibrium close to synchronisation, and blue areas the 3: 2 spin-orbit resonance. In order to avoid confusion when looking at these maps, one should bear in mind that they show both the low-eccentricity regime, which is the subject of the analytical theory detailed in Sects. 4 and 5, and the high-eccentricity regime, where the derived results do not apply since it is beyond the scope of this work.

We shall also stress here a somehow counter-intuitive 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 spin-orbit 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 low-eccentricity 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. Auclair-Desrotour 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 eAR (top panels) and associated normalised final-rotation period Prot∕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 final-rotation period satisfies the condition Prot ∕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 non-resonant regime, where the critical eccentricity is high (log10 eAR ≈−0.6). Conversely, the blue shades designate the regions where eAR is low owing to the resonances of the oceanic tidal response. In the final-rotation maps, the yellow colour designates states that are close to synchronisation, while the dark blue colour corresponds to the 3/2 spin-orbit resonance, induced by the degree-3 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 Earth-sized planet orbiting the TRAPPIST-1 star with a six-day 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⁻⁵ s⁻¹ and σ_R = 10⁻⁶ s⁻¹, 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 spin-orbit 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 Bulirsch-Stoer 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³ 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 spin-orbit 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⁻⁶ s⁻¹, 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 spin-orbit 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 non-synchronised 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 non-resonant level.

In the third regime (e = 0.3, right panels) the contribution of non-resonant 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 spin-orbit 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 self-consistent 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, self-attraction, 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 high-frequency 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 quasi-adiabatic 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 super-Earth. 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 spin-orbit 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 degree-3 (or eastward propagating) eccentricity term, which is the largest one in the low-eccentricity 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 low-eccentricity 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 spin-orbit synchronous rotation; and (iii) the high-eccentricity regime (e ≳ 0.3), where the eccentricity is strong enough to make planets converge towards the 3: 2 spin-orbit 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 TRAPPIST-1 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 close-in, 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 (Anglada-Escudé 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 spin-orbit 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 TRAPPIST-1 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 non-linear 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 quasi-adiabatic regime, below σ_R ~ 10⁻⁶ s⁻¹. A second limitation is related to the shallow water approximation, which eludes the role played by the ocean-internal structure on the tidally dissipated energy. Although the formalism can easily be extended to stably stratified deep oceans if the effects of ocean loading, self-attraction, and deformation of the solid regions are ignored (Auclair-Desrotour 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 free-surface 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 Prot∕P⋆ as function of time (yr) in logarithmic scale. This is given in the case of an Earth-sized planet orbiting the TRAPPIST-1 star with a six-day 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 visco-elastic 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)

and (C.4)

Appendix D Nomenclature

References

All Tables

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 Rc and a thin uniform ocean of external radius Rp 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 super-synchronous 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 degree-0 mode, of frequency σ0, in the non-rotating 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 super-synchronous regime (σ2,2 > 0), downwardblue peaks tend to drive the planet towards tidal locking in spin-orbit 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 Earth-sized planet orbiting TRAPPIST-1 (M⋆ = 0.09 M⊙) with a six-day 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) spin-orbit synchronousrotation.

In the text

Fig. 6 Normalised final-rotation period Prot∕P⋆ as a function of eccentricity (horizontal axis) and ocean depth (vertical axis) in logarithmic scales. Two planets are considered: an Earth-sized planet with six-day (top panels) and one-day (middle panels) orbital periods, and a 10 M⊕ -super-Earth with a six-day 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 six-day-period Earth-sized 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 spin-orbit 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 spin-orbit 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 (Mp = M⊕, P⋆ = 6 days, and σR = 10−6 s−1).

In the text

Fig. 8 Logarithm of critical eccentricity eAR (top panels) and associated normalised final-rotation period Prot∕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 final-rotation period satisfies the condition Prot ∕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 non-resonant regime, where the critical eccentricity is high (log10 eAR ≈−0.6). Conversely, the blue shades designate the regions where eAR is low owing to the resonances of the oceanic tidal response. In the final-rotation maps, the yellow colour designates states that are close to synchronisation, while the dark blue colour corresponds to the 3/2 spin-orbit resonance, induced by the degree-3 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 Prot∕P⋆ as function of time (yr) in logarithmic scale. This is given in the case of an Earth-sized planet orbiting the TRAPPIST-1 star with a six-day 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