© ESO 2019

1 Introduction

The inner planets of the solar system, like the majority of the main satellites, currently present a rotation state that is different from what is believed to have been the initial state (e.g., Goldreich & Soter 1966). Just after their formation, they are supposed to have rotated much faster (e.g., Kokubo & Ida 2007), but subsequently slowed down owing to tidal interactions with the central more massive body (e.g., MacDonald 1964).

The last stage for tidal evolution is synchronous rotation with the orbital mean motion, since it corresponds to the minimum of energy (Hut 1980; Adams & Bloch 2015). This equilibrium, also known as the 1/1 spin-orbit resonance, was first observed for the Moon. However, other possibilities exist, such as the 3/2 spin-orbit resonance of Mercury (Pettengill & Dyce 1965), which is possible because the eccentricity is different from zero (Colombo 1965; Goldreich & Peale 1966; Correia & Laskar 2004; Noyelles et al. 2014). Another interesting example is the chaotic rotation of Hyperion (Wisdom et al. 1984; Harbison et al. 2011), which also results from nonzero eccentricity.

When the eccentricity is nearly zero, asynchronous spin-orbit resonances are still possible, provided that a companion is present in the system. The rotation of Venus is also tidally evolved and not synchronous (Goldstein 1964; Carpenter 1964); it is close to a −5/1 spin-orbit resonance with the Earth, which suggests that the rotation can be trapped at this ratio¹ (Goldreich & Peale 1966, 1967). Indeed, one possibility for asynchronous rotation is a spin-orbit resonance with a conjunction frequency. The minor satellites of the Pluto–Charon system may also have ended up trapped in a resonance of this kind (Correia et al. 2015) or even present chaotic rotation (Showalter & Hamilton 2015). Another possibility, which is very important for planets involved in mean motion resonances, is resonant rotation with the orbital libration frequency (Correia & Robutel 2013; Leleu et al. 2016; Delisle et al. 2017).

Despite the proximity of solar system bodies, the determination of their rotational periods has only been achieved in the second half of the 20th century. Mercury and Venus using radar ranging (Goldstein 1964; Carpenter 1964; Pettengill & Dyce 1965), Hyperion after the Voyager 2 flyby (Smith et al. 1982), and the minor satellites of the Pluto with images from the Hubble Space Telescope (Showalter & Hamilton 2015). We thus do not expect that it will be easy to observe the rotation of recently discovered exoplanets. Nevertheless, many of these planets are close to their host star, and we can assume that their spins have already undergone enough dissipation and evolved into a final equilibrium. Therefore, we can try to anticipate all the possible scenarios for the rotation of these planets. In particular, we aim to determine when asynchronous rotation is possible, which is important for habitability studies.

In this paper we generalize previous studies to any kind of orbital perturbations. In Sect. 2, we write the equation for the rotation of a planet that is perturbed by other bodies, and show that the orbital forcing can be decoupled from the rotation. In Sect. 3, we determine the possible configurations for the rotation of a rigid body by analyzing the orbital forcing. In Sect. 4, we consider that the planet can deform under the action of tides, and present a new method for determining the equilibrium shape and the forced oscillations in resonance. In Sect. 5, using a Maxwell model, we obtain the tidal torque as a function of the orbital forcing and the corresponding capture probabilities in resonance. In Sect. 6, we perform complete numerical simulations including tidal dissipation for the KOI-1599 system. Finally, we discuss our results in Sect. 7.

2 Orbital forcing

We consider a system composed by a central star of mass M, a companion planet of mass m, and additional companions in coplanar orbits. The planet is considered an oblate ellipsoid with mean radius R and gravity field coefficients given by C₂₂ and S₂₂. We furthermore assume that the spin axis of the planet is orthogonal to the orbital plane (zero obliquity). The torque of the star on the planet changes the rotation angle, θ, whose variations are given by (e.g., Correia et al. 2014): $$\ddot{\theta }=-\frac{6GMm{R}^2}{C{r}^3}[C_{22}\sin 2(\theta -f)+S_{22}\cos 2(\theta -f)],$$(1)

where G is the gravitational constant, C is the moment of inertia with respect to the rotation axis, r is the distance and f is the true longitude of the center of mass of the planet with respect to the star. We neglect terms in ${(R/r)}^3$ (quadrupolar approximation) and the torques of the companions on the planet’s spin.

The above expression can be rewritten in a more compact and usable form with complex numbers as $$\ddot{\theta }=\mathcal{B}\Im \left(Zx{e}^{-2i\theta }\right)\text{},$$(2)

with $$Z=C_{22}-i{S}_{22},$$(3) $$x={\left(\frac{a}{r}\right)}^3{e}^{2if},$$(4) $$\mathcal{B}=\frac{6GMm{R}^2}{C{a}^3}\approx 18n^2,$$(5)

where a is a scale constant with units of distance, that is usually made equal to a mean semimajor axis, n is the mean motion, and C∕mR² ≈ 1∕3, which is theusual value for rocky planets (Yoder 1995). This way, we are able to isolate in Eq. (1) all the important quantities that change the rotation, namely, the shape (or deformation) of the planet, Z, the orbital forcing, x, and the rotation angle, e^−2iθ. This factorization allows us to treat each contribution separately.

The orbital and spin evolutions are coupled, such that the total angular momentum is conserved. Therefore, a change in the rotation angle (Eq. (1)) also modifies the position r (see Eqs. (15) and (16) in Correia et al. 2014). However, the corrections in the orbital motion are on the order of ${(R/r)}^2\left|Z\right|\ll 1$ (the orbital angular momentum is much larger than the rotational angular momentum), so their contribution can usually be neglected. For planets in multi-body systems, the important contributions to the orbital forcing thus come from the mutual gravitational perturbations. We let $$X=r{e}^{if}=(r\cos f,r\sin f)$$(6)

be the position vector of the planet, whose spin we want to study. The evolution of X can be obtained by integrating a system of N point-mass bodies (the star, the planet and the remaining companions). The orbital forcing on the spin (Eq. (4)) can then be directly computed from X as $$x=\frac{a^3{X}^2}{|X{|}^5}.$$(7)

Assuming quasi-periodic orbital motion, we can expand the previous expression as a Fourier series $$x={\displaystyle \sum _k}{x}_k{e}^{ik\cdot \nu t},$$(8)

where ν is the vector of fundamental frequencies, and the amplitudes |x_k| decay with ||k||, such that ∑_k|x_k| is convergent. The orbital forcing on the spin, x, can then be characterized by the main terms in the above series.

For the shape of the planet, Z, we can split it in two main contributions: one rigid or permanent, Z₀, due to inhomogeneous mass distribution inside the planet, and another that accounts for the deformation, Z_t, due to external tidal forcing, such that Z = Z₀ + Z_t. In the solar system, we observe that Z₀ is a function of the planet mass. We have |Z₀|~ 0.1 for bodies with diameters smaller than 100 km (e.g., Correia 2018), but larger bodies evolve into spherical shapes (|Z₀ |→ 0); we have |Z₀ | ~ 10⁻⁵ for Mercury (Genova et al. 2013) and |Z₀|~ 10⁻⁶ for the Earth and Venus (Yoder 1995). On the other hand, $\left|Z_t\right|~\mathcal{A}$ (Eq. (16)), that is, it depends on the mass and distance of the star. For Mercury we have |Z_t | ~ 10⁻⁷ ≪|Z₀|, so we can assume Z ≃ Z₀. However, for the Galilean satellites of Jupiter, we have |Z_t|~ 10⁻⁴ ≫|Z₀|, so we can assume Z ≃ Z_t (Correia & Rodríguez 2013).

3 Rigid bodies

As in most previous studies on spin-orbit coupling, let us assume for now that Z = Z₀, either because |Z_t | ≪|Z₀| (case of Mercury) or because Z_t is frozen from a time the planet was more fluid (case of the Moon). Then, Eq. (2) becomes² $$\ddot{\theta }=\mathcal{B}\Im \left(Z_0{\displaystyle \sum _k}{x}_k{e}^{i(k\cdot \nu t-2\theta )}\right).$$(9)

