Final spin states of eccentric ocean planets

TL;DR

This paper models how eccentricity tides influence the spin states of ocean-bearing rocky planets orbiting cool dwarf stars, revealing how oceanic resonances can lead to asynchronous rotation at lower eccentricities.

Contribution

It develops a comprehensive global model combining oceanic and solid-body tides to analyze resonance effects on planetary rotation states.

Findings

Resonances significantly affect the critical eccentricity for asynchronous rotation.

Critical eccentricity can be reduced from ~0.3 to ~0.06 or lower due to oceanic tidal resonances.

The model provides analytical and numerical tools to explore tidal responses across parameter space.

Abstract

Eccentricity tides generate a torque that can drive an ocean planet towards asynchronous rotation states of equilibrium when enhanced by resonances associated with the oceanic tidal modes. We investigate the impact of eccentricity tides on the rotation of rocky planets hosting a thin uniform ocean and orbiting cool dwarf stars such as TRAPPIST-1, with orbital periods ~1-10 days. Combining the linear theory of oceanic tides in the shallow water approximation with the Andrade model for the solid part of the planet, we develop a global model including the coupling effects of ocean loading, self-attraction, and deformation of the solid regions. We derive from this model analytic solutions for the tidal Love numbers and torque exerted on the planet. These solutions are used with realistic values of parameters provided by advanced models of the internal structure and tidal oscillations of solid bodies to explore the parameter space both analytically and numerically. Our model allows us to fully characterise the frequency-resonant tidal response of the planet, and particularly the features of resonances associated with the oceanic tidal modes (eigenfrequencies, resulting maxima of the tidal torque and Love numbers) as functions of the planet parameters (mass, radius, Andrade parameters, ocean depth and Rayleigh drag frequency). Resonances associated with the oceanic tide decrease the critical eccentricity beyond which asynchronous rotation states distinct from the usual spin-orbit resonances can exist. We provide an estimation and scaling laws for this critical eccentricity, which is found to be lowered by roughly one order of magnitude, switching from ~0.3 to ~0.06 in typical cases and to ~0.01 in extremal ones.

Tables

Table 1. Table 1 : Parameters of the Andrade model derived from the elasto-gravitational theory (Tobie et al., 2005 ; Breton et al., 2018 ) .

Parameter	Earth	Super-Earth
Planet mass $M_p$ ($M_{♁}$)	1.0	10
Planet radius $R_p$ ($R_{♁}$)	1.0	1.86
Shear modulus $\mu$ (GPa)	14.6	75.2
Rheological exponent $\alpha$	0.25	0.25
Maxwell time ${\tau }_M$ (yr)	$741$	$1.87\times {10}^5$
Andrade time ${\tau }_A$ (yr)	$2.19\times {10}^4$	$1.66\times {10}^6$

Table 2. Table 2: List of the notations used along this work and their designations in order of appearance in the text.

Symbol	Definition	Reference
$e$	Planet eccentricity	Sect. 1
$e_{\mathrm{AR}}$	Critical eccentricity for asynchronous rotation	Sect. 1
$M_p$	Planet mass	Sect. 2
$R_p$	Planet radius	Sect. 2
$H_{\mathrm{oc}}$	Ocean depth	Sect. 2
${\rho }_{\mathrm{oc}}$	Ocean density	Sect. 2
$n_{\star }$	Planet mean motion	Sect. 2
${𝒏}_{\star }$	Orbital angular momentum of the planet	Sect. 2
${ℛ}_G$	Planeto-centric reference frame	Sect. 2
${ℛ}_p$	Reference frame co-rotating with the planet	Sect. 2
$𝛀$	Spin vector	Sect. 2
$\Omega$	Rotation rate of the planet	Sect. 2
$t$	Time	Sect. 2
$r$	Radial coordinate (spherical coordinates)	Sect. 2
$\theta$	Colatitude (spherical coordinates)	Sect. 2
$\phi$	Longitude (spherical coordinates)	Sect. 2
$({𝐞}_r,{𝐞}_{\theta },{𝐞}_{\phi })$	Spherical unit-vector basis	Sect. 2
$g$	Planet surface gravity	Sect. 2
${\rho }_p$	Planet mean density	Sect. 2
$M_{\mathrm{oc}}$	Mass of the ocean	Sect. 2
$M_c$	Mass of the solid part	Sect. 2
$R_c$	Radius of the solid part	Sect. 2
${\rho }_c$	Mean density of the solid part	Sect. 2
${\sigma }_R$	Rayleigh drag frequency	Sect. 2
${\tau }_R$	Typical timescale of core-ocean coupling by viscous friction	Sect. 2
$r_{\star }$	Star-planet distance	Eq. 1
$M_{\star }$	Mass of the host star	Eq. 1
${\hat{U}}_{\star }$	Gravitational potential of the host star	Eq. 1
$\hat{}$	Symbol used to highlight real quantities	Eq. 1
$\stackrel{~}{}$	Symbol used to highlight complex quantities	Eq. 1
${\hat{U}}_T$	Real tidal gravitational potential	Eq. 2
$\Re$	Real part of a complex number	Eq. 3
$\Im$	Imaginary part of a complex number	Eq. 3
$i$	Imaginary number	Eq. 3
$U_T$	Complex tidal gravitational potential	Eq. 3
$l$	Latitudinal degree (spherical harmonics)	Eq. 3
$m$	Longitudinal degree (spherical harmonics)	Eq. 3
$s$	Eccentricity degree	Eq. 3
$\sigma$	Forcing frequency	Eq. 3
${\sigma }_{m,s}$	Forcing tidal frequency of the mode associated with the doublet $(m,s)$	Eq. 3
$P_l^m$	Normalised associated Legendre functions	Eq. 3
$U_{T;l}^{m,s}$	$(l,m,s)$-component of the forcing tidal gravitational potential at the planet surface	Eq. 3
$a$	Semi-major axis of the planet	Eq. 4
$A_{l,m,s}$	Dimensionless eccentricity functions	Eq. 4
${\delta }_{l,k}$	Kronecker symbol	Eq. 5
	Coefficient equal to 1 for , 0 otherwise	Eq. 5