Each term in this expression individually behaves like a pendulum, with islands of rotational libration centered at $$\langle \dot{\theta }\rangle ={\nu }_k=\frac{k\cdot \nu }{2},$$(10)

and libration width $${\beta }_k=\sqrt{2\mathcal{B}\left|Z_0{x}_k\right|}\approx 6n\sqrt{\left|Z_0\right|}\sqrt{\left|x_k\right|}.$$(11)

In general, the rotation can be trapped in any of these islands, provided they are far apart. However, when the libration widths of some individual resonant islands overlap, the rotation can be chaotic (Chirikov 1979; Morbidelli 2002), i.e., the rotation exhibits random variations in short periods of time. Overlap between two resonances with ${\nu }_{k_1}<{\nu }_{k_2}$ then occurs for ${\nu }_{k_1}+{\beta }_{k_1}\approx {\nu }_{k_2}-{\beta }_{k_2}$, that is, when $$\overline{\beta }={\beta }_{k_1}+{\beta }_{k_2}\approx {\nu }_{k_2}-{\nu }_{k_1}=\Delta \nu .$$(12)

This simple condition allows us to differentiate easily between the different dynamical regimes. Indeed, when Z₀ is known for a given kind of body, we can infer its dynamics simply by analyzing the orbital forcing (Eq. (8)). For $\overline{\beta }\ll \Delta \nu$ the rotation can be trapped in any isolated island with $\dot{\theta }={\nu }_k$, for $\overline{\beta }~\Delta \nu$ the rotation is chaotic, while for $\overline{\beta }\gg \Delta \nu$ the rotation can be stable again because the ν_k are engulfedin a single libration island with only a small chaotic zone at the separatrix.

3.1 Single planet case

In the classical case of a single planet in an elliptical orbit, the orbital forcing (Eq. (4)) is given by $$x={\displaystyle \sum _k}{\chi }_k^{-3,2}(e)e^{2i\varpi }{e}^{iknt},$$(13)

where ${\chi }_k^{l,m}(e)$ are the Hansen coefficients (Hughes 1981) and ϖ is the longitude of the pericenter. Thus, resonant islands occur for ν_k = kn∕2, with amplitude $\left|x_k\right|={\chi }_k^{-3,2}(e)$. The two largest terms correspond to k = 2 and k = 3, for which ${\chi }_2^{-3,2}(e)\approx 1$ and ${\chi }_3^{-3,2}(e)\approx 7e/2$, respectively (we neglect terms in e²). Then $$\overline{\beta }/\Delta \nu \approx 12\sqrt{\left|Z_0\right|}(1+\sqrt{7e/2}).$$(14)

For Mercury (|Z₀|≈ 10⁻⁵, e = 0.206), we get $\overline{\beta }/\Delta \nu \approx$ 0.07 ≪ 1, so we can expect stable rotation. Indeed, the planet is observed in a stable spin-orbit resonance with k = 3 (Colombo 1965; Goldreich & Peale 1966). On the other hand, for Saturn’s moon Hyperion (|Z₀ |≈ 0.03, e = 0.123), we get $\overline{\beta }/\Delta \nu \approx 3.5~1$. According to condition (12), resonance overlap occurs for Hyperion, which shows chaotic rotation (Harbison et al. 2011).

The single planet case can be treated analytically. However, we can also study it numerically, in order to understand more complex situations. We can integrate the orbital motion of the two-body problem with an eccentricity e = 0.165 (mid-value of the orbits of Mercury and Hyperion), which provide us a time series for X. We then compute a fast Fourier transform of the quantity x (Eq. (7)). The frequencies ν_k∕n = k∕2 immediately show up as well as the corresponding amplitudes |x_k|. Using expression (11) we obtain the rotational libration width corresponding to the resonance islands $\dot{\theta }={\nu }_k$ for a given value of |Z₀|. In Fig. 1 we horizontally draw the width of the libration islands centered at each ν_k for the |Z₀| values of Mercury and Hyperion. This plot clearly shows why Mercury can be trapped in isolated spin-orbit resonances, while Hyperion shows chaotic motion.

Fig. 1 Fast Fourier transform of the orbital forcing x (Eq. (7)) for the two-body problem with eccentricity e = 0.165. The horizontal lines give the amplitudes βk for |Z0 | = 1.5 × 10−5, which corresponds to Mercury’s value (Genova et al. 2013), and |Z0| = 2.9 × 10−2, which corresponds to Hyperion’s value (Thomas et al. 2007).

3.2 Multi-planet systems

In general, planets are not alone in their orbits. For widely spaced orbits, as in the solar system, planet–planet mutual interactions are relatively weak, and only modify the orbital parameters on secular timescales. As a result, the single planet case from the previous section is a good approximation. However, a large number of exoplanets has been observed in compact configurations and/or showing resonant interactions. For these systems, the mutual perturbations are strong and need to be taken into account in the spin-orbit coupling analysis.

One possibility is to expand the disturbing function in its main perturbing terms (see Chap. 6, Murray & Dermott 1999), and use it to obtain the dominating contributions to x (Eq. (8)). Closely packed systems usually have small eccentricities, so a common approach is to truncate the series for some power in the eccentricities (e.g., Goldreich & Peale 1966, 1967; Correia & Robutel 2013; Correia et al. 2015; Leleu et al. 2016; Delisle et al. 2017). Although modern series manipulators allow us to perform this calculation rapidly, the number of terms to retain is large (even for a two-planet system), and their amplitudes depend on the masses and semimajor axes. Moreover, for resonant systems the eccentricities can be high and the libration amplitude is also important. Therefore, the dominating terms vary from system to system and it is difficult to predict the terms to be retained.

A powerful alternative to analytical studies is to obtain the series expansion for x numerically. This method is exact (without approximations), fast, and valid for an unlimited number of planets in the system. We simply need to run a simulation of the full system and then perform a frequency analysis of the orbital motion. The only requirement is that we integrate the system for a time long enough to extract the frequencies whose amplitudes may be dominating. However, very small frequencies (corresponding to very long periods) have no perceptible effect on the spin, since they are engulfed by the synchronous resonance (case where $\overline{\beta }\gg \Delta \nu$). Thus, in general a short integration of the system (10⁴ times the longest orbital period) is enough to characterize the spin dynamics.

As an example, we apply the numerical study to the KOI-1599 system, whose planets were detected using the transit method (Steffen & Hwang 2015) and validated with transit time variations (Panichi et al. 2019). This system is a good candidate for a rich spin dynamics because it shows some transit time variations (Delisle et al. 2017). The best stable fits yield a two super-Earth system in nearly circular orbits and involved in a 3/2 mean motion resonance (Table 1). We integrated the KOI-1599 system for 1000 years, and then performed a frequency analysis of the forcing function x (Laskar 1993). In Table 2, we show the results for the outer orbit because Δν∕n is larger, i.e., the individual resonant islands are more distant. We observe that the classic perturbing terms at half-integers of the mean motion (Sect. 3.1) are still present, but with very low amplitude, since the eccentricity is always very small (e < 0.025). The three dominating terms correspond to harmonics of the mean motion, n, and the orbital libration frequency, υ, with similar amplitudes.

We need to estimate the rotational libration width β_k for each frequency in Table 2 to determine the behavior of the spin. Since we know the orbital forcing amplitudes, the only uncertainty is related to the permanent deformation of the planet, Z₀ (Eq. (11)). We can adopt the Earth’s value, |Z₀| = 1.5 × 10⁻⁶ because KOI-1599 planets are likely rocky, as their densities are similar to that of the Earth. However, their surface gravity is a bit stronger, so maybe the inhomogeneities in the mass distribution are smaller, such that |Z₀ |~ 10⁻⁷. These two values can be seen as maximal and minimal estimations for |Z₀|, respectively. In Fig. 2 we horizontally draw the libration width of the islands centered at each ν_k for both |Z₀ | values over the Fourier spectra of the orbital forcing x. We observe that, for |Z₀|~ 10⁻⁷, each island is isolated, so the rotation can be trapped in any of them, which includes the synchronous resonance, along with the super- and sub-synchronous resonances. On the other hand, for |Z₀ |~ 10⁻⁶, the libration widths overlap, so we can expect chaotic rotation.

In order to study in more detail the spin dynamics, we have to run direct numerical simulations of the rotation angle of the outer planet. We then integrated Eq. (1) together with the orbital solution of the KOI-1599 system given by a classical N−body integrator. In Fig. 3 we show a frequency map analysis of the rotation angle (Laskar 1993; Leleu et al. 2016) for the two |Z₀| values in the neighborhood of the synchronous resonance. We conclude that our predictions for the spin of the outer planet based on the Chirikov (1979) general criterion are correct: for |Z₀ |~ 10⁻⁷ the planet can present stable asynchronous rotation, while for |Z₀|~ 10⁻⁶ it likely presents chaotic rotation.

Fig. 2 Fast Fourier transform of the orbital forcing x (Eq. (7)) for the outer planet in the KOI-1599 system. The horizontal lines give the amplitudes βk for |Z0 |= 1.5 × 10−6 and |Z0 | = 1.0 × 10−7.

Fig. 3 Frequency map analysis of KOI-1599 outer planet rotation angle assuming a rigid body with |Z0 | = 1.5 × 10−6 (top) or |Z0| = 1.0 × 10−7 (bottom). Each dot corresponds to a grid with 181 × 201 different initial conditions of γ(0) and $\dot{\gamma }(0)$, respectively, with γ(t) = θ(t) − λ(0). The color index gives the derivative $\partial \eta /\partial \dot{\gamma }(0)$, where η is the main frequency of γ(t). The blue dots correspond to resonant motion, green dots to nonresonant regular motion, and red dots to chaotic motion.

4 Tidal deformation and forced oscillations

Planets are not completely rigid bodies. They deform under tidal interactions with the parent star because the mass distribution inside the planet adjusts to the tidal potential (Darwin 1879, 1880). Viscoelastic rheologies have been shown to reproduce the main features of tidal effects. The planet can respond as an elastic solid or as a viscous fluid, depending on the frequency of the perturbation (for a review see Henning et al. 2009).

One of the simplest viscoelastic models is to consider that the planet behaves like a Maxwell material, which is represented by a purely viscous damper and a purely elastic spring connected in series (e.g., Turcotte & Schubert 2002). More complex models, such as the Andrade model (Andrade 1910), can be thought of as the Maxwell model equipped with extra terms describing hereditary reaction of strain to stress (Efroimsky 2012). Therefore, the main features of tidal response are similar to the Maxwell model.

For simplicity, we adopt the Maxwell viscoelastic rheology. The viscous response of the planet to the tidal excitation is modeled by the parameter τ, which corresponds to the relaxation time of the planet, while the elastic response in modeled by the parameter τ_e. The deformation Z follows (Correia et al. 2014) $$Z+\tau \dot{Z}=Z_0+Z_{\text{e}}+{\tau }_{\text{e}}{\dot{Z}}_{\text{e}},$$(15)

with $$Z_{\text{e}}=\mathcal{A}\overline{x}{e}^{2i\theta },\text{\hspace{1em}where\hspace{1em}}\mathcal{A}=\frac{k_{\text{f}}}{4}\frac{M}{m}{\left(\frac{R}{a}\right)}^3,$$(16)

k_f is the fluid second Love number for potential, and $\overline{x}$ is the complex conjugate of x. The function Z_e is the equilibrium shape, i.e., the maximal deformation that the planet can attain owing to the tidal perturbation.

We denote the average rotation rate of the planet by $\omega =\langle \dot{\theta }\rangle$. For simplicity, we assume that θ has a small forced amplitude of libration due to orbital forcing and a zero free amplitude of libration, such that $$\begin{array}{lll}\hfill & \hfill & \theta =\omega t+{\displaystyle \sum _k}{\theta }_k{e}^{ik\cdot \nu t},\hfill \\ \hfill & \hfill & e^{2i\theta }=e^{2i\omega t}{\displaystyle \sum _{j\ge 0}}\frac{{\left(i{\displaystyle \sum _k}2{\theta }_k{e}^{ik\cdot \nu t}\right)}^j}{j!}={\displaystyle \sum _l}{y}_l{e}^{i(2\omega +l\cdot \nu )t},\hfill \end{array}$$

where all θ_k are small, and the series is absolutelyconvergent (∑_k|θ_k| converges). Each coefficient y_l can be expressed as a power series of the coefficients θ_k, by identifying terms on both sides of Eq. (18).

From expressions (8), (16), and (18), we obtain $$Z_{\text{e}}=\mathcal{A}{\displaystyle \sum _k}\left({\displaystyle \sum _l}{\overline{x}}_{k+l}{y}_l\right)e^{i(2\omega -k\cdot \nu )t},$$(19)

and replacing this in expression (15), we get $$Z=Z_0+Z_t+C{e}^{-t/\tau },$$(20)

with $$Z_t=\mathcal{A}{\sum }_k\left({\displaystyle \sum _l}{\overline{x}}_{k+l}{y}_l\right)\frac{1+i(2\omega -k\cdot \nu ){\tau }_{\text{e}}}{1+i(2\omega -k\cdot \nu )\tau }{e}^{i(2\omega -k\cdot \nu )t}.$$(21)

In the following, we neglect the term Ce^−t∕τ because it vanishes after some time, and we only consider the solution Z = Z₀ + Z_t.

4.1 No resonant rotation

We assume that ω≠ν_k = k ⋅ν∕2 for all k (no resonant rotation). Then, over one rotational period we have (Eq. (21)) $$\langle Z_t\rangle =0,$$(22)

that is, there is no permanent tidal deformation.

4.2 Resonant rotation without forced oscillations

The situation changes if the rotation is locked in some spin-orbit resonance, that is, we assume that ω = ν_r = r ⋅ν∕2 (Eq. (10)). We can then rewrite expression (20) as $$Z={\displaystyle \sum _k}{z}_k{e}^{i(r-k)\cdot \nu t},$$(23)

with $$z_k=Z_0{\delta }_{r,k}+\mathcal{A}{b}_{r-k}{\displaystyle \sum _l}{\overline{x}}_{k+l}{y}_l,$$(24)

and $$b_k=\frac{1+i(k\cdot \nu ){\tau }_{\text{e}}}{1+i(k\cdot \nu )\tau }.$$(25)

The b_k function depends on the tidal model, in this case, the Maxwell model. Different tidal models would imply different expressions for b_k, but for all models we get b₀ = 1.

In the case of very weak or no forced oscillations, expressions (17) and (18) simplify as $$\begin{array}{lll}\hfill & \hfill & \theta ={\theta }_0+\frac{r\cdot \nu }{2}t,\hfill \\ \hfill & \hfill & e^{2i\theta }=e^{2i{\theta }_0}{e}^{ir\cdot \nu t},\hfill \end{array}$$

that is, y₀ = e^2iθ₀, and y_k≠0 = 0. Replacing in expression (24), we get for the deformation of the planet $$z_k=Z_0{\delta }_{r,k}+\mathcal{A}{b}_{r-k}{e}^{2i{\theta }_0}{\overline{x}}_k.$$(28)

Finally, to get the permanent tidal deformation, we can average over one rotational period $$\langle Z_t\rangle =\langle Z\rangle -Z_0=\mathcal{A}{e}^{2i{\theta }_0}{\overline{x}}_r,$$(29)

or $$|\langle Z_t\rangle |=\mathcal{A}\left|x_r\right|.$$(30)

For instance, in the classical case of a single planet in an elliptical orbit (Sect. 3.1), we have (Eq. (13)) $${\overline{x}}_k={\chi }_k^{-3,2}(e)e^{-2i\varpi },$$(31)

thus, for a planet trapped in the ν_k∕n = k∕2 resonance $$|\langle Z_t\rangle |=\mathcal{A}{\chi }_k^{-3,2}(e).$$(32)

For the outer planet in the KOI-1599 system (Sect. 3.2), we have for the synchronous resonance (Table 2) $$\left|x_{[2,0]}\right|\approx 0.93\Rightarrow |\langle Z_t\rangle |\approx 2.5\times {10}^{-6},$$(33)