${\overline{P}}_l^m$	Unnormalised associated Legendre functions	Eq. 5
$X_s^{l,m}$	Hansen coefficients	Eq. 5
${\sigma }_{2,s}$	Forcing frequency of the degree-$s$ eccentricity term	Eq. 6
$U_{T;2}^{2,s}$	Degree-$s$ component of the perturbing tidal potential	Eq. 7
$U_{T;2}^{2,2}$	Semidiurnal component of the perturbing tidal potential	Eq. 8
$U_{T;2}^{2,1}$	Degree-1 eccentricity component of the perturbing tidal potential	Eq. 9
${\sigma }_{2,1}$	Forcing frequency of the degree-1 eccentricity component	Eq. 9
$U_{T;2}^{2,3}$	Degree-3 eccentricity component of the perturbing tidal potential	Eq. 10
${\sigma }_{2,3}$	Forcing frequency of the degree-3 eccentricity component	Eq. 10
$e_{\mathrm{trans}}$	Transition eccentricity	Sect. 3.1
$\mu$	Effective unrelaxed shear modulus of the solid core	Sect. 3.2
$\underset{¯}{\underset{¯}{𝜻}}$	Stress tensor	Eq. 12
$\underset{¯}{\underset{¯}{𝜺}}$	Strain tensor	Eq. 12
$\tilde{\mu }$	Complex shear modulus of the solid part	Eq. 13
$\Gamma$	Gamma function	Eq. 14
$\alpha$	Rheological exponent of the material in the Andrade model	Eq. 14
${\tau }_M$	Maxwell relaxation time of the material	Eq. 14
$\eta$	Viscosity of the material	Eq. 14
${\tau }_A$	Andrade time of the material in the Andrade model	Eq. 14
$P_T$	Tidal period	Sect. 3.2
$\beta$	Andrade parameter	Sect. 3.2
$k_l^{\sigma }$	Gravitational Love number of the degree-$l$ mode	Eq. 15
$h_l^{\sigma }$	Displacement Love number of the degree-$l$ mode	Eq. 15
${𝔨}_l^{\sigma }$	Gravitational load Love number of the degree-$l$ mode	Eq. 15
${𝔥}_l^{\sigma }$	Displacement load Love number of the degree-$l$ mode	Eq. 15
${\tilde{\mu }}_l^{\sigma }$	Dimensionless effective rigidity	Eq. 16
$A_l$	Dimensionless rigidity coefficient associated with the degree-$l$ component	Eq. 17
${\xi }_c$	Vertical displacement of the oceanic floor	Sect. 3.3
$𝑽$	Velocity vector of the horizontal component of tidal flows	Sect. 3.3
$V_{\theta }$	Latitudinal component of the velocity field	Sect. 3.3
$V_{\phi }$	Longitudinal component of the velocity field	Sect. 3.3
${\xi }_{\mathrm{oc}}$	Vertical displacement of the ocean surface	Sect. 3.3
$\eta$	Variation of the ocean thickness	Sect. 3.3
${\partial }_X$	Partial derivative with respect to $X$	Eq. 18
$\Psi$	Perturbed potential encompassing the gravitational forcing and coupling with the solid part	Eq. 18
${\nabla }_h$	Horizontal gradient operator	Eq. 20
${\nabla }_h\cdot$	Horizontal divergence operator	Eq. 21
${\Psi }^{m,\sigma }$	Fourier component of the perturbed potential	Eq. 22
${𝑽}^{m,\sigma }$	Fourier component of the velocity vector	Eq. 22
${\eta }^{m,\sigma }$	Fourier component of the vertical variation of the ocean depth	Eq. 22
${\Psi }_l^{m,\sigma }$	Degree-$l$ component of ${\Psi }^{m,\sigma }$ expanded in series of associated Legendre functions	Eq. 23
${𝑽}_l^{m,\sigma }$	Degree-$l$ component of ${𝑽}^{m,\sigma }$ expanded in series of associated Legendre functions	Eq. 23
${\eta }_l^{m,\sigma }$	Degree-$l$ component of ${\eta }^{m,\sigma }$ expanded in series of associated Legendre functions	Eq. 23
$\tilde{\sigma }$	Complex tidal frequency	Eq. 24
$\tilde{\nu }$	Complex spin parameter	Eq. 24
${ℒ}^{m,\tilde{\nu }}$	Laplace’s tidal operator	Eq. 27
$n$	Degrees of Hough functions	Eq. 29
${\Psi }_n^{m,\sigma }$	Degree-$n$ component of ${\Psi }^{m,\sigma }$ expanded in series of Hough functions	Eq. 29
${\Theta }_n^{m,\tilde{\nu }}$	Degree-$n$ Hough function	Eq. 29
${\Lambda }_n^{m,\tilde{\nu }}$	Eigenvalue associated with the degree-$n$ Hough function	Eq. 30
$h_n^{m,\sigma }$	Equivalent depth of the degree-$n$ Hough mode	Eq. 31
$A_{n,l}^{m,\tilde{\nu }}$	Coefficients of Hough functions expanded in series of the $P_l^m$	Eq. 32
$B_{l,n}^{m,\tilde{\nu }}$	Coefficients of the $P_l^m$ expanded in series of Hough functions	Eq. 33
$C_{l,n,k}^{m,\tilde{\nu }}$	Overlap coefficients weighting the degree-$n$ component of the oceanic tidal torque	Eq. 34
${\gamma }_{U;l}^{\sigma }$	Tilt factor associated with the tidal gravitational forcing of the perturber	Eq. 35
${\gamma }_{\xi ;l}^{\sigma }$	Tilt factor associated with the distortion of the oceanic layer	Eq. 36
${\sigma }_{n,k}$	Complex characteristic frequencies of the planet	Eq. 41
$k_{h;n}$	Horizontal wavenumber of the degree-$n$ Hough mode	Eq. 41
$F_n^{m,\sigma }$	Component of the force vector inducing the tidal perturbation	Eq. 42
${𝒯}_p$	Tidal torque exerted on the planet	Eq. 44
$k_{p;l}^{m,\sigma }$	Quadrupolar component of the effective gravitational Love number of the planet	Eq. 45
${𝒯}_{\mathrm{sol}}$	Tidal torque exerted on the solid part in absence of oceanic layer	Eq. 46
$k_2^{\sigma }$	Quadrupolar tidal Love number of the solid part in absence of oceanic layer	Eq. 47
$A_2$	Dimensionless rigidity coefficient associated with the $l=2$ component	Eq. 48
${𝒯}_{\mathrm{oc}}$	Tidal torque exerted on the ocean in case of solid part of infinite rigidity	Eq. 49
${\tilde{\sigma }}_n^{+}$	Complex eigenfrequency of the degree-$n$ Hough mode in the positive-frequency range	Eq. 50
${\tilde{\sigma }}_n^{-}$	Complex eigenfrequency of the degree-$n$ Hough mode in the negative-frequency range	Eq. 50
${\sigma }_n$	Characteristic frequency of the degree-$n$ surface gravity mode	Eq. 54
${\sigma }_{\max ;n}$	Frequency at which the resonant peak of a mode reaches a maximum	Eq. 56
${\Im \left\{k_{\mathrm{oc};2}^{2,\sigma }\right\}|}_{\max ;n}$	Maximum of $\Im \left\{k_{\mathrm{oc};2}^{2,\sigma }\right\}$ reached by the resonance peak	Eq. 58
${\sigma }_{2,s}^{▴}$	Degree-$s$ eccentricity frequency associated with a peak	Sect. 4.3
${\Omega }_{\mathrm{eq};s}$	Rotation rate corresponding to the degree-$s$ spin-orbit resonance	Sect. 4.3
$P_{\star }$	Orbital period of the planet	Sect. 4.4
$M_{☉}$	Mass of the Sun	Sect. 4.4
${\tilde{\nu }}_{2,2}$	Spin parameter associated with the semidiurnal tidal component	Eq. 62
${\gamma }_s$	Parameter controlling the existence of asynchronous states in the case of pure oceanic tide	Eq. 65
$X$	Normalised frequency	Eq. 69
$𝒜$	Dimensionless constant	Eq. 69
$u_{\alpha }$	Maximum of the polynomial function defined by Eq. 69	Eq. 71
$P_{\mathrm{rot}}$	Rotation period of the planet	Sect. 6.1
$H_{\mathrm{oc};s}$	Ocean depth associated with the resonance of the degree-$s$ eccentricity component	Sect. 6.1
${ℬ}_c$	Constant factor in the expression of the solid tidal torque in the non-resonant regime	Eq. 73
${𝒯}_p^{\mathrm{tide}}$	Tidal torque exerted on the planet	Eq. 75
${𝒯}_p^{\mathrm{tri}}$	Torque due to triaxiality	Eq. 75
$\dot{}$	Time derivative	Eq. 75
$A,B,C$	Principal moments of inertia of the planet (in increasing magnitude)	Eq. 76
${\gamma }_{2,s}$	angle associated with the degree-$s$ perturbation component	Eq. 76
$\vartheta$	Planet rotation angle	Eq. 76
$M$	Planet mean anomaly	Eq. 76
$\ddot{}$	Second order time derivative	Eq. 78
$\Delta$	Width of the resonance	Eq. 79
${𝒯}_{\mathrm{even}}^{\mathrm{tide}}$	Even component of the tidal torque	Eq. 81
${𝒯}_{\mathrm{odd}}^{\mathrm{tide}}$	Odd component of the tidal torque	Eq. 82
${𝒫}_{\mathrm{cap};s}$	Capture probability of the $s/2$ spin-orbit resonance	Eq. 83
${\Delta }_P$	Normalised width of the $1:1$ spin-orbit resonance	Eq. 84
$\gamma$	Logarithm of the ocean depth	Sect. 8
$\check{C}$	Normalised moment of inertia	Sect. 8
${\tau }_{\mathrm{ev}}$	Evolution timescale of the planet spin rotation	Eq. 89
$g_l^m$	Eccentricity and mean anomaly function	Eq. 94
$v$	True anomaly	Eq. 94
$N$	Order of the fast Fourier transform (FFT) in the calculation of Hansen coefficients	Appendix B
$K$	Truncation index of the series of Hansen coefficients calculated using the FFT method	Appendix B
${\widetilde{𝒯}}_p^{\mathrm{tri}}$	Complex triaxial torque	Eq. 100

Equations

$${\hat{U}}_{\star }\left(r,r_{\star }\right)=\frac{G{M}_{\star }}{\left|r-r_{\star }\right|}-\frac{G{M}_{\star }}{r_{\star }^2}r\cos \theta ,$$

$${\hat{U}}_{\mathrm{T}}\left(\theta ,\phi ,r_{\star }\right)\equiv {\hat{U}}_{\star }\left(R_{\mathrm{p}},\theta ,\phi ,r_{\star }\right)-\frac{G{M}_{\star }}{r_{\star }},$$

$$U_{\mathrm{T}}=\sum _{l=2}^{+\infty }\sum _{m=0}^l\sum _{s=-\infty }^{+\infty }{U}_{\mathrm{T};l}^{m,s}\left(R_{\mathrm{p}}\right)P_l^m\left(\cos \theta \right){\mathrm{e}}^{i\left[{\sigma }_{m,s}t+m\phi \right]},$$

$$U_{\mathrm{T};l}^{m,s}\left(R_{\mathrm{p}}\right)=\frac{G{M}_{\star }}{a}{\left(\frac{R_{\mathrm{p}}}{a}\right)}^l{A}_{l,m,s}\left(e\right),$$

$$A_{l,m,s}\left(e\right)=$$

$$\times \sqrt{\frac{2\left(l-m\right)!}{\left(2l+1\right)\left(l+m\right)!}}{\bar{P}}_l^m\left(0\right)X_s^{-\left(l+1\right),m}\left(e\right).$$

$${\sigma }_{2,s}=2\Omega -s{n}_{\star },$$

$$U_{\mathrm{T};2}^{2,s}=\sqrt{\frac{3}{5}}\frac{G{M}_{\star }}{a}{\left(\frac{R_{\mathrm{p}}}{a}\right)}^2{X}_s^{-3,2}\left(e\right),$$

$$U_{\mathrm{T};2}^{2,2}=\sqrt{\frac{3}{5}}\frac{G{M}_{\star }}{a}{\left(\frac{R_{\mathrm{p}}}{a}\right)}^2.$$

$$U_{\mathrm{T};2}^{2,1}=-\frac{e}{2}\sqrt{\frac{3}{5}}\frac{G{M}_{\star }}{a}{\left(\frac{R_{\mathrm{p}}}{a}\right)}^2,$$

$$U_{\mathrm{T};2}^{2,3}=\frac{7e}{2}\sqrt{\frac{3}{5}}\frac{G{M}_{\star }}{a}{\left(\frac{R_{\mathrm{p}}}{a}\right)}^2,$$

$$\begin{array}[]{ll}\displaystyle\frac{s}{2}n_{\star}<\Omega\leq n_{\star}&\mbox{if}\ s\leq 1,\\ \displaystyle n_{\star}\leq\Omega<\frac{s}{2}n_{\star}&\mbox{if}\ s\geq 3.\end{array}$$

$$\underline{\underline{\mathbf{\epsilon }}}=\frac{1}{2}\left[\nabla \mathbf{\xi }+{\left(\nabla \mathbf{\xi }\right)}^{\mathrm{T}}\right],$$

$$\underline{\underline{\mathbf{\zeta }}}=2\tilde{\mu }\underline{\underline{\mathbf{\epsilon }}},$$

$$\tilde{\mu }=\frac{\mu }{1+{\left(i\sigma {\tau }_{\mathrm{A}}\right)}^{-\alpha }\Gamma \left(1+\alpha \right)+{\left(i\sigma {\tau }_{\mathrm{M}}\right)}^{-1}},$$

$$\left\{k_l^{\sigma },h_l^{\sigma },{\mathfrak{k}}_l^{\sigma },{\mathfrak{h}}_l^{\sigma }\right\}=\frac{1}{1+{\tilde{\mu }}_l^{\sigma }}\left\{\frac{3}{2\left(l-1\right)},\frac{2l+1}{2\left(l-1\right)},-1,-\frac{2l+1}{3}\right\},$$

$${\tilde{\mu }}_l^{\sigma }\equiv A_l\frac{\tilde{\mu }}{\mu },$$

$$A_l\equiv \frac{4\left(2l^2+4l+3\right)\pi R_{\mathrm{c}}^4\mu }{3lG{M}_{\mathrm{c}}^2}.$$

$${\partial }_t\mathbf{V}+2\mathbf{\Omega }\wedge \mathbf{V}+{\sigma }_{\mathrm{R}}\mathbf{V}={\nabla }_{\mathrm{h}}\Psi ,$$

$${\partial }_t\eta +H_{\mathrm{oc}}{\nabla }_{\mathrm{h}}\cdot \mathbf{V}=0,$$

$${\nabla }_{\mathrm{h}}\equiv R_{\mathrm{p}}^{-1}\left[{\mathrm{e}}_{\theta }{\partial }_{\theta }+{\mathrm{e}}_{\phi }{\left(\sin \theta \right)}^{-1}{\partial }_{\phi }\right],$$