while for the super- and sub-synchronous resonances $$\left|x_{[2,\pm 1]}\right|\approx 0.25\Rightarrow |\langle Z_t\rangle |\approx 7.0\times {10}^{-7}.$$(34)

We note that ⟨Z_t⟩ (Eq. (29)) does not depend on the tidal model function b_k (Eq. (25)). Indeed, ⟨Z_t⟩ is independent of the tidal model, as it corresponds to the maximal deformation that the planet can achieve while trapped in a given resonance.

4.3 Resonant rotation with forced oscillations

In case of resonant motion (ω = ν_r), replacing the deformation of the planet, Z (Eq. (23)), in expression (2), we get $$\ddot{\theta }=\mathcal{B}\Im \left({\displaystyle \sum _k}\left({\displaystyle \sum _{l,m}}{x}_{k+l+m}{\overline{y}}_l{z}_m\right)e^{ik\cdot \nu t}\right),$$(35)

and from expression (17) $$\ddot{\theta }=-{\displaystyle \sum _k}{(k\cdot \nu )}^2{\theta }_k{e}^{ik\cdot \nu t}.$$(36)

Then, by identifying terms in both series, we find $${\theta }_k=\frac{1}{2i}\left({\alpha }_k-{\overline{\alpha }}_{-k}\right),$$(37)

where $${\alpha }_k=\frac{\mathcal{B}}{{(k\cdot \nu )}^2}{\displaystyle \sum _{l,m}}{\overline{x}}_{l+m-k}{y}_l{\overline{z}}_m.$$(38)

We also note that ${\theta }_{-k}={\overline{\theta }}_k$. Therefore, with ϕ_k = arg α_k, we can rewrite expression (17) as $$\theta ={\nu }_{r}t+{\theta }_0+{\displaystyle \sum _{k\ne 0}}\left|{\alpha }_k\right|\sin (k\cdot \nu t+{\varphi }_k).$$(39)

We recall that the coefficients y_l are expressed as power series of the coefficients θ_k (Eq. (18)), and that the coefficients z_k are linear combinations of the coefficients y_l (Eq. (24)). Equations (18), (24), and (37) thus define a closed system of equations, which can be solved for the amplitudes of the forced oscillations θ_k and for the deformation of the planet z_k.

We can estimate the order of magnitude of the expected forced oscillations. Since |x_k | ≤ 1, |y_k | ≤ 1 and $\left|z_k\right|\le \mathcal{Z}=\max (|Z_0|,\mathcal{A})$, we deduce fromexpression (37) $$\left|{\theta }_k\right|\le \frac{\mathcal{B}\mathcal{Z}}{{(k\cdot \nu )}^2}\approx 18\mathcal{Z}{\left(\frac{n}{2{\nu }_k}\right)}^2.$$(40)

For an Earth-size planet and a perturbing frequency on the order of the orbital frequency (ν_k ~ n), we get |θ_k | ≲ 10⁻⁵, that is, the amplitude of the forced oscillations remains very small. We then linearize expression (18) and find $$y_0=e^{i2{\theta }_0},\text{\hspace{1em}and\hspace{1em}}{y}_{l\ne 0}\approx 2i{\theta }_l{e}^{i2{\theta }_0}=\left({\alpha }_l-{\overline{\alpha }}_{-l}\right)e^{i2{\theta }_0}.$$(41)

The coefficients z_k thus become linear in θ_k, so we only need to solve a system of linear equations to get the amplitudes of forced libration, θ_k, and the corresponding corrections to the deformation of the planet, z_k. For planets involved in mean-motion resonances, we may have ν_k∕n ~ 0.01 (Delisle et al. 2017), so the amplitude of forced oscillations is expected to be stronger (|θ_k |≲ 0.1). In this case, higher order terms should be included and we would need to solve a system of polynomial equations (see Delisle et al. 2017, Appendix D).

For the sake of simplicity and readability, we restrict our study to the first order case in the following, as it already provides a good approximation in most cases. Then, from expressions (38) and (41) we have $${\alpha }_k\approx \frac{18n^2{e}^{i2{\theta }_0}}{{(k\cdot \nu )}^2}{\displaystyle \sum _m}{\overline{x}}_{m-k}{\overline{z}}_m,$$(42)

where ${\overline{z}}_m$ is given by expression (24). In order to be captured in a given resonance (ω = ν_r), we have that x_r ≠ 0 (Sect. 3). Since we always have k≠0 (Eq. (39)), it is needed that at least another x_m≠r≠0 to observe forced librations. Indeed, if there is only one orbital forcing term, as in the case of a single planet in a circular orbit, all α_k = 0, and hence there are no forced librations.

4.3.1 Single planet case

We first consider a single planet in an elliptical orbit with small eccentricity. The two dominating terms of the orbital forcing are x₂ ≈ 1 and x₃ ≈ 7e∕2 (Sect. 3.1), so we neglect the remaining terms. Applying to the case of Mercury, we have r = 3, ν = n, and $\left|Z_0\right|\gg \mathcal{A}$. Thus, expression (42) becomes $${\alpha }_k\approx \frac{18e^{i2{\theta }_0}}{k^2}{\overline{Z}}_0{\overline{x}}_{3-k}.$$(43)

The only α_k≠0 is for x_3−k = x₂, that is, for k = 1. Then, the rotation angle becomes (Eq. (39)) $$\theta \approx \frac{3}{2}nt+{\theta }_0+18\left|Z_0{x}_2\right|\sin (nt+{\varphi }_1).$$(44)

This simple approximation gives a forced libration period equal to the orbital period, and a forced libration amplitude |α₁ |≈ 18 |Z₀| = 56″, in perfect agreement with a full numerical simulation for Mercury (see Fig. 2 in Peale et al. 2007).

4.3.2 Multi-planet systems

We considerthe outer planet in the KOI-1599 system, for which the dominating terms of the orbital forcing are x_[2,0] = 0.93, x_[2,1] = 0.26, and x_[2,−1] = 0.25 (Table 2), and we neglect the remaining terms. In this case, we have $\left|Z_0\right|\ll \mathcal{A}$, so replacing expression (24) in (42) gives $${\alpha }_k\approx \frac{18\mathcal{A}{n}^2{e}^{i2{\theta }_0}}{{(k\cdot \nu )}^2}{\displaystyle \sum _m}{b}_{m-r}{\overline{x}}_{m-k}{\displaystyle \sum _l}{x}_{m+l}{\overline{y}}_l,$$(45)

which depends on the tidal model through b_m−r (Eq. (25)). For simplicity, we further assume that τ_e = 0 and τν_k ≫ 1 (high frequency regime, see Sect. 5), such that b_m−r ≈ δ_r,m. Then $${\alpha }_k\approx \frac{18\mathcal{A}{n}^2{e}^{i2{\theta }_0}}{{(k\cdot \nu )}^2}{\overline{x}}_{r-k}{\displaystyle \sum _l}{x}_{r+l}{\overline{y}}_l.$$(46)

Let us consider, for example, that the planet is captured in the r = [2, 1] spin-orbit resonance, with frequency ν_[2,1] = n + υ∕2. The only two α_k≠0 are for ${\overline{x}}_{[2,1]-k}={\overline{x}}_{[2,0]}$ and ${\overline{x}}_{[2,-1]}$, that is, for k = [0, 1] and [0, 2], respectively, with associated frequencies υ and 2υ, where υ is the orbital libration frequency $${\alpha }_{\upsilon }\approx 18\mathcal{A}{e}^{i2{\theta }_0}{\left(\frac{n}{\upsilon }\right)}^2{\overline{x}}_{[2,0]}{\displaystyle \sum _l}{x}_{[2,1]+l}{\overline{y}}_l,$$(47) $${\alpha }_{2\upsilon }\approx 18\mathcal{A}{e}^{i2{\theta }_0}{\left(\frac{n}{2\upsilon }\right)}^2{\overline{x}}_{[2,-1]}{\displaystyle \sum _l}{x}_{[2,1]+l}{\overline{y}}_l.$$(48)

Since |α_2υ∕α_υ|≈ 0.067, we can also neglect the contribution of α_2υ. Thus, the rotation angle (Eq. (39)) simply becomes $$\theta \approx (n+\frac{1}{2}\upsilon )t+{\theta }_0+\left|{\alpha }_{\upsilon }\right|\sin (\upsilon t+{\varphi }_{\upsilon }),$$(49)

and the rotation rate $$\dot{\theta }\approx n+\upsilon \left[\frac{1}{2}+\left|{\alpha }_{\upsilon }\right|\cos (\upsilon t+{\varphi }_{\upsilon })\right].$$(50)

From expression (41) we additionally have $$\begin{array}{lll}{\displaystyle \sum _l}{x}_{[2,1]+l}{\overline{y}}_l\hfill & =\hfill & x_{[2,1]}{\overline{y}}_{[0,0]}+x_{[2,0]}{\overline{y}}_{[0,-1]}+x_{[2,-1]}{\overline{y}}_{[0,-2]}\hfill \\ \hfill & \approx \hfill & e^{-2i{\theta }_0}\left(x_{[2,1]}-x_{[2,0]}{\alpha }_{\upsilon }\right).\hfill \end{array}$$(51)

Inserting in expression (47) and solving for α_υ we finally get $${\alpha }_{\upsilon }\approx \frac{18\mathcal{A}{\left(n/\upsilon \right)}^2{\overline{x}}_{[2,0]}{x}_{[2,1]}}{1+18\mathcal{A}{\left(n/\upsilon \right)}^2|x_{[2,0]}{|}^2}.$$(52)

With $\mathcal{A}=2.7\times {10}^{-6}$ (Eq. (16)), and (n∕υ) ≈ 150 (Table 2), we estimate |α_υ|≈ 0.14 rad. From expression (50) we also estimate $\Delta \dot{\theta }/n=(\upsilon /n)\left|{\alpha }_{\upsilon }\right|\approx 9\times {10}^{-4}$. In Fig. 4 we show the forced librations of the rotation rate for the r = [2, 1] spin-orbit resonance, which is obtained performing a numerical simulation of the KOI-1599 system with τn = 10⁴ (that is, τυ ~ 10²). We observe that the libration frequency of the spin is indeed υ, and that the forced libration amplitude also matches the estimated value.

We can now estimate the correction that the forced libration introduces in the tidal deformation (Eq. (24)) as follows: $$\begin{array}{lll}\langle Z_t\rangle \hfill & \approx \hfill & \mathcal{A}\left({\overline{x}}_{[2,1]}{y}_{[0,0]}+{\overline{x}}_{[2,0]}{y}_{[0,-1]}\right)\hfill \\ \hfill & =\hfill & \mathcal{A}{e}^{2i{\theta }_0}\left({\overline{x}}_{[2,1]}-{\overline{x}}_{[2,0]}{\overline{\alpha }}_{\upsilon }\right)\hfill \\ \hfill & =\hfill & \mathcal{A}{e}^{2i{\theta }_0}{\overline{x}}_{[2,1]}\left(1-\frac{18\mathcal{A}{\left(n/\upsilon \right)}^2|x_{[2,0]}{|}^2}{1+18\mathcal{A}{\left(n/\upsilon \right)}^2|x_{[2,0]}{|}^2}\right).\hfill \end{array}$$(53)

We thus obtain for the outer planet in the KOI-1599 system, |⟨Z_t ⟩|≈ 3.6 × 10⁻⁷, which is about half of the value without considering the forced librations (Eq. (34)); this value corresponds to what is observed in the numerical simulations (Fig. 4). We conclude that the forced libration contributes to decreasing the tidal deformation (see also Delisle et al. 2017, Appendix D).

Fig. 4 Forced libration of the rotation rate for the super-synchronous resonance (νr = n+ υ∕2) of the outer planet in the KOI-1599 system. We show the evolution of $(\dot{\theta }-{\nu }_r)/n$ adopting τe = 0 and τn = 104.

5 Tidal dissipation and capture probabilities

The tidal deformation of the planet is not instantaneous; there is a delay between the perturbation and maximal deformation that leads to a nonzero net torque on the spin of the planet. As a consequence, tidal effects slowly modify the average rotation rate of the planet and may drive it into the different configurations described in Sect. 3.

5.1 Tidal dissipation

Replacing the tidal deformation of the planet Z_t (Eq. (21)) in the general expression of the torque (Eq. (2)), and averaging over the fundamental frequencies of the system we get $$\begin{array}{lll}\mathcal{T}\hfill & =\hfill & \langle \mathcal{B}\Im \left(Z_{t}x{e}^{-2i\theta }\right)\rangle \hfill \\ \hfill & =\hfill & \mathcal{A}\mathcal{B}\Im \left({\displaystyle \sum _{k,l,m}}{x}_{k+l}{\overline{x}}_{k+m}{y}_m{\overline{y}}_l\frac{1+i(2\omega -k\cdot \nu ){\tau }_{\text{e}}}{1+i(2\omega -k\cdot \nu )\tau }\right).\hfill \end{array}$$(54)

In addition, assuming very small or no forced oscillations, that is, θ ≈ ⟨θ⟩ = θ₀ + ω t (Eq. (17)), we have y₀ = e^2iθ₀ and y_k≠0 = 0, so the average net torque can be simplified as $$\mathcal{T}(\omega )=-\mathcal{K}{\displaystyle \sum _k}|x_k{|}^2\frac{(2\omega -k\cdot \nu )\tau }{1+{(2\omega -k\cdot \nu )}^2{\tau }^2},$$(55)

where $$\mathcal{K}=\mathcal{A}\mathcal{B}\left(1-\frac{{\tau }_{\text{e}}}{\tau }\right)=k_{\text{f}}\frac{3G{M}^2{R}^5}{2C{a}^6}\left(1-\frac{{\tau }_{\text{e}}}{\tau }\right).$$(56)

The initial rotation rate of the planets in unknown. A small number of large impacts at the end of the formation process can significantly change the spin rate or direction (Dones & Tremaine 1993; Kokubo & Ida 2007). However, since planets contract from an initial larger protoplanet, in general it is expected that they rotate fast. Indeed, based on the solar system planets, MacDonald (1964) proposed an empirical relation ω₀ ∝ m^4∕5∕R², which leads to about 12 h for Earth-size planets. Therefore, in general we expect ω₀ ≫k ⋅ν, and since τ_e < τ, the initial rotation rate can only decrease (Eq. (55)).

As the planet arrives in a regime where ω ~k ⋅ν the situationchanges. Depending on τ and the amplitudes |x_k|, the torque may become positive, and hence create stable equilibria for the rotation rate. In Fig. 5 we draw the net torque as a function of τ for two different orbital forcing functions. The left-hand column represents the case of Mercury, for which $\left|x_k\right|={\chi }_k^{-3,2}(e)$ with e = 0.206, while the right-hand column is for the outer planet in the KOI-1599 system, for which |x_k | are given in Table 2.

We observe that, for small τ values, only one equilibrium point near the synchronous rotation is possible. However, for large τ values, several equilibria are possible, and they coincide exactly with the resonant islands identified in Sect. 3. This result was already known for the classical case of a single planet in an eccentric orbit (see Correia et al. 2014), but it is remarkable that it prevails for any kind of orbital forcing. This property of the tidal torque enhances the chances of observing the rotation of close-in planets trapped in some asynchronous resonance.

5.2 Capture probabilities

The total torque on the rotation (Eq. (2)) can be written as $$\ddot{\theta }=\mathcal{B}\Im \left(Z_{0}x{e}^{-2i\theta }\right)+\mathcal{B}\Im \left(Z_{t}x{e}^{-2i\theta }\right),$$(57)

since Z = Z₀ + Z_t. The average of the second term is the tidal torque $\mathcal{T}(\omega )$ (Eq. (55)). The average of the first term is different from zero when the rotation rate is close to some resonant island, ω ≈ν_r, with ν_r = r ⋅ν∕2 (Eq. (9)): $$\langle \mathcal{B}\Im \left(Z_{0}x{e}^{-2i\theta }\right)\rangle =\mathcal{B}\Im \left(Z_0{x}_r{e}^{-2i(\gamma -\varphi )}\right),$$(58)