$${\nabla }_{\mathrm{h}}\cdot \mathbf{V}\equiv {\left(R_{\mathrm{p}}\sin \theta \right)}^{-1}\left[{\partial }_{\theta }\left(\sin \theta V_{\theta }\right)+{\partial }_{\phi }{V}_{\phi }\right],$$

$$\begin{array}[]{ll}\displaystyle\Psi=\sum_{m,\sigma}\Psi^{m,\sigma}\left(\theta\right){\rm e}^{i\left(\sigma t+m\varphi\right)},&\displaystyle\boldsymbol{V}=\sum_{m,\sigma}\boldsymbol{V}^{m,\sigma}\left(\theta\right){\rm e}^{i\left(\sigma t+m\varphi\right)},\\ \displaystyle\eta=\sum_{m,\sigma}\eta^{m,\sigma}\left(\theta\right){\rm e}^{i\left(\sigma t+m\varphi\right)}.&\end{array}$$

$$\begin{array}[]{ll}\displaystyle\Psi^{m,\sigma}=\sum_{l\geq m}\Psi_{l}^{m,\sigma}P_{l}^{m}\left(\cos\theta\right),&\displaystyle\boldsymbol{V}^{m,\sigma}=\sum_{l\geq m}\boldsymbol{V}_{l}^{m,\sigma}P_{l}^{m}\left(\cos\theta\right),\\ \displaystyle\eta^{m,\sigma}=\sum_{l\geq m}\eta_{l}^{m,\sigma}P_{l}^{m}\left(\cos\theta\right),&\end{array}$$

$$\begin{array}[]{lrl}\tilde{\sigma}\equiv\sigma-i\sigma_{\rm R},&\mbox{and}&\displaystyle\tilde{\nu}\equiv\frac{2\Omega}{\tilde{\sigma}},\end{array}$$

$$V_{\theta }^{m,\sigma }=-\frac{i}{R_{\mathrm{p}}\tilde{\sigma }\left(1-{\tilde{\nu }}^2{\cos }^2\theta \right)}\left(\frac{\mathrm{d}}{\mathrm{d}\theta }+m\tilde{\nu }\cot \theta \right){\Psi }^{m,\sigma },$$

$$V_{\phi }^{m,\sigma }=\frac{1}{R_{\mathrm{p}}\tilde{\sigma }\left(1-{\tilde{\nu }}^2{\cos }^2\theta \right)}\left(\tilde{\nu }\cos \theta \frac{\mathrm{d}}{\mathrm{d}\theta }+\frac{m}{\sin \theta }\right){\Psi }^{m,\sigma }.$$

$${\mathcal{L}}^{m,\tilde{\nu }}\left({\Psi }^{m,\sigma }\right)=\frac{R_{\mathrm{p}}^2\sigma \tilde{\sigma }}{H_{\mathrm{oc}}}{\eta }^{m,\sigma },$$

$${\mathcal{L}}^{m,\tilde{\nu }}\equiv$$

$$-\frac{1}{1-{\tilde{\nu }}^2{\cos }^2\theta }\left(m\tilde{\nu }\frac{1+{\tilde{\nu }}^2{\cos }^2\theta }{1-{\tilde{\nu }}^2{\cos }^2\theta }+\frac{m^2}{{\sin }^2\theta }\right).$$

$$\begin{array}[]{ll}\displaystyle\Psi^{m,\sigma}=\sum_{n}\Psi_{n}^{m,\sigma}\Theta_{n}^{m,\tilde{\nu}}\left(\theta\right),&\displaystyle\eta^{m,\sigma}=\sum_{n}\eta_{n}^{m,\sigma}\Theta_{n}^{m,\tilde{\nu}}\left(\theta\right),\end{array}$$

$${\mathcal{L}}^{m,\tilde{\nu }}{\Theta }_n^{m,\tilde{\nu }}=-{\Lambda }_n^{m,\tilde{\nu }}{\Theta }_n^{m,\tilde{\nu }}.$$

Full text

*Aims. *We investigate the impact of eccentricity tides on the rotation of rocky planets hosting a thin uniform ocean and orbiting cool dwarf stars such as TRAPPIST-1, with orbital periods $\sim 1-10$ ${\rm days}$.

*Methods. *Combining the linear theory of oceanic tides in the shallow water approximation with the Andrade model for the solid part of the planet, we develop a global model including the coupling effects of ocean loading, self-attraction, and deformation of the solid regions. We derive from this model analytic solutions for the tidal Love numbers and torque exerted on the planet. These solutions are used with realistic values of parameters provided by advanced models of the internal structure and tidal oscillations of solid bodies to explore the parameter space both analytically and numerically.