where γ = θ − ν_r t + ϕ, and $\varphi =\frac{1}{2}\arg (Z_0{x}_r)$ is a constant phase. Thus, the average of the total torque (Eq. (57)) becomes $$\ddot{\gamma }=-\frac{{\beta }_r^2}{2}\sin 2\gamma +\mathcal{T}({\nu }_r+\dot{\gamma }),$$(59)

where ${\beta }_r^2=2\mathcal{B}\left|Z_0{x}_r\right|$ (Eq. (11)). If Z₀ = 0, capture could only be possible at the equilibrium of the tidal torque, that is, when $\mathcal{T}=0$. However, in a more general situation, the tidal torque can be counterbalanced by the conservative torque resulting from β_r. Therefore, even when $\mathcal{T}\ne 0$, each time the planet crosses a resonant island, there is some chance of being captured.

Goldreich & Peale (1966) computed a first estimation of the capture probability P_cap, and subsequently more detailed studies proved their expression to be essentially correct (for a review, see Henrard 1993). We follow the same approach described in Sect. 5.3 of Correia & Laskar (2010). Expression (59) can be rewritten as $$\frac{\text{d}h}{\text{d}t}=2\mathcal{T}\frac{\text{d}\gamma }{\text{d}t},$$(60)

where $$h={\dot{\gamma }}^2-\frac{{\beta }_r^2}{2}\cos 2\gamma$$(61)

is an integral of motion in the conservative problem (Sect. 3), which can be related to the total energy. The separatrix equation between dynamical regimes is $h={\beta }_r^2/2$, where $h>{\beta }_r^2/2$ gives the trajectories in the circulation zone (outside the resonance) and $h<{\beta }_r^2/2$ the trajectories in the libration zone (captured in resonance). The “energy” variation of the planet after a cycle around the resonance is then given by the function $$\Delta h({\gamma }_1)=2{\displaystyle {\int }_{{\gamma }_1}^{{\gamma }_2}\mathcal{T}}\text{\hspace{0.17em}d}\gamma +2{\displaystyle {\int }_{{\gamma }_2}^{{\gamma }_3}\mathcal{T}}\text{\hspace{0.17em}d}\gamma ,$$

where γ₁ is the γ value when the planet crosses the separatrix between the circulation and libration zones. The values γ₂ and γ₃ are the first two following γ values corresponding to $\dot{\gamma }=0$. In order to be captured, after a cycle inside the resonance, the planet must remain within the libration zone, that is, the “energy” variation after a cycle must be negative Δh(γ₁) < 0. Assuming that γ₁ is uniformly distributed between − π∕2 and π∕2, capture inside the resonance occurs for “energy” variations within $$\left\{E:\Delta h(-\pi /2)<\Delta h({\gamma }_1)<0\right\}\text{},$$(62)

from a total of possibilities $$\left\{E_T:\Delta h(-\pi /2)<\Delta h({\gamma }_1)<\Delta h(\pi /2)\right\}\text{},$$(63)

that is, $$P_{\text{cap}}=\frac{{\displaystyle {\int }_E\text{d}}{\gamma }_1}{{\displaystyle {\int }_{E_T}\text{d}}{\gamma }_1}=\frac{\Delta h(-\frac{\pi }{2})}{\Delta h(-\frac{\pi }{2})-\Delta h(\frac{\pi }{2})}.$$(64)

Usually the tidal torque is very weak, that is, $|\mathcal{T}(\dot{\gamma })|\ll {\beta }_r^2$, which gives γ₃ ≃−π∕2, γ₂ ≃ π∕2, and $h\approx {\beta }_r^2/2$. Then, from expression (61) we have $$\dot{\gamma }=\text{sign}(\dot{\gamma }){\beta }_r\cos \gamma ,$$(65)

and replacing it in expression (64) we finally get $$P_{\text{cap}}=1-\frac{T_{-}}{T_{+}},\text{\hspace{1em}with\hspace{1em}}{T}_{\pm }={\displaystyle {\int }_{-\frac{\pi }{2}}^{\frac{\pi }{2}}\mathcal{T}}\left({\nu }_r\pm {\beta }_r\cos \gamma \right)\text{d}\gamma .$$(66)

This expression is valid for any kind of tidal torque. For the Maxwell rheology that we adopt in this paper, we use expression (55) for $\mathcal{T}$, which leads to $$T_{\pm }=-\mathcal{K}{\displaystyle \sum _k}|x_k{|}^2\left[\frac{1}{{\Delta }_k^{+}}\ln \left(\frac{1+{\delta }_k^{+}}{1-{\delta }_k^{+}}\right)-\frac{1}{{\Delta }_k^{-}}\ln \left(\frac{1+{\delta }_k^{-}}{1-{\delta }_k^{-}}\right)\right],$$(67)

with ${\Delta }_k^{\pm }=\sqrt{{\eta }^2-{({\lambda }_k\pm i)}^2}$, ${\delta }_k^{\pm }=(\pm (\eta -{\lambda }_k)-i)/{\Delta }_k^{\pm }$, η = ±2τβ_r, and λ_k = τ(r −k) ⋅ν. The ± in η corresponds to T_±. Since both η and λ_k depend on τ and |Z₀ | (through β_r), the chances of capture in each resonant island change with these two parameters.

In Fig. 6, we show the capture probabilities in the presently observed 3/2 spin-orbit resonance for Mercury (e = 0.206 and |Z₀ |= 1.25 × 10⁻⁵). We see there is a very good agreement between the theoretical approximation (Eq. (66)) and the results of numerical simulations (Eqs. (1) and (15)). For τn < 0.1, the probability is constant and equal to about 8.5%, the same value obtained with a linear tidal model (Correia & Laskar 2004). For τn ~ 1, the probability slightly decreases and then increases very fast. For τn > 10, the probability is always 100%. The transition occurs because for τn ~ 1 the tidal torque changes from a regime proportional to the forcing frequency to a regime inversely proportional (Eq. (55)). As a result, for τn > 1, multiple equilibria for the tidal torque are possible (Fig. 5). One of this equilibrium points is placed exactly at ω∕n= 3∕2, so the rotation rate of the planet is naturally driven to the corresponding resonant island and capture is unavoidable. Indeed, for τn > 10, the planet can remain in the 3/2 spin-orbit equilibria even with Z₀ = 0 (Correia et al. 2014).

In Fig. 7, we show the theoretical capture probabilities (Eq. (66)) for the super-synchronous resonance ν_r = n + ν∕2 for the outer planet in the KOI-1599 system (r = [2, 1], Table 2). We adopt |Z₀| = 10⁻⁷, since for larger values of Z₀ the libration widths of the individual resonances overlap, and thus it is no longer possible to capture the rotation in individual resonances (see Fig. 2). As for the case of Mercury, we observe that for τn ≪ 1 the probability stabilizes at a constant value (54%), while for τn ≫ 1 the probability is 100%. However, there are two transitions of regime. The main transition occurs around τn ~ 10² ⇔ τυ ~ 1 because the dominating terms of the orbital forcing are those involving the orbital libration frequency, υ (Table 2). If we only consider the main three terms of the orbital forcing, this transition, for which the capture probability drops to zero, would be the only observed (green dashed line). The other transition occurs around τn ~ 1 and results from terms that involve the precession rate of the orbit, g. These terms have smaller amplitudes than in the case of Mercury because e ≈ 0.02 (Table 1), but they contribute to increase the capture probabilities in the super-synchronous resonance.

In Fig. 7, we additionally show the results from numerical simulations (Eqs. (1) and (15)). Contrary to the case of Mercury (Fig. 6), in this case the agreement between the theoretical and numerical capture probabilities is not so good. Forτn ≪ 10², this is because the theoretical probabilities (Eq. (66)) were estimated using the approximation (59), which assumes that the effect of the nearby resonances can be averaged to zero. However, although isolated resonant islands exist for |Z₀ |= 10⁻⁷, there is already a small chaotic zone around these islands (Fig. 3, bottom). Therefore, in the numerical simulations the rotation is chaotic for some time before it is stabilized in one of these islands (see Sect. 6) and the probabilities are altered. For τn ≫ 10², in addition to the previous effect, there is a significant increase in the permanent tidal deformation ⟨Z_t ⟩ > 10⁻⁷ (Eqs. (33) and (34)) that extends the size of the chaotic zone (Fig. 3, top). As a result, the rotation may remain permanently chaotic or evolve into some totally unexpected spin-orbit resonance (see next section for more details).

Fig. 5 Normalized net torque $\mathcal{T}/\mathcal{K}$ (Eq. (55)) for different τ values and two orbital forcing functions: Mercury, for which $\left|x_k\right|={\chi }_k^{-3,2}(e)$ and e = 0.206 (left column); the outer planet in the KOI-1599 system, for which the |xk| are given in Table 2 (right column).

Fig. 6 Capture probabilities in the presently observed 3/2 spin-orbit resonance for Mercury (e = 0.206 and |Z0 | = 1.25 × 10−5). The solid line gives the theoretical approximation given by expression (66), while the dots give the results of numerical simulations with Eqs. (1) and (15). We run 180 initial conditions with θ0 differing by one degree.

Fig. 7 Capture probabilities for the super-synchronous resonance νr = n + ν∕2 for the outer planet in the KOI-1599 system with |Z0| = 10−7. The solid line gives the theoretical approximation given by expression (66) with all terms from Table 2, while the dashed line is obtained with only the first three dominating terms. The dots give the results of numerical simulations with Eqs. (1) and (15). We run 180 initial conditions with θ0 differing by one degree.

6 Complete evolution

The spin evolution of a planet is governed by Eq. (1). For a preliminary analysis, this equation can be decoupled in the shape of the planet Z (Eq. (3)) and in the orbital forcing function x (Eq. (4)). The orbital forcing determines the equilibrium states for the spin, and for a given shape we can predict the possible dynamical regimes. However, when we take into account tidal deformation, the shape itself depends on the orbital forcing (Eq. (20)). As a result, the exact behavior of the spin depends on the tidal model that we adopt.

From the analysis of the orbital forcing, we can anticipate which systems can present a rich spin dynamics, such as the KOI-1599 system. However, a deep understanding of all possible scenarios can only be achieved by performing numerical simulations of the complete equations of motion, namely the orbital equations of an N-body system, which give the solution for x(t), together with Eq. (1) for the spin and Eq. (15) for the deformation (tidal model). At the end of the previous section, we already did that for estimating the numerical probabilities of capture in a specific spin-orbit resonance (Figs. 6 and 7).

In this section, we study the complete behavior of the outer planet in the KOI-1599 system (Table 1). This planet has an estimated mass of 4.3 M_⊕, that is, it can be classified as a super-Earth. An upper bound of ten Earth masses is commonly accepted to warrant that the planet is mostly rocky (e.g., Valencia et al. 2007). Therefore, the internal structure and rheology of this planet is likely similar to the terrestrial planets in the solar system. We adopt for the residual deformation Z₀ = 1.0 × 10⁻⁷. For the tidal deformation we use a Maxwell model (Eq. (15)), with k_f = 0.9, τ_e = 0, and C∕mR² = 1∕3. For rocky planets, such as the Earth and Mars, we have a dissipation factor Q = 10 (Dickey et al. 1994) and Q = 80 (Lainey et al. 2007), respectively. We then compute for the Earth τ = 1.6 day, and for Mars τ = 14.7 day. However, in the case of the Earth, the present Q−factor is dominated by the oceans, the Earth’s solid body Q is estimated to be 280 (Ray et al. 2001), which increases the relaxation time by more than one order of magnitude (τ = 46 day). Although these values provide a good estimate for the average present dissipation ratios, they appear to be inconsistent with the observed deformation of the planets. Indeed, in the case of the Earth, the surface post-glacial rebound from the last glaciation about 10⁴ yr ago is still going on, suggesting that the Earth’s mantle relaxation time is something like τ = 4400 yr (Turcotte & Schubert 2002). Since τ can range from a few minutes up to thousands of years, we adopt − 6 ≤ log₁₀τ_yr ≤ 4, which is equivalent to − 4 ≤ log₁₀τn ≤ 6.

In our numerical simulations we use the initial conditions from Table 1. The initial rotation period is set at 19.5 day, which corresponds to $\dot{\theta }/n\approx 1.05$. Since the libration frequency of the outer planet is υ∕n ≈ 0.0066 (Table 2), the initial rotation rate is completely outside the resonant area, which is bounded by $\dot{\theta }/n\lesssim 1.01$ (Fig. 3). For any τ value, the rotation rate of the planet decreases owing to tidal effects until it approaches the region where multiple spin-orbit resonances are present. The shape of the tidal torque depends on the τ value (Fig. 5), so at this stage different scenarios can occur. In addition, capture in resonance is a stochastic process. Therefore, for each τ value we have run 180 simulations with θ₀ differing by one degree. In Table 3 we list the final distribution of the spin in the different resonances for each τ value.

For τn ≪ 10², the spin is in the “low frequency regime” (τυ ≪ 1), which is usually known as the viscous or linear regime (e.g., Singer 1968; Mignard 1979). Because the relaxation time is small, the deformation of the planet is very close to the equilibrium deformation, Z ≈Z₀ + Z_e (Eq. (15)). As a consequence, ⟨Z⟩≈ Z₀, that is, the behavior of the spin follows that of a rigid body (Sect. 3). We can thus predict the dynamics from the frequency map analysis shown in Fig. 3 (bottom). We see that the rotation can be trapped in multiple spin-orbit resonances. Close to the synchronous resonance, there is a small chaotic zone where the spin can wander for some time before it finds a path into a stable resonant island. In Fig. 8 we show three examples of capture in different spin-orbit resonances. From Table 3 we conclude that the most likely resonance is the super-synchronous resonance with ω =n + υ∕2. In the linear regime, the value of τ only impacts the damping timescale, so the statistics do not change much with this parameter.

For τn ~ 10², the spin is in a “frequency transition regime” (τυ ~ 1). In this case, the deformation of the planet is still ⟨Z⟩≈ Z₀, but the tidal torque has an important modification. In the low frequency regime it decreases linearly, allowing capture in all spin-orbit resonances. In the transition regime, however, the torque drives the rotation away from all spin-orbit equilibria except for $|\dot{\theta }/n-1|<0.002$, where there is an inversion (Fig. 5). As a result, only the synchronous resonance can be reached (Table 3). The impact of this inversion in the tidal torque is also pretty clear in Fig. 7, the capture probability in the super-synchronous resonance sharply decreases for τn ~ 10².

For τn ≫ 10², the spin is in the “high frequency regime” (τυ ≫ 1). In this regime, the relaxation time is large and the planet is thus able to retain some tidal deformation ⟨Z_t ⟩≠0 (Sect. 4). For the outer planet in the KOI-1599 system we have $\mathcal{A}=2.7\times {10}^{-6}\gg Z_0$, so for the main spin-orbit resonances the tidal deformation dominates over the permanent deformation. As a result, the phase space of the spin is modified with respect to the phase space drawn for Z₀ = 10⁻⁷, it resembles more to that for Z₀ ~ 10⁻⁶ (Fig. 3, top). We therefore expect chaotic motion for the spin in this regime. However, in our simulations we verify that, although chaotic rotation is indeed a possibility for some τ values, stable rotation is also possible in some asynchronous spin-orbit resonances (Table 3). One reason for this is that in this regime the tidal torque facilitates the capture in all these resonances (Fig. 5). The other reason is because the tidal deformation varies from resonance to resonance (Eq. (30)). For the synchronous resonance |⟨Z_t ⟩|~ 10⁻⁶, so even if the rotation is somehow captured there, it becomes chaotic after some time. For the super-synchronous resonance |⟨Z_t ⟩|≈ 3.6 × 10⁻⁷ (Eq. (53)), so it is not impossible to remain trapped there, provided that the spin finds a way in.