*Results. *Our model allows us to fully characterise the frequency-resonant tidal response of the planet, and particularly the features of resonances associated with the oceanic tidal modes (eigenfrequencies, resulting maxima of the tidal torque and Love numbers) as functions of the planet parameters (mass, radius, Andrade parameters, ocean depth and Rayleigh drag frequency). Resonances associated with the oceanic tide decrease the critical eccentricity beyond which asynchronous rotation states distinct from the usual spin-orbit resonances can exist. We provide an estimation and scaling laws for this critical eccentricity, which is found to be lowered by roughly one order of magnitude, switching from $\sim 0.3$ to $\sim 0.06$ in typical cases and to $\sim 0.01$ in extremal ones.

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 the 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 preliminary constrain their rotation. As tides drive the long term evolution of planetary systems, this requires to characterise 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 well this configuration since its seven rocky planets orbit within a disk of radius 0.065 ${\rm 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 an important 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 by the semidiurnal tide solely 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 a 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 the Laplace’s tidal equation (e.g. Longuet-Higgins, 1968; Webb, 1980). In deeper oceans, stable stratification leads in addition 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 account in 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 atmospheric tides (e.g. Lindzen & Chapman, 1969), where the body is torqued away from synchronization (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 synchronization (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 idealized ocean planet with an eccentric orbit and no obliquity, the central object of the system being 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 into account self-consistently. We shall emphasize here that we do not include the coupling between the internal tidal heating of the planet and its structure and surface conditions, these laters being fixed in the model. Therefore, the planetary system is essentially parametrized 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 conclusions.

To describe the tidal response of the solid part, we use 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 inputs 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 use 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_{\rm AR}$ (’AR’ refers to ’Asynchronous Rotation’). As we focus on the frequency interval bounded by the $1{:}1$ (synchronization) and $3{:}2$ spin-orbit resonances, $e<e_{\rm AR}$ means that the planet is driven towards the spin-orbit synchronous rotation, while $e>e_{\rm 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 the 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 term 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 examine the simplified case of a single terrestrial planet orbiting its host star with an eccentric orbit. The planet, of mass $M_{\rm p}$ and radius $R_{\rm 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_{\rm oc}\ll R_{\rm p}$ and density $\rho_{\rm oc}$. Its orbital motion is described by its mean motion $n_{\star}$ and eccentricity $e$. The orbital angular momentum vector is denoted $\boldsymbol{n_{\star}}$. We introduce the planeto-centric referential $\mathcal{R}_{{\rm G}}{:}\left(O,\boldsymbol{X}_{{\rm G}},\boldsymbol{Y}_{{\rm G}},\boldsymbol{Z}_{{\rm G}}\right)$, where $\boldsymbol{Z}_{{\rm G}}\equiv\boldsymbol{n_{\star}}/n_{\star}$ (the symbol $\equiv$ meaning ’defined by’), $\boldsymbol{X}_{{\rm G}}$ and $\boldsymbol{Y}_{{\rm G}}$ designate the directions of two distant stars defining the orbital plane, and $O$ the planet gravity center.

We assume the absence of obliquity, which allows us to define the reference frame co-rotating with the planet as $\mathcal{R}_{{\rm p}}{:}\left(O,\boldsymbol{X}_{{\rm p}},\boldsymbol{Y}_{{\rm p}},\boldsymbol{Z}_{{\rm p}}\right)$ with $\boldsymbol{Z}_{{\rm p}}=\boldsymbol{Z}_{{\rm G}}$. The vectors $\boldsymbol{X}_{{\rm p}}$ and $\boldsymbol{Y}_{{\rm p}}$ define the equatorial plane of the planet, which is also its orbital plane in the present case. The rotating motion of $\mathcal{R}_{{\rm p}}$ with respect to $\mathcal{R}_{{\rm G}}$ is defined by the spin vector $\boldsymbol{\Omega}=\Omega\boldsymbol{Z}_{{\rm p}}$, $\Omega$ being the rotation rate of the planet. Hence, denoting $t$ the time, $\boldsymbol{X}_{{\rm p}}=\cos\left(\Omega t\right)\boldsymbol{X}_{{\rm G}}+\sin\left(\Omega t\right)\boldsymbol{Y}_{{\rm G}}$ and $\boldsymbol{Y}_{{\rm p}}=-\sin\left(\Omega t\right)\boldsymbol{X}_{{\rm G}}+\cos\left(\Omega t\right)\boldsymbol{Y}_{{\rm G}}$. We finally need to introduce the spherical coordinates, namely the radius $r$, colatitude $\theta$, and longitude $\varphi$, and the associated unit-vector basis $\left(\boldsymbol{{\rm e}}_{r},\boldsymbol{{\rm e}}_{\theta},\boldsymbol{{\rm e}}_{\varphi}\right)$.

The surface planet gravity is denoted $g$. We introduce the planet mean density $\rho_{{\rm p}}\equiv M_{\rm p}/\left(\frac{4\pi}{3}R_{\rm p}^{3}\right)$ and the mass of the ocean $M_{\rm oc}\equiv 4\pi R_{\rm p}^{2}H_{\rm oc}\rho_{\rm oc}$. Besides, the mass and radius of the solid part are denoted $M_{\rm c}\equiv M_{\rm p}-M_{\rm oc}\approx M_{\rm p}$ and $R_{\rm c}\equiv R_{\rm p}-H_{\rm oc}\approx R_{\rm p}$, respectively. To simplify calculations, we assume in the following that this layer can be treated as a homogeneous body of density $\rho_{\rm c}\equiv M_{\rm c}/\left(\frac{4\pi}{3}R_{\rm c}^{3}\right)\approx\rho_{{\rm p}}$.

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, similarly as the core and the mantle in rocky planets (e.g. Correia & Laskar, 2009).

In the framework of Rayleigh drag approximation chosen for this study (e.g. Vallis, 2006), the typical timescale of the core-ocean coupling by viscous friction is $\tau_{\rm R}\equiv\sigma_{\rm R}^{-1}$, where $\sigma_{\rm 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 $\sigma_{\rm R}=10^{-5}$ ${\rm s^{-1}}$ (Webb, 1980), $\tau_{\rm R}\approx 30$ hr, 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 enable thereby 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 classical tidal theory to describe analytically the dynamics of the planet tidal response, by treating successively the solid part and the oceanic layer.

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_{\star}$, that is $r_{\star}\gg R_{\rm p}$ typically, this gravitational potential is expressed in the accelerated frame of the planet as

[TABLE]

the first term of the right-hand side of the equation corresponding to the attraction by the host star, of mass $M_{\star}$, and the second term to the centrifugal force due to the orbital motion. We note that the symbol $\hat{\leavevmode\nobreak\ }$ is employed here and all along the article to highlight real quantities with respect to complex ones, the laters 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 $\tilde{\leavevmode\nobreak\ }$ 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_{\rm p}$),

[TABLE]

where we have removed the constant component as it does not contribute to the tidal force. The associated complex gravitational potential $U_{\rm T}$, such that $\hat{U}_{\rm T}=\real{U_{\rm T}}$, 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

[TABLE]

where $l$ and $m$ designate the latitudinal and longitudinal degrees, $s$ an integer, $\sigma_{m,s}=m\Omega-sn_{\star}$ the forcing tidal frequency of the mode associated with the doublet $\left(m,s\right)$, $P_{l}^{m}$ the normalised associated Legendre function associated with the doublet $\left(l,m\right)$ (see Appendix A), and $U_{{\rm T};l}^{m,s}$ the surface tidal gravitational potential associated with the triplet $\left(l,m,s\right)$. This later is given by

[TABLE]

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 we get different numerical factors). Let us introduce the Kronecker symbol $\delta_{l,k}$, such that $\delta_{l,k}=1$ if $l=k$ and $\delta_{l,k}=0$ otherwise. By analogy with the Kronecker symbol, we define by $\delta_{s<0}$ the coefficient such that $\delta_{s<0}=1$ if $s<0$ and [math] otherwise. The $A_{l,m,s}$ coefficients are thus expressed as

[TABLE]

In the above expression, the $\bar{P}_{l}^{m}$ are the unnormalised associated Legendre functions, and the eccentricity functions $X_{s}^{l,m}$ are the so-called Hansen 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 $U_{{\rm T};l}^{m,s}\propto\left(r/a\right)^{l}$, 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

[TABLE]

and expressed as

[TABLE]

with $s\in\mathbb{Z}$. In the case of a circular orbit ($e=0$), the tidal potential reduces to the semidiurnal component

[TABLE]

At small eccentricities ($e\ll 1$), eccentricity tides are predominantly forced by terms associated with $s=1$ (westward propagating potential) and $s=3$ (eastward propagating potential), given by (e.g. Ogilvie, 2014, Table 1)

[TABLE]

which are associated with the tidal frequencies $\sigma_{2,1}=2\Omega-n_{\star}$ and $\sigma_{2,3}=2\Omega-3n_{\star}$, respectively. At high eccentricities ($e\lesssim 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 $U_{{\rm T};2}^{2,3}\approx U_{{\rm T};2}^{2,2}$, that is at the transition eccentricity $e_{\rm trans}\approx 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\lesssim e_{\rm trans}$, meaning that the critical eccentricity beyond which asynchronous final rotation states of equilibrium can exist is $e_{\rm AR}=e_{\rm trans}$ in this case. As shown in the following, $e_{\rm 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, are traveling 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 $n_{\star}\leq\Omega<\left(3/2\right)n_{\star}$. 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 ($\Omega\leq n_{\star}$), the eccentricity component of degree $s=3$ torques the planet towards synchronisation, as the semidiurnal one does. The desynchronising component is that of degree $s=1$ in this case, from the moment that $\Omega>\left(1/2\right)n_{\star}$. The role played the two predominating eccentricity components in the approximation of small eccentricities may be generalised to other ones ($s\in\mathbb{Z}$). The degree-$s$ eccentricity component tends to drive the planet away from synchronisation provided that $\Omega$ satisfies the conditions

[TABLE]

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.

3.2 Tidal response of the solid part

As a first step, let us 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 as 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 to preliminary assume 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 $\mu$. This body undergoes a time-periodic tidal force of frequency $\sigma$. Thus, in the linear approximation, any quantity $\hat{q}$ describing the distortion can be written as $\hat{q}=\real{q^{\sigma}{\rm e}^{i\sigma t}}$, where $q^{\sigma}$ stands for the associated complex quantity, which depends on the forcing frequency. We introduce the stress tensor $\underline{\underline{\boldsymbol{\zeta}}}$ and the strain tensor in the linear approximation

[TABLE]

the notation $\boldsymbol{\xi}$ designating the displacement vector. The rheology of the solid is determined by the relationship between $\underline{\underline{\boldsymbol{\zeta}}}$ and $\underline{\underline{\boldsymbol{\varepsilon}}}$.

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

[TABLE]

where $\tilde{\mu}$ is the complex shear modulus accounting for the material rheology. We note that the imaginary part of $\tilde{\mu}$ 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, $\tilde{\mu}$ 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 $\tilde{\mu}$ as an explicit function of the forcing frequency,

[TABLE]

where we have introduced the gamma function $\Gamma$ (Abramowitz & Stegun, 1972), the dimensionless rheological exponent $\alpha$ (the values of $\alpha$ are determined experimentally and typically fall within the interval $0.2-0.4$ for olivine-rich rocks, see Castelnau et al., 2008), the Maxwell relaxation time $\tau_{\rm M}$ associated with the viscous response of the material ($\tau_{\rm M}$ is defined by the ratio of viscosity $\eta$ to unrelaxed rigidity, $\tau_{\rm M}\equiv\eta/\mu$, see e.g. Castillo-Rogez et al., 2011), and the Andrade – or anelastic – time $\tau_{\rm A}$ related to anelasticity. This latter parameter is the characteristic timescale associated with the material creep.

By setting $\tau_{\rm A}=+\infty$, 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_{{\rm T}}$ and the Maxwell time. In the zero-frequency limit ($P_{{\rm T}}\gg\tau_{\rm M}$), the response is dominated by viscous friction, leading to the purely imaginary complex shear modulus $\tilde{\mu}\sim i\mu\tau_{\rm M}\sigma$ (we remind here that $\sigma$ designates the forcing frequency of the tidal mode). In the high-frequency regime ($P_{{\rm T}}\ll\tau_{\rm M}$), the behaviour of the solid part tends to be purely elastic, with $\tilde{\mu}\sim\mu$. In the case of telluric planets, the order of magnitude of the Maxwell time is of several centuries (typically $\tau_{\rm M}\sim 10^{4}-10^{5}$  ${\rm days}$; for the Earth $\tau_{\rm M}$ is about 500 years, see e.g. Efroimsky, 2012), that is far greater than usual tidal periods ($P_{{\rm T}}\sim 10^{0}-10^{2}$  ${\rm days}$). As a consequence, the imaginary part of $\tilde{\mu}$ drops rapidly as $\sigma$ increases in the high-frequency regime, which leads to underestimate the energy 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 tidally dissipated within planetary interiors with a minimal number of additional parameters, namely $\alpha$ and $\tau_{\rm 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 $\alpha$ (this is illustrated by Fig. 3, discussed further, 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 $\tau_{\rm A}$ and $\tau_{\rm M}$ to be comparable in the general case since they are related to two different physical mechanisms. The order of magnitude of $\tau_{\rm 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 $\beta=\mu^{-1}\tau_{\rm A}^{-\alpha}$.

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 – detailed in Sotin et al. (2007) – that solves numerically 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 as 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 this later 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 homogeneous in density presenting 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. 2019, in preparation).

The variation of mass distribution associated with the tidal distortion is described by Love numbers. In the linear approach, each degree-$l$ mode of the expansion in spherical harmonics (see Eq. 3) has its associated Love numbers, denoted by $k_{l}^{\sigma}$, $h_{l}^{\sigma}$, $\mathfrak{k}_{l}^{\sigma}$, and $\mathfrak{h}_{l}^{\sigma}$ in the present work. The first two parameters, $k_{l}^{\sigma}$ and $h_{l}^{\sigma}$, 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 $\mathfrak{k}_{l}^{\sigma}$ and $\mathfrak{h}_{l}^{\sigma}$ 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 are given by (e.g. Munk & MacDonald, 1960)

[TABLE]

where we have introduced the dimensionless effective rigidity

[TABLE]

with

[TABLE]

The contribution of the degree-$l$ component to the tidal torque exerted on the solid part is proportional to $\Im\left\{k_{l}^{\sigma}\right\}$ (e.g. Makarov, 2012). For comparison, the imaginary part of the quadrupole Love number $k_{2}^{\sigma}$ 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 $\sigma\tau_{\rm M}$. The behaviour of the two models does not differ in the low-frequency regime, where $\Im\left\{k_{2}^{\sigma}\right\}$ increases with $\sigma$ linearly. A maximum is reached for $P_{{\rm T}}\sim\tau_{\rm M}$. In the high-frequency regime, the torque decays as $\sigma$ 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 ($\alpha=0.25$), and $-1$ per decade for the Maxwell model (Fig. 3).

3.3 Tidal response of the ocean

We now consider the tidal response of the ocean. Similarly as 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 $\xi_{\rm 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 supposed incompressible. Mean flows are ignored (the whole planet rotates as a solid, 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 $\boldsymbol{V}=V_{\theta}\boldsymbol{{\rm e}}_{\theta}+V_{\varphi}\boldsymbol{{\rm e}}_{\varphi}$ and the variation of the oceanic depth $\eta=\xi_{\rm oc}-\xi_{\rm c}$, the notation $\xi_{\rm oc}$ designating the oceanic surface displacement (see Fig. 1). In the frame rotating with the planet ($\mathcal{R}_{{\rm p}}$), the dynamics of the oceanic tidal response is governed by the momentum equation (Matsuyama, 2014),

[TABLE]

and the equation of mass conservation,

[TABLE]

where we have introduced the horizontal gradient operator in spherical coordinates,

[TABLE]

the horizontal component of the velocity divergence,

[TABLE]

the Rayleigh drag frequency $\sigma_{\rm 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 $\Psi=\Psi\left(t,\theta,\varphi\right)$ and defined further, which encompasses both the tidal gravitational forcing and coupling with the solid part. We note that the radial component of the Coriolis acceleration $-2\Omega V_{\varphi}\boldsymbol{{\rm 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 $2\boldsymbol{\Omega}\wedge\boldsymbol{V}=2\Omega\cos\theta\left(V_{\theta}\boldsymbol{{\rm e}}_{\varphi}-V_{\varphi}\boldsymbol{{\rm e}}_{\theta}\right)$.

In the linear approximation, the response of a given mode is proportional to the corresponding component of the forcing. It follows that $\Psi$, $\boldsymbol{V}$, and $\eta$ can be expanded in Fourier series of $t$ and $\varphi$ similarly as the tidal gravitational potential (Eq. 3),

[TABLE]

The Fourier coefficients $\Psi^{m,\sigma}$, $\boldsymbol{V}^{m,\sigma}$, and $\eta^{m,\sigma}$ may themselves be expanded as series of the normalised associated Legendre functions $P_{l}^{m}$,

[TABLE]

where $\Psi_{l}^{m,\sigma}$, $\boldsymbol{V}_{l}^{m,\sigma}$ and $\eta_{l}^{m,\sigma}$ designate the degree-$l$ components of the expansion.

By introducing the complex tidal frequency $\tilde{\sigma}$ and spin parameter $\tilde{\nu}$, defined by

[TABLE]

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 $\Psi$,

[TABLE]

Then, substituting $V_{\theta}^{m,\sigma}$ and $V_{\varphi}^{m,\sigma}$ by the above expressions in Eq. 19, we end up with

[TABLE]

the notation $\mathcal{L}^{{m,\tilde{\nu}}}$ referring to the Laplace’s tidal operator, defined by (e.g. Lee & Saio, 1997)

[TABLE]

We recognise in Eq. 27 the Laplace’s Tidal Equation, which determines the horizontal structure of the fluid tidal response (e.g. Lindzen & Chapman, 1969). If we assume that the Fourier components $\Psi^{m,\sigma}$ and $\eta^{m,\sigma}$ may be written as

[TABLE]

the parameter $n$ being an integer, and consider that the latitudinal functions $\Theta_{n}^{m,\tilde{\nu}}$ are bounded at the poles, Eq. 27 defines an eigenfunction-eigenvalue problem. Hence, the set of eigenfunctions $\left\{\Theta_{n}^{m,\tilde{\nu}}\right\}$, called Hough functions after Hough’s work (Hough, 1898), is associated with a set of eigenvalues $\left\{\Lambda_{n}^{m,\tilde{\nu}}\right\}$ through the relationship

[TABLE]

One shall also introduce here the corresponding complex equivalent depths $h_{n}^{m,\sigma}$, which are defined by analogy with real equivalent depths (e.g. Taylor, 1936; Lindzen & Chapman, 1969), by

[TABLE]

The Hough function $\Theta_{n}^{m,\tilde{\nu}}$ may be written as a combination of the normalised associated Legendre functions,

[TABLE]

Conversely, the normalised associated Legendre function $P_{l}^{m}$ may be written as

[TABLE]

In the preceding expressions, the $A_{n,l}^{m,\tilde{\nu}}$ and $B_{l,n}^{m,\tilde{\nu}}$ are complex overlap coefficients. They are computed in the meantime as the eigenvalues $\Lambda_{n}^{m,\tilde{\nu}}$ using the standard method detailed by Wang et al. (2016), which is based on series of the normalised associated Legendre functions. Besides, we introduce here the overlap coefficients

[TABLE]

which will be used further to weight the degree-$n$ components of the oceanic tidal torque.

In the absence of friction ($\sigma_{\rm R}=0$), $h_{n}^{m,\sigma}=H_{\rm oc}$ 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 here into details. We just recall the main aspects of this regime.

When $\sigma_{\rm R}=0$, the complex quantities defined by Eq. 24 reduce to $\tilde{\sigma}=\sigma$ and $\tilde{\nu}=2\Omega/\sigma$, which are the real forcing frequency and spin parameter usually met in literature. The set of Hough functions divides into two families. A first family, refereed to as gravity modes or g-modes (see e.g. Lee & Saio, 1997), is defined for $\tilde{\nu}\in\mathbb{R}$. This family corresponds to the ordinary spherical modes modulated by rotation and naturally reduce to the set of the associated Legendre functions in the non-rotating case. For $\left|\tilde{\nu}\right|>1$, an other 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 confined by 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 $\left|\sigma\right|\gg\sigma_{\rm R}$. For $\left|\sigma\right|\lesssim\sigma_{\rm R}$, the friction affects the behaviour of tidal modes. While the ratios $\left|\sigma\right|/\sigma_{\rm R}$ and $\left|2\Omega\right|/\sigma_{\rm R}$ 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 ($\left|\sigma\right|/\sigma_{\rm R}\rightarrow 0$), 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 $\left|\tilde{\nu}\right|\rightarrow 0$ (i.e. $\Omega\rightarrow 0$ or $\sigma\rightarrow+\infty$). This is no longer the case when a drag is taken into account.

We now come back to $\Psi$, which still has 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 $\gamma_{U;l}^{\sigma}$ and $\gamma_{\xi;l}^{\sigma}$, respectively. In the framework of the thin shell approximation111In the case of a thick oceanic layer, the ocean depth intervenes in the expressions of tilt factors as well as in those of the solid Love numbers (Dermott, 1979; Remus et al., 2012)., these factors are defined by

[TABLE]

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 $\Psi^{m,\sigma}$ is thus defined as

[TABLE]

By expanding $\Psi^{m,\sigma}$ and $\eta^{m,\sigma}$ in series of Hough functions, we transform Eq. 27 into

[TABLE]

the $\Psi_{n}^{m,\sigma}$ being given by

[TABLE]

In practice, the series of Hough functions and associated Legendre functions given by Eqs. 23 and 29 are truncated in numerical calculations. We denote by $N$ the number of functions used to approximate series, which is arbitrarily chosen large enough to be associated with negligible overlap coefficients. It follows that the $\eta_{n}^{m,\sigma}$ of Eq. 29 are the solutions of an algebraic system of the form

[TABLE]

In Eq. 40, the $\sigma_{n,k}$ are the complex characteristic frequencies defined by

[TABLE]

where $k_{{\rm h};n}\equiv\sqrt{\Lambda_{n}^{m,\tilde{\nu}}}/R_{\rm p}$ designates the horizontal wavenumber of the degree-$n$ mode. The components $F_{n}^{m,\sigma}$ of the force vector in Eq. 40 are defined by

[TABLE]

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_{\rm p}\ll 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 $\eta_{n}^{m,\sigma}$ by using the overlap coefficients introduced in Eq. 32,

[TABLE]

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. This is the object of this section. The obtained analytic formulae are used further to 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_{\rm p}\ll a$), the tidal torque is thus given by

[TABLE]

where $k_{{\rm p};2}^{2,\sigma}$ is the quadrupolar component of the effective gravitational Love number of the planet, $k_{{\rm p};l}^{m,\sigma}$, defined in the thin shell approximation by

[TABLE]

In this expression, the degree-$l$ component of the tidal gravitational potential $U_{{\rm T};l}^{m,\sigma}$ and variation of ocean thickness $\eta_{l}^{m,\sigma}$ are given by Eqs. 4 and 43, respectively.

4.2 Pure solid tidal response

In the absence of oceanic layer, $R_{\rm c}=R_{\rm 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)

[TABLE]

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 – for instance – Makarov (2012) (Eq. 6). As discussed in Sect. 3.2, the Andrade and Maxwell times generally far exceed typical tidal periods. As a consequence the imaginary part of $k_{2}^{\sigma}$ may be approximated by

[TABLE]

with

[TABLE]

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 ($\mu\rightarrow+\infty$), so that the Love numbers defined by Eq. 15 annihilate. Besides, the Cowling approximation (e.g. Cowling, 1941; Unno et al., 1989) is assumed, which means that the term in $\rho_{\rm oc}/\rho_{\rm c}$ resulting from the self-attraction of the ocean in $\Psi$ is ignored. In this case, the oceanic tidal torque may be written similarly as that of the solid part (Eq. 46),

[TABLE]

where $k_{{\rm oc};2}^{2,\sigma}$ designates the degree-2 tidal gravitational Love number of the ocean.

Owing to the Cowling approximation, $\sigma_{n,k}=0$ for $n\neq 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),

[TABLE]

the notation $\tilde{\sigma}_{n}^{-}$ and $\tilde{\sigma}_{n}^{+}$ designating the complex eigenfrequencies of the degree-$n$ Hough mode, defined by (see e.g. Webb, 1980, Eq. (2.12))

[TABLE]

In this case, the quadrupolar Love number introduced in Eq. 49 is simply expressed as

[TABLE]

Using the mean density of the planet ($\rho_{{\rm p}}=3M_{\rm p}/\left(4\pi R_{\rm p}^{3}\right)\approx\rho_{\rm c}$), it may also be put into the form

[TABLE]

which highlights the fact that a resonance occurs when the equivalent depth of a given mode $h_{n}^{m,\sigma}=R_{\rm p}^{2}\sigma^{2}/\left(\Lambda_{n}^{m,\tilde{\nu}}g\right)$ equalises the ocean depth in the absence of friction.

The friction of tidal flows with the oceanic floor affects both the overlap coefficients $C_{2,n,2}^{2,\tilde{\nu}}$ and the equivalent depths of Hough modes. Since these parameters are complex in the general case, characterising the dependence of the tidal torque on $\sigma_{\rm R}$ is not straightforward. However, this may be done in the quasi-adiabatic asymptotic regime ($\sigma_{\rm R}\ll\left|\sigma\right|$), where the imaginary part of $C_{2,n,2}^{2,\tilde{\nu}}$ is negligible, and for g-modes, as these laters are associated with positive $\Lambda_{n}^{m,\tilde{\nu}}$.

By introducing the characteristic frequency of the degree-$n$ surface gravity wave,

[TABLE]

we express the imaginary part of the oceanic Love number as

[TABLE]

The preceding expression characterises the shape of a resonant peak in the frequency-spectra of the tidal torque. The frequency at which the resonant peak of a mode reaches a maximum is expressed as

[TABLE]

In the quasi-adiabatic regime ($\sigma_{\rm R}\ll\left|\sigma\right|$), it reduces to

[TABLE]

The corresponding maximum of $\Im\left\{k_{{\rm oc};2}^{2,\sigma}\right\}$ for the resonance associated with the degree-$n$ Hough mode is denoted by $\left.\Im\left\{k_{{\rm oc};2}^{2,\sigma}\right\}\right|_{{\rm max};n}$ and obtained by substituting $\sigma$ by Eq. 57 in Eq. 55. If the contributions of non-resonant components are neglected, we thus end up with

[TABLE]

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 $\left.\Im\left\{k_{{\rm oc};2}^{2,\sigma}\right\}\right|_{{\rm max};n}\propto\sigma_{n}\propto\sqrt{H_{\rm oc}}$ 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, $\sigma_{2,2}=2\left(\Omega-n_{\star}\right)$, in the super-synchronous regime ($\Omega>n_{\star}$). In this figure, the tidal response of the solid part is ignored to focus on the oceanic tidal torque, which is reduced to the contribution of the degree-0 Hough mode for simplification. For a pedagogical purpose, 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 $\pm\sigma_{0}$ and the resulting tidal torque contributions, which both drive the system towards synchronism, are represented schematically by blue cones in Fig. 4. As Eq. 54 shows, the frequencies of these intrinsic modes of the ocean increase with the ocean depth $H_{\rm oc}$. Indeed, we can see in Fig. 4 that the blue cones’ positions shift towards the higher frequencies when $H_{\rm oc}$ increases from panels a to f.

One of these peaks is located in the super-synchronous frequency range ($\sigma_{2,2}>0$) and the other one in the sub-synchronous frequency range ($\sigma_{2,2}<0$), symmetrically (this symmetry with respect to the synchronisation 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 ($\sigma_{2,s}=-\sigma_{0}$). This configuration corresponds for instance to case 1 of Fig. 2, where the degree-3 eccentricity potential is propagating eastward.

In this framework, Fig. 4 simultaneously shows the forcing frequencies $\sigma_{2,s}=2\Omega-sn_{\star}$ 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<\sigma_{2,2}<n_{\star}$. These particular frequencies are superscripted by the symbol $\blacktriangleup$, and the created peak is designated by a red cone. Hence, the contribution of one of the eccentricity components is potentially leading to asynchronous equilibrium states located at $\sigma_{2,2}\approx\sigma_{2,2}^{\blacktriangleup}$. 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 $\sigma_{2,2}$.

For the cases a to c of Fig. 4, the assumption made about $\sigma_{2,2}^{\blacktriangleup}$ (i.e. $0<\sigma_{2,2}^{\blacktriangleup}<n_{\star}$) 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 $\sigma_{2,3}=-\sigma_{0}=\sigma_{2,3}^{\blacktriangleup}$. When the depth of the ocean increases (from d to f), keeping $\sigma_{2,2}^{\blacktriangleup}$ in the range 0 - $n_{\star}$ leads to an excitation of the ocean resonant mode only with the degree-4 eccentricity mode (for $\sigma_{2,4}=-\sigma_{0}=\sigma_{2,4}^{\blacktriangleup}$). 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 $\sigma_{0}=n_{\star}/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 $U_{{\rm T};2}^{2,3}\ll U_{{\rm T};2}^{2,2}$. As a consequence, the oceanic tide leads the planet toward synchronisation ($\Omega=n_{\star}$) in this configuration, and cannot generate an asynchronous rotation state of equilibrium in the range $-n_{\star}<\sigma_{2,2}<n_{\star}$.

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 $\propto\sqrt{H_{\rm oc}}$ (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 – gets also closer to synchronisation until being annihilated by the solid tidal torque, which outweighs the oceanic torque when $P_{{\rm T}}\sim\tau_{\rm 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 is 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\gtrsim e_{\rm trans}$, non-linearities lead eccentricity components of higher degrees to predominate; see e.g. Ogilvie, 2014, Fig. 3).

Similarly as the peak created by the degree-3 component, the peak associated with the degree-4 eccentricity term moves towards the synchronisation while the ocean depth increases, as shown by configurations (e) and (f). When $\sigma_{0}$ becomes greater than $2n_{\star}$, 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 $\Omega_{{\rm eq};s}\approx\left(s/2\right)n_{\star}$, as discussed in Sect. 5.

In reality, the simplified oceanic tidal response is not solely due to the degree-0 Hough mode of frequency $\sigma_{0}$ since $\left|\Omega\right|\sim n_{\star}\sim\left|\sigma_{2,s}\right|$ 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), these laters 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 $\sigma_{\rm R}\rightarrow 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 to $\sigma_{\rm R}\approx 10^{-5}$ ${\rm s^{-1}}$ (e.g. Webb, 1980), which is comparable with typical tidal frequencies. In this case, the enhanced oceanic tidal response resulting from a resonance remains in the framework of the linear theory, and the associated amplification of the tidal torque generally does not exceed a decade in logarithmic scale (see 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^{-5}$ ${\rm s^{-1}}$ is an upper estimation for $\sigma_{\rm R}$. It seems likely that $\sigma_{\rm R}$ takes smaller values in the case of a planet hosting a uniform ocean. To take these discrepancies into account, we consider the values $10^{-7}$ and $10^{-6}$ ${\rm s^{-1}}$ in addition to $10^{-5}$ ${\rm s^{-1}}$ 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.

4.4 Frequency behaviour of the tidal torque

The tidal torque resulting from the planet tidal response is characterised by a complex frequency behaviour, which slightly diverges from the idealised picture showed 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_{\star}\sim 1-10$ days) since this regime corresponds to $H_{\rm oc}\sim 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 days orbital period (Gillon et al., 2017).

We thus consider along the whole study the case of an Earth-sized planet orbiting TRAPPIST-1 ($M_{\star}=0.09\,M_{\sun}$, see Van Grootel et al., 2018) with a 6-days 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 purpose, 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<\sigma_{2,2}<n_{\star}$. We set the typical frequency of the Rayleigh friction to $\sigma_{\rm R}=1.0\times 10^{-6}$ ${\rm s^{-1}}$. Finally, the response of the solid part is described by the values of Andrade parameters given by Table 1, as mentioned in Sect. 3.2.

The tidal torque exerted on the planet is plotted in Fig. 5 as a function of the normalised semidiurnal tidal frequency $\sigma_{2,2}/\left(2n_{\star}\right)=\left(\Omega-n_{\star}\right)/n_{\star}$, 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_{\rm oc}=0$), the case of a pure oceanic tide ($\mu=+\infty$), and the general case, respectively. In the three cases, the tidal torque is calculated using the expression given by Eq. 44, with $H_{\rm oc}=0$ or $\mu=+\infty$ to obtain the pure solid and oceanic response, respectively. To highlight the effect of eccentricity components, regimes where the planet is driven towards (away from) the spin-orbit synchronous rotation state of equilibrium are designated by solid (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 $\mathcal{T}_{\rm sol}\propto\left|\sigma_{2,2}\right|^{-\alpha}$ 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 $\Omega_{{\rm eq};s}\approx\left(s/2\right)n_{\star}$ (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 $\left(\Omega-n_{\star}\right)/n_{\star}\approx 0.4$ clearly corresponds to the configuration illustrated by Fig. 4 (panel (b)), 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 ($\sigma_{\rm R}\lesssim\left|\sigma\right|<\left|2\Omega\right|$), the tidal forcing tends to couple with an infinite number of Hough modes as $\left|\tilde{\nu}\right|$ increases (this corresponds to the regime of sub-inertial waves, defined by $\left|\tilde{\nu}\right|>1$, and discussed by, e.g. Lee & Saio, 1997). The method adopted to solve the Laplace’s tidal equation (see Eq. 30) in this study cannot treat correctly 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 $\left|\tilde{\nu}\right|$ increases.

As a result, the oceanic tidal response is poorly described by the theory for $\left|\tilde{\nu}\right|\gg 1$, 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 to explore 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 $\left|2\Omega\right|\lesssim\sigma_{\rm R}$, which is the case in this study if $\sigma_{\rm R}\gtrsim 10^{-6}$ ${\rm s^{-1}}$.

We finally consider the tidal torque in the general case, where the tidal responses of the solid and oceanic layers are coupled together by gravitational and pressure forces (Fig. 5, bottom panel). We mainly retrieve here 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 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.

5 Analytical estimation of the 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\gtrsim e_{\rm trans}\sim 0.25$ typically, as showed in Sect. 3. The critical eccentricity $e_{\rm 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 as 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\ll 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\geq 3$, as 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

[TABLE]

where $\sigma_{2,s}=2\Omega-sn_{\star}$ is the degree-$s$ eccentricity frequency defined by Eq. 6 and $\sigma_{n}$ the eigenfrequency of the degree-$n$ Hough mode given by Eq. 54. By combining Eqs. 11 and 59, and introducing the notation $\tilde{\nu}_{2,s}\equiv\tilde{\nu}\left(\sigma_{2,s}\right)$, we hence derive from (i) and (ii) a first condition on the ocean depth ($H_{\rm oc}$), which must satisfy the inequality

[TABLE]

This inequality determines the range of depths for which a resonance may occur. We observe that the range of $H_{\rm 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_{\rm oc}$ increases with $n_{\star}$, which favours resonances in the case of close-in terrestrial planets. Conversely, this limit is inversely proportional to the eigenvalue $\Lambda_{n}^{m,\tilde{\nu}}$ of the mode. In the quasi-adiabatic regime, the $\Lambda_{n}^{m,\tilde{\nu}}$ are sorted in ascending order (e.g. Lee & Saio, 1997). Thus, the upper limit of $H_{\rm oc}$ is lower for high degrees $n$ than for low degrees. Basically, this means that the resonance of the degree-[math] mode will express 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 the planet self-attraction (this is the Cowling approximation mentioned above). Moreover, as we focus on the case of small eccentricities ($e\ll 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_{\star}/2$ in the regime where $P_{{\rm T}}\ll\tau_{\rm 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 $\Omega_{{\rm eq};s}\approx\left(s/2\right)n_{\star}$. For small parameters $\alpha$ of the Andrade rheology (typically $\alpha=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_{\star}$) 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

[TABLE]

where $\sigma_{2,s}=2\Omega-sn_{\star}=-\sigma_{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 $\sigma_{2,2}=2\left(\Omega-n_{\star}\right)$ 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 ($\mu\rightarrow+\infty$). In this case, $k_{{\rm p};2}^{2,\sigma}\approx k_{{\rm oc};2}^{2,\sigma}$. The associated Legendre function $P_{2}^{2}$ is generally well coupled with the degree-0 Hough mode, which means that $C_{2,n,2}^{2,\tilde{\nu}}\ll 1$ 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

[TABLE]

To obtain this expression, we have assumed that the degree-$s$ eccentricity component is resonant, that is $\sigma_{2,s}=-\sigma_{0}$ and $\sigma_{2,2}=\left(s-2\right)n_{\star}-\sigma_{0}$ (Eq. 59), made use of Eqs. 55 and 58 for the degree-0 Hough mode, and introduced the notation $\tilde{\nu}_{2,2}\equiv\tilde{\nu}\left(\sigma_{2,2}\right)$.

The condition on the eigenfrequency $\sigma_{0}$ resulting from the preceding expression is implicit in the general case since the overlap coefficients ($C_{2,0,2}^{2,\tilde{\nu}_{2,2}}$ and $C_{2,0,2}^{2,\tilde{\nu}_{2,s}}$) and eigenvalue ($\Lambda_{0}^{2,\tilde{\nu}_{2,s}}$) associated with the degree-0 Hough mode both depend on the forcing frequencies $\sigma_{2,s}$ and $\sigma_{2,2}$. This dependence can be neglected in the regime of rapid rotation ($\left|\Omega\right|\gg n_{\star}$) since spin parameters determining the solutions of Laplace’s tidal equation (Eq. 30) hardly vary in this case. In the general case ($\left|\Omega\right|\sim n_{\star}$), 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 to 8.

Nevertheless, the fact that analytic solutions of the Laplace’s tidal equation do not exist in the general case leads us to ignore this dependence here, and we obtain thereby that condition (iii) is satisfied for

[TABLE]

or

[TABLE]

where we have introduced the supposed constant parameter

[TABLE]

Eqs. 63 and 64 determine the interval of $H_{\rm 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 $\gamma$, which compares the Rayleigh friction frequency to a characteristic frequency proportional to the degree-$s$ eccentricity component of the tidal gravitational potential. By reminding us Eq. 60, we hence observe that the degree-$s$ eccentricity term cannot drive the planet away from spin-orbit synchronous rotation if $\gamma_{s}\gg 1$ since the torque induced by the semidiurnal component is systematically stronger.

Asynchronous rotation rates may exist only if $\gamma_{s}\lesssim 1$. In the asymptotic limit ($\gamma_{s}\rightarrow 0$), the eccentricity component always predominates when its forcing frequency is equal to the eigenfrequency of the degree-[math] Hough mode ($\sigma_{0}$). In this case, the only condition that has to be satisfied to make possible the existence of asynchronous final rotation rates is the condition given by Eq. 60. This highlights the criterion for which resonances may enable eccentricity terms to predominate at small eccentricities,

[TABLE]

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,

[TABLE]

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_{\star}=6$ ${\rm days}$, $n_{\star}=1.21\times 10^{-5}$ ${\rm s^{-1}}$. This leads to $e_{\rm AR}\approx 0.1$ for $\sigma_{\rm R}=10^{-5}$ ${\rm s^{-1}}$ and $e_{\rm AR}\approx 0.01$ for $\sigma_{\rm R}=10^{-6}$ ${\rm s^{-1}}$. We retrieve these orders of magnitude in the numerical results detailed in Sect. 6.

We investigated above the case where the oceanic tidal torque always predominates over that 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

[TABLE]

By assuming that $\left|\sigma_{2,2}\right|\gg\max\left\{\tau_{\rm M}^{-1},\tau_{\rm A}^{-1}\right\}$, 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

[TABLE]

where we have introduced the normalised frequency $X=\sigma_{0}/n_{\star}$ and the dimensionless constant

[TABLE]

Similarly as what we did to establish Eq. 67, we suppose that the oceanic tidal response behaves as in the non-rotating case, and ignore the implicit dependence of the overlap coefficient $C_{2,0,2}^{2,\tilde{\nu}_{2,s}}$ on the spin rotation ($\Omega$) and forcing ($\sigma_{2,s}$) frequencies through the spin parameter ($\tilde{\nu}_{2,s}$). We thus assume that $C_{2,0,2}^{2,\tilde{\nu}_{2,s}}\approx 1$. 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 $\alpha$ is real. However, the calculation of the maximum of the polynomial function corresponding to the left-hand member of the equation is straightforward. Denoting by $u_{\alpha}$ this maximum, we obtain for $s=3$,

[TABLE]

Since the inequality given by Eq. 69 can be satisfied only if $u_{\alpha}>\mathcal{A}$, it follows that an estimation of the critical eccentricity can be derived by substituting $X_{3}^{-3,2}=\left(7/2\right)e$ in this later inequality. We end up with

[TABLE]

The comparison between the two obtained expressions of $e_{\rm 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_{\rm AR}\,\propto\,\sigma_{\rm R}$. When the solid tide determines the non-resonant background level, this scaling law switches to $e_{\rm AR}\,\propto\,\sigma_{\rm R}^{1/2}$. As the non-resonant background level of the oceanic component of the tidal torque scales as $\,\propto\,\sigma_{\rm R}$ (see e.g. Auclair-Desrotour et al., 2018, Fig. 5), this means that the critical eccentricity ceases to decay proportionally with $\sigma_{\rm R}$ from the moment that $\sigma_{\rm 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 the planet final rotation with the 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_{\Earth}$ super-Earth orbiting the TRAPPIST-1 dwarf star ($M_{\star}=0.09\ M_{\sun}$). In the first case, we consider two orbital periods: $P_{\star}=1$ day and $P_{\star}=6$ days. In the second case, the super-Earth orbits the star with a 6 days orbital period. As the Rayleigh drag frequency cannot be accurately specified for an exoplanet of unknown topography, we perform the calculations for $\sigma_{\rm R}=10^{-7},10^{-6},10^{-5}$ ${\rm s^{-1}}$, the lower and upper values corresponding to weak and strong drags, respectively. The values of Andrade parameters used for the two planets are those specified in Table 1.

For given eccentricity and ocean depth, the final rotation of a planet is determined by successive iterations from an initial value, specified thereafter, with a frequency step scaling as $\,\propto\,\sigma_{\rm R}$ in order to take the dependence of the width of resonance peaks on $\sigma_{\rm R}$ into account. 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 of $\mathcal{T}_{\rm p}$. 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<\sigma_{2,2}<n_{\star}$ (see Fig. 4), which corresponds to the super-synchronous frequency range ($\Omega>n_{\star}$). In this interval, the tidal torque generated by the solid part in the absence of rotation reaches its maximal values for $\left|\sigma_{2,2}\right|\sim\tau_{\rm M}^{-1}$ and $\left|\sigma_{2,3}\right|\sim\tau_{\rm M}^{-1}$, that is in the vicinity of the interval bounds. As a consequence, if the initial value of the semidiurnal tidal frequency is too close from 0 or $n_{\star}$, 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, $\Omega_{{\rm eq};3}\approx\left(3/2\right)n_{\star}$ (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 $\sigma_{2,2}=n_{\star}/2$. This corresponds to the initial rotation rate $\Omega=\left(5/4\right)n_{\star}$.

Results are plotted in Fig. 6. Each panel of the figure represents the normalised final rotation period $P_{\rm rot}/P_{\star}$ 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_{\star}$), and the Rayleigh drag frequency ($\sigma_{\rm R}$). We note that, although global oceans of depths $\lesssim 0.1$ km are very unlikely owing to the extremely flat topography it would suppose, the lower bound of ${\log_{10}}\left(H_{\rm oc}\right)$ is set to ${\log_{10}}\left(H_{\rm oc}\right)=-3$ 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 $n_{\star}<\Omega_{{\rm eq}}<\left(3/2\right)n_{\star}$. The six top and middle panel correspond to the Earth, the three bottom ones to the 10-$M_{\Earth}$ super-Earth.

Numerical calculations highlight the mechanism discussed in the previous sections: a resonance associated with an oceanic mode can decrease the lower eccentricity of the region where asynchronous final states exist. This region approximately corresponds to ${\log_{10}}e\gtrsim-0.6$ (i.e. $e\approx 0.25$) in absence of ocean, that is in the asymptotic limit $H_{\rm oc}\rightarrow 0$. We note that this value is approximately the transition eccentricity $e_{\rm 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_{{\rm oc};s}$. Moreover, the asynchronous rotation state of equilibrium resulting from a resonance tends to get closer to the synchronisation while $H_{\rm oc}$ increases, since it follows the decay of the resonant eigenfrequency (see Fig. 4).

Owing to the scaling law $H_{{\rm oc};s}\,\propto\,\left(s-2\right)^{2}n_{\star}^{2}$ given by Eq. 60, the resonant ocean depths increase when the orbital period switches from 6 to 1 day. 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_{\rm AR}\approx 0.25$ to $e_{\rm AR}\approx 0.015$ (see e.g. the top right and middle panels of Fig. 6). In the most extremal treated case ($P_{\star}=1$ ${\rm day}$ and $\sigma_{\rm R}=10^{-7}$ ${\rm s^{-1}}$, middle right panel), the lowest critical eccentricity reached is $e_{\rm AR}\approx 0.006$. However, this later case is unlikely since it clearly goes beyond the scope of the linear theory.

By considering the cases defined by a 6-days orbital period and the Rayleigh frequencies $10^{-5}$ and $10^{-6}$ ${\rm s^{-1}}$ (top left and middle panels), we observe the transition from the frictional regime ($\sigma_{\rm R}\sim\left|\sigma\right|$) to the quasi-adiabatic regime ($\sigma_{\rm R}\ll\left|\sigma\right|$). In the first configuration, the ocean depth minimising the critical eccentricity is affected by the drag, while it is not the case for $\sigma_{\rm R}\lesssim 10^{-6}$ ${\rm s^{-1}}$.

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 $\sqrt{g/R_{\rm p}^{2}}$ and $R_{\rm p}^{3}/M_{\rm p}$ 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_{\rm oc}\rightarrow 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_{\rm AR}\approx 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

[TABLE]

and write the tidal torque exerted on the solid body as

[TABLE]

As may be noticed, the sign of the torque in this equation only depends on the rheological exponent $\alpha$, 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_{\rm oc}\rightarrow 0$) in Fig. 4 (top and bottom panels).

The above observations suggest that final rotation states of equilibrium distinct from spin-orbit resonances ($\Omega_{{\rm eq};s}=\left(s/2\right)n_{\star}$) 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 ($\alpha$) in the high-eccentricity regime. The important parameters for these states are the ocean parameters (depth, density, Rayleigh drag frequency) and the orbital period.

6.2 Capture in the 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 – denoted $\mathcal{T}_{{\rm p}}^{\rm tri}$ – adds to the tidal torque – superscripted ’${\rm tide}$’ here to avoid confusion with the total torque – in the right-hand member of the equation of the spin evolution,

[TABLE]

This torque is expressed in the general case as (see Appendix C)

[TABLE]

In the preceding equations, the symbol \leavevmode\nobreak\ \dot{}\leavevmode\nobreak\ is the time derivative, $A$, $B$, and $C$ the principal moments of inertia of the planet (in increasing magnitude), and $\gamma_{2,s}\equiv\vartheta-\left(s/2\right)M$ the angles associated with the eccentricity frequencies $\sigma_{2,s}$, the notation $\vartheta$ and $M$ designating the planet rotation angle and mean anomaly, respectively. The triaxiality, $\left(B-A\right)/C$, is the dimensionless parameter that quantifies the degree of permanent non-rotationally symmetric deformation of the planet. It annihilates 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 by222In the vicinity of the spin-orbit resonance $s/2$, the angle $\gamma_{2,s}$ is small. Therefore, integrating Eq. 76 over an orbital period makes that all eccentricity terms vanish except the degree-$s$ term, which reduces Eq. 75 to the equation of motion of a simple pendulum in the absence of tides (see e.g., Eq. (5) in Goldreich, 1966). However, we note that $\sin\left(\gamma_{2,s}\right)$ vanishes at the resonance ($\vartheta/M=s/2$).

[TABLE]

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

[TABLE]

where we have made use of the relationship $\sigma_{2,s}=2\dot{\gamma_{2,s}}$ and denoted the second order time derivative by \leavevmode\nobreak\ \ddot{}\leavevmode\nobreak\. Assuming that the triaxial torque predominates in the vicinity of the resonance, we neglect the tidal torque ($\left|\mathcal{T}_{\rm p}^{\rm tid}\right|\ll\left|\mathcal{T}_{{\rm p}}^{\rm tri}\right|$). It follows that the maximum absolute value that $\dot{\gamma_{2,s}}$ can reach inside the resonance is given by (Goldreich & Peale, 1966, Eq. (17))

[TABLE]

Following Ribas et al. (2016), $\Delta$ is called the width of the resonance in the following, since the separatrix between the librating (trapped) and circulating states is defined as

[TABLE]

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 $\left|\sigma_{2,s}\right|<2\Delta$.

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

[TABLE]

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)

[TABLE]

where the integrals should be performed over the separatrix between the librating and circulating states, given by Eq. 80. We note that $\mathcal{P}_{{\rm cap;}s}>1$ 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,

[TABLE]

where $P_{\rm rot}^{\Delta}\equiv 2\pi/\left(\Delta+n_{\star}\right)$ designates the rotation period at the bound of the resonance, that is such that $\Omega-n_{\star}=\Delta$.

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): ${\sim}1.4\times 10^{-4}$ for Mercury, as derived by Ribas et al. (2016) from the gravity moments measured by Smith et al. (2012), ${\sim}2\times 10^{-5}$ for the Earth, and ${\sim}6\times 10^{-6}$ for Venus (Yoder, 1995). Similarly, the length of the Solar day at the bound of the resonance ($\dot{\gamma_{2,s}}=\Delta$) is plotted as a function of the orbital period and triaxiality in logarithmic scales (Fig. 7, middle panel). In this two plots, we assume $X_{2}^{-3,2}\left(e\right)=1$, since this Hansen coefficient hardly varies for $e\lesssim 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, $P_{\rm sol}^{\Delta}\equiv P_{\rm rot}^{\Delta}P_{\star}/\left(P_{\star}-P_{\rm rot}^{\Delta}\right)=2\pi/\Delta$ varies from ${\sim}50$ to ${\sim}3000$ Earth days in the range $P_{\star}\sim 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,

[TABLE]

where $\Omega_{{\rm eq}}$ designate 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 6 day orbital period with $\sigma_{\rm R}=10^{-6}$ ${\rm s^{-1}}$ (Fig. 6, top middle panel) as a function of the logarithms of the eccentricity and ocean depth (Fig. 7, right panel).

The very small triaxialities obtained outside of oceanic tidal resonances for $e\lesssim 0.25$ (blue to orange areas) 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.

We characterised in this section 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 the spin-orbit synchronous rotation. We aim to identify the regions of the parameter space where this eccentricity is the smallest.

7 Evolution of the critical eccentricity with the orbital period and ocean depth

The critical eccentricity is the parameter that determines whether asynchronous rotation states of equilibrium may exist outside of spin-orbit resonances ($\Omega_{{\rm eq}}{:}n_{\star}=3{:}2,2{:}1,5{:}2,\ldots$; e.g. Makarov & Efroimsky, 2013; Correia et al., 2014) or not. It provides a qualitative information on the rotation state at which a planet can be found. Thus, we examine in this section how it 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 ($\Omega_{{\rm eq}}\approx n_{\star}$). 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 in this section the critical eccentricity as the minimal eccentricity such that the final rotation period satisfies the condition $P_{\rm rot}/P_{\star}\leq 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 ${\sim}10^{-3}$ (see Fig. 7, left panel).

In Fig. 8, we plot the logarithm of the critical eccentricity calculated using the above definition (top panels) and the associated normalised final rotation period $P_{\rm rot}/P_{\star}$ (bottom panels) as functions of the logarithms of the orbital period (horizontal axis) and ocean depth (vertical axis), for a strong drag ($\sigma_{\rm R}=10^{-5}$ ${\rm s^{-1}}$, left) and a weak drag ($\sigma_{\rm R}=10^{-6}$ ${\rm s^{-1}}$, right). To obtain these maps, we performed the calculations detailed in Sect. 6 for each of the sampled orbital periods with the same parametrisation.

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_{10}}e_{\rm AR}\approx-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

[TABLE]

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

[TABLE]

respectively. We remind here that these scaling laws where derived in the framework of approximations such as the quasi-adiabatic ($\sigma_{\rm R}\ll\left|\sigma\right|$) and non-rotating (no Coriolis effects) approximations, which leads to differences with the results of numerical calculations, and particularly in the strong drag configuration ($\sigma_{\rm R}=10^{-5}$ ${\rm s^{-1}}$).

We now move to the maps representing the final rotation rates associated with the computed critical eccentricity (bottom panels). Yellow areas designate sates 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 emphasized 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. This is the object of the next section.

8 Evolution timescale

We return to the case of the Earth-sized planet orbiting the TRAPPIST-1 star with a 6-days 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 $\sigma_{\rm R}=10^{-5}$ ${\rm s^{-1}}$ and $\sigma_{\rm R}=10^{-6}$ ${\rm s^{-1}}$, respectively. Three orbits of various eccentricity are studied. In the first one, $e=0.063$, which is the value of the critical eccentricity predicted by the theory in the weak drag configuration (see Fig. 6, top middle panel). In the second orbit, $e$ is set to 0.2, that is just below the value of the critical eccentricity obtained in the absence of ocean. In the third orbit, it 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\leq\gamma\leq 0$, where $\gamma\equiv{\log_{10}}\left(H_{\rm oc}\right)$ with $H_{\rm oc}$ given in kilometres. As mentioned in Sect. 6, the range of values $H_{\rm oc}\lesssim 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 bound $\gamma=-2$ is intended to emphasise the continuous transition between the asymptotic regime of dry planets ($H_{\rm oc}\rightarrow 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 $\Omega=\left(5/4\right)n_{\star}$, 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

[TABLE]

which is integrated over time using a variable step size Bulirsch-Stoer method (Press et al., 2007). The normalised moment of inertia $\check{C}=C/\left(M_{\rm p}R_{\rm p}^{2}\right)$ is set to $\check{C}=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_{\rm rot}/P_{\star}$ (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 $\tau_{\rm ev}$ can be estimated by considering the contribution of the solid part in the absence of ocean in the circular case. By assuming that $\sigma_{2,2}\approx n_{\star}/2=\pi/P_{\star}$, we thus obtain

[TABLE]

where the expression of $\Im\left\{k_{2}^{\sigma}\right\}$ as a function of the Andrade parameters is given by Eq. 47. The substitution of parameters by their numerical values leads to $\tau_{\rm ev}\sim 10^{3}$ ${\rm 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 $\sim 5000$ ${\rm yr}$ whatever the depth of their oceanic layer. For $\sigma_{\rm R}=10^{-6}$ ${\rm s^{-1}}$, two of the studied cases avoid the synchronisation owing to the action of oceanic tides and evolve towards asynchronous rotation states located between the $1{:}1$ and $3{:}2$ 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 describe properly 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 explore numerically 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 $\sim\!0.3$ to $\sim\!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 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\ll 0.2$), where planets are driven towards synchronisation except in case of strong resonance of the oceanic tidal response, (ii) the transition regime ($e\approx 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\gtrsim 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 laters implies 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 to quantify the critical eccentricity in the quasi-adiabatic regime, below $\sigma_{\rm R}\sim 10^{-6}$ ${\rm s^{-1}}$. 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 account complicates the analytic treatment of the problem significantly.

Finally, we would like to draw the attention of the reader on 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, one should investigate in future studies 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).

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. 679030/WHIPLASH and 771620/EXOKLEIN), and has been carried out within the framework of the NCCR PlanetS supported by the Swiss National Science Foundation. S. Mathis acknowledges funding by the European Research Council through the ERC grant No. 647383/SPIRE and by the PLATO CNES grant at CEA Saclay. E. Bolmont and P. Auclair-Desrotour 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 $x\in\left[-1,1\right]$ and $l\geq\left|m\right|$ as (Abramowitz & Stegun 1972; Arfken & Weber 2005)

[TABLE]

the $\bar{P}_{l}$ designating the Legendre polynomials, given by

[TABLE]

The normalised associated Legendre functions $P_{l}^{m}$ are characterised by

[TABLE]

where the Kronecker symbol $\delta_{l,k}$ is such that $\delta_{l,k}=1$ if $l=k$ and $\delta_{l,k}=0$ otherwise. Hence they are expressed as

[TABLE]

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)

[TABLE]

where $v$ designate the true anomaly and $r_{\star}$ the star planet distance, these two quantities being themselves expressed as

[TABLE]

Hansen coefficients are thus defined by (e.g. Hughes 1981)

[TABLE]

We note that the $X_{s}^{l,m}$ are real since $g_{l}^{m}$ is an even function of $M$.

By changing the sign of $s$ and using the symmetry property $X_{-s}^{l,-m}=X_{s}^{l,m}$, we obtain

[TABLE]

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 $g_{l}^{m}$, one can efficiently get a numerical estimate of $2^{N}+1$ coefficients for $-K\leq s\leq K$ with $K=2^{N-1}$ by means of a fast Fourier transform (FFT) (Adel Sharaf & Hassan Selim 2010; Correia et al. 2014). The integer $N$ shall be chosen appropriately so that $\max\left\{\left|U_{{\rm T};l}^{m,s}\right|\right\}_{s=\pm K}\ll\max\left\{\left|U_{{\rm T};l}^{m,s}\right|\right\}_{-K\leq s\leq K}$, 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 ($\vartheta$) and true anomaly ($v$) as (e.g. Goldreich & Peale 1966, Eq. (2))

[TABLE]

which can be rewritten as $\mathcal{T}_{{\rm p}}^{\rm tri}=\Im\left\{\tilde{\mathcal{T}}_{{\rm p}}^{\rm tri}\right\}$, where

[TABLE]

By using Eq. 98 and the notation $\gamma_{2,s}=\vartheta-\left(s/2\right)M$, we thus obtain

[TABLE]

and

[TABLE]

Appendix D Nomenclature