In Fig. 9 we show an example for τn = 10⁶, which clarifies the behavior of the spin in the high frequency regime. In this figure we show the evolution of the rotation rate (top) together with the evolution of the shape (bottom). Before the spin reaches the resonant area, we have ⟨Z⟩ = Z₀ = 10⁻⁷, so in principle the rotation can be trapped in a stable spin-orbit resonance (Fig. 3, bottom). As soon as the spin enters the resonant area, the deformation increases and can experience values |Z| > 10⁻⁶. Then, the rotation becomes chaotic (Fig. 3, top). However, just outside the chaotic zone there exists the sub-synchronous resonance with ω = n − υ. Somehow, after some time in the chaotic region, the spin finds a path into this stable region. Once there, the deformation of the planet decreases again to a value |Z|~ 10⁻⁷. Thus, the sub-synchronous resonance is no longer connected with the chaotic region and the spin can remain there forever (Fig. 3, bottom). In this example, it is remarkable that although the spin starts with an initial rotation rate faster than the synchronous rotation, it finally ends up with a rotation slower than that.

Fig. 8 Some examples of the spin evolution of the outer planet in the KOI-1599 system for τn = 10−4 (low frequency regime). The blue lines give the position of the main spin-orbit resonances n, n ± υ∕2, and n ± υ.

Fig. 9 Spin evolution of the outer planet in the KOI-1599 system (top) and the corresponding deformation |Z| (bottom) for τn = 106 (high frequency regime). The blue lines give the position of the main spin-orbit resonances n, n ± υ∕2, and n ± υ.

7 Discussion

In this paper we have revisited the spin-orbit coupling of close-in planets. These bodies are assumed to be tidally evolved, and thus can be trapped in some equilibrium states. We derived our results using a very general approach, that is suitable to apply to multi-planetsystems in various orbital configurations (e.g., very eccentric orbits or involved in orbital resonances). We introduced a simple way to probe the spin dynamics from the orbital perturbations (Eq. (12)), a new method for computing the forced librations (Eq. (35)) and the permanent tidal deformation (Eq. (53)), and a general expression for the tidaltorque (Eq. (55)) and for the capture probabilities in resonance (Eq. (67)).

We considered that the spin of the planet is only subject to tidal torques from the parent star (Eq. (1)). Tidal torques between planets are also present, but their contribution is very small for the main spin-orbit resonances (see Table II in Goldreich & Peale 1966), and accounts for less than a few percent of the total tidal dissipation (Hay & Matsuyama 2019). In addition, the positions of the spin-orbit resonances do not change because they depend on the fundamental frequencies of the system; the additional torques only modify the libration width of the resonant islands. However, if we consider a system with two or morestars, our model should be corrected to include the torques from all bodies on the spin (Correia et al. 2015; Leleu 2017).

In our model we assumed for simplicity that the spin axis of the planet is normal to the orbit, since this is the final outcome of tidal evolution (Hut 1980). For multi-planet systems, the obliquity does not evolve exactly to zero degrees, it is rather captured in an equilibrium configuration known as Cassini state (Colombo 1966; Peale 1969). Although high obliquity states exist, they are usually unstable under strong tidal interactions. Therefore, the only stable equilibrium state has in general very small obliquity (less than 1°), for which our model applies. We note, however, that this may not be the case for planetary systems showing large mutual inclinations (Correia 2015).

We also assumed coplanar orbits. This restriction was used to simplify the orbital forcing function x (Eq.(6)), since this function has only two dimensions in the planar case. If the system is not coplanar we can still use our model because the forcing function captures the effect of the inclination. It is nevertheless convenient to align the inertial plane where we compute the forcing function with the mean orbital plane of the planet whose rotation we want to study. This alignment is also important to ensure that the obliquity remains close to zero.

In our model we only considered gravitational tidal torques. For Earth-size or larger planets, we can expect that a dense atmosphere is present. As a result, thermal atmospheric tides may give rise to additional equilibrium states (Correia & Laskar 2001). However, it has been shown that gravitational tides dominate for close-in planets with a ≲ 0.3 au (Correia et al. 2008; Leconte et al. 2015). Thus, as long as we stay in this regime, the presence of an atmosphere should not modify the global picture drawn in this paper much.

We have shown that the spin dynamics of Earth-size planets in compact systems can be very rich. In particular, we have shown that synchronous rotation is not the only possibility; the rotation is allowed to evolve into some unexpected spin-orbit resonances or even present chaotic behavior. In the high frequency regime, we have shown that the tidal torque has equilibrium points at exactly the same position as the spin-orbit resonances, which enhances the chances of capture in these equilibria. Future studies can take into account all the neglected contributions in our model, which may increase the diversity of spin states for these planets even more.

Acknowledgements

We thank G. Boué, A. Leleu, and P. Robutel for discussions. A.C. acknowledges support by CFisUC strategic project (UID/FIS/04564/2019), ENGAGE SKA (POCI-01-0145-FEDER-022217), and PHOBOS (POCI-01-0145-FEDER-029932), funded by COMPETE 2020 and FCT, Portugal. J.B.D. acknowledges support by the Swiss National Science Foundation (SNSF). This work has, in part, been carried out within the framework of the National Centre for Competence in Research PlanetS supported by SNSF.

References

All Tables

All Figures

	Fig. 1 Fast Fourier transform of the orbital forcing x (Eq. (7)) for the two-body problem with eccentricity e = 0.165. The horizontal lines give the amplitudes βk for |Z0 | = 1.5 × 10−5, which corresponds to Mercury’s value (Genova et al. 2013), and |Z0| = 2.9 × 10−2, which corresponds to Hyperion’s value (Thomas et al. 2007).
In the text

	Fig. 2 Fast Fourier transform of the orbital forcing x (Eq. (7)) for the outer planet in the KOI-1599 system. The horizontal lines give the amplitudes βk for |Z0 |= 1.5 × 10−6 and |Z0 | = 1.0 × 10−7.
In the text

Fig. 3 Frequency map analysis of KOI-1599 outer planet rotation angle assuming a rigid body with |Z0 | = 1.5 × 10−6 (top) or |Z0| = 1.0 × 10−7 (bottom). Each dot corresponds to a grid with 181 × 201 different initial conditions of γ(0) and $\dot{\gamma }(0)$, respectively, with γ(t) = θ(t) − λ(0). The color index gives the derivative $\partial \eta /\partial \dot{\gamma }(0)$, where η is the main frequency of γ(t). The blue dots correspond to resonant motion, green dots to nonresonant regular motion, and red dots to chaotic motion.

In the text

	Fig. 4 Forced libration of the rotation rate for the super-synchronous resonance (νr = n+ υ∕2) of the outer planet in the KOI-1599 system. We show the evolution of $(\dot{\theta }-{\nu }_r)/n$ adopting τe = 0 and τn = 104.
In the text

	Fig. 5 Normalized net torque $\mathcal{T}/\mathcal{K}$ (Eq. (55)) for different τ values and two orbital forcing functions: Mercury, for which $\left|x_k\right|={\chi }_k^{-3,2}(e)$ and e = 0.206 (left column); the outer planet in the KOI-1599 system, for which the |xk| are given in Table 2 (right column).
In the text

	Fig. 6 Capture probabilities in the presently observed 3/2 spin-orbit resonance for Mercury (e = 0.206 and |Z0 | = 1.25 × 10−5). The solid line gives the theoretical approximation given by expression (66), while the dots give the results of numerical simulations with Eqs. (1) and (15). We run 180 initial conditions with θ0 differing by one degree.
In the text

Fig. 7 Capture probabilities for the super-synchronous resonance νr = n + ν∕2 for the outer planet in the KOI-1599 system with |Z0| = 10−7. The solid line gives the theoretical approximation given by expression (66) with all terms from Table 2, while the dashed line is obtained with only the first three dominating terms. The dots give the results of numerical simulations with Eqs. (1) and (15). We run 180 initial conditions with θ0 differing by one degree.

In the text

	Fig. 8 Some examples of the spin evolution of the outer planet in the KOI-1599 system for τn = 10−4 (low frequency regime). The blue lines give the position of the main spin-orbit resonances n, n ± υ∕2, and n ± υ.
In the text

	Fig. 9 Spin evolution of the outer planet in the KOI-1599 system (top) and the corresponding deformation |Z| (bottom) for τn = 106 (high frequency regime). The blue lines give the position of the main spin-orbit resonances n, n ± υ∕2, and n ± υ.
In the